# Diffuse emission from black hole remnants

Sina Kazemian<EMAIL_ADDRESS>Dept. of Physics & Astronomy, Western
University, N6A 3K7, London ON, Canada  Mateo Pascual<EMAIL_ADDRESS>Dept. of
Physics & Astronomy, Western University, N6A 3K7, London ON, Canada  Carlo
Rovelli<EMAIL_ADDRESS>Dept. of Philosophy and Rotman Institute, Western
University, N6A 3K7, London ON, Canada  Aix Marseille University, Université
de Toulon, CNRS, CPT, 13288 Marseille, France  Perimeter Institute, 31
Caroline Street North, N2L 2Y5 Waterloo ON, Canada  Francesca Vidotto
<EMAIL_ADDRESS>Dept. of Physics & Astronomy, Western University, N6A 3K7,
London ON, Canada  Dept. of Philosophy and Rotman Institute, Western
University, N6A 3K7, London ON, Canada

(August 29, 2024)

###### Abstract

At the end of its evaporation, a black hole may leave a remnant where a large
amount of information is stored. We argue that the existence of an area gap as
predicted by Loop Quantum Gravity removes a main objection to this scenario.
Remnants should radiate in the low-frequency spectrum. We model this emission
and derive properties of the diffuse radiation emitted by a population of such
objects. We show that the frequency and energy density of this radiation,
which are measurable in principle, suffice to estimate the mass of the parent
holes and the remnant density, if the age of the population is known.

##  I Introduction

The possibility that black holes tunnel into long-living remnants at the end
of their evaporation J.H.MacGibbon1987 ; Giddings1992a ; haggard_black_2015 ;
haggard_quantum-gravity_2015 ; bianchi_white_2018 ; Rovelli2018f  has recently
been receiving renewed attention Adler2001 ; Nozari2008 ; Inomata2021 ;
DiGennaro2021 ; Green2021 ; EslamPanah2020 ; chen_black_2003 ; Chen2004 . Here
we study the emission that we may expect from such remnants.

As a black hole approaches the end of its evaporation, it enters a Planckian
regime, since the curvature of its surroundings reaches the Planck scale. A
number of recent results from non-perturbative quantum gravity as well as from
classical General Relativity have revived the old idea J.H.MacGibbon1987 ;
Giddings1992a  that the end of the evaporation could leave long living Planck
scale remnants. These results are: (i) _Classical_ General Relativity provides
a surprisingly natural model for such remnants: white holes with a small
horizon but very large interior haggard_black_2015 . These are exact solutions
of Einstein’s equations. (ii) _Classical_ General Relativity allows for the
existence of spacetimes where such white holes are in the future of the parent
black hole, after a quantum tunnelling transition localised in space and time
haggard_black_2015  (the black to white tunnelling transition is itself an old
idea narlikar_high_1974 ; hajicek_singularity_2001 ; ambrus_quantum_2005 ;
olmedo_black_2017 ).  (iii) Non-perturbative calculations in Loop Quantum
Gravity (LQG) show that the tunnelling transition is permitted, and
increasingly probable towards the end of the evaporation christodoulou_planck-
star_2016 ; christodoulou_characteristic_2018 ; DAmbrosio2021 ; soltani2021 .
(iv) Planck scale white holes can be stabilised against instability by quantum
gravity rovelli_small_2018 .  (v) A number of objections that made the remnant
idea unconvincing a few decades ago banks_are_1992 ; giddings_dynamics_1992 ;
banks_black_1993 ; giddings_constraints_1994 ; banks_lectures_1995  have now
been shown not to apply to this scenario bianchi_white_2018 . We also give
below a further argument removing previous objections to the idea of remnants.

Thanks to these results, long living remnants are again plausible outcomes for
the end of the evaporation. Notice that the alternative idea that black holes
could magically pop out of existence is a scenario not directly predicted or
supported by any quantum gravity theory, and hard to harmonize with the
persistence of the large volume inside evaporating black holes
Christodoulou2015 ; de_lorenzo_improved_2016 .

If a black hole ends up tunneling into a white hole, its horizon is not an
event horizon because its interior is causally connected with future null
infinity. On the other hand, causality prevents the black hole from being an
ergodic system, because energy cannot freely traverse the horizon in both
directions. For a non-ergodic system, the von Neumann entanglement entropy can
be higher than the thermodynamic entropy, because all degrees of freedom can
contribute to the first, but they do not necessarily contribute to the second
Rovelli2017e ; Rovelli2019a . Hence the von Neumann entanglement entropy
across the horizon of a black hole can remain high even when the Bekenstein-
Hawking entropy decreases throughout the evaporation. This is a possibility
which is manifestly distinct from the popular strong version of the holography
‘dogma’ Maldacena2020 , but does not contradict any known physics. In this
scenario, the interior state of the black hole can still be highly entangled
with the emitted Hawking radiation when the black hole reaches Planck mass and
this is sufficient to purify the Hawking radiation. Information has no reason
to start escaping at Page time and remains stored in the (vast) hole’s
interior when the horizon tunnels from black to white.

But eventually information has to come out, before the final dissipation of
the white hole. The minuteness of the white hole horizon’s area and energy
implies that this can only happen slowly preskill_black_1992 . Bringing out a
large amount of information involving only little energy is what gives rise to
the low energy radiation that we model here.

Here we study the characteristics of this radiation and in particular the
diffused background produced by a population of remnants. We use a crude
model, with only a few parameters, that provides a quantitative estimate of
the property of this radiation as a function of the original mass $m$ of the
black hole at its time of formation. We stress the fact that what we give is
not a first-principles derivation of the radiation profile (as Hawking did
with the black hole radiation), but only an estimate of its properties based
on a simple modelization and exploiting basic conservation principles.

Our results confirm previous estimates that black hole remnants take a long
time to radiate ($t\sim m^{4}$ in natural units). We also estimate the
frequency and amplitude of this diffuse background radiation. These determine
the number density of the remnants, as a function of the characteristic mass
and time of formation of the parent black hole population. We also discuss the
quantum field theoretical picture of the remnant-radiation coupling, showing
that a key result in non-perturbative Loop Quantum Gravity resolves a
potential objection to this scenario.

An interesting possibility is that a population of small long-living black
hole remnants produced by primordial black holes Carr1974a , born either in
the early universe or before a big-bounce, could be a component of dark matter
Rovelli2018f ; chen_black_2003 ; Chen2004 ; Green2021 ; Inomata2021 ;
DiGennaro2021 ; EslamPanah2020 . At first sight, Planck size white holes
stabilized by quantum gravity seem natural cold dark matter candidates, since
they interact almost only gravitationaly, behave like a rarefied dust of
micro-gram size grains, and their existence is predicted by current
fundamental theory (General Relativity plus quantum tunnelling) without any
new physics. A preliminary exploration of the version of this scenario with
post big-bang black holes derome_improved_2003  finds difficulties with it,
but not with the bouncing scenario. A white-hole-remnant dark matter component
would be hard to detect directly. Perhaps the background low frequency
emission we describe here could help in this regard. Here, however, we only
study the physics of the emission by a population of these objects in flat
space: we do not delve into any question concerning its cosmological evolution
or concrete possibility of observation.

##  II The model

The hypothesis of the model we study is the following. A black hole of initial
mass $m$ evaporates via Hawking evaporation, leaving behind a remnant of
Planckian mass which contains an amount of information that is sufficient to
purify its Hawking radiation, namely of order

$S\sim\frac{A}{4}=4\pi m^{2}$ (1)

in natural units $\hbar=G=c=k=1$. Here $A$ and $m$ are the area of the horizon
and the mass of the hole at its formation, not to be confused with area and
mass at the time of tunnelling, which are presumably approximately Planckian.
Our assumption is that the information associated to $S$ is later emitted in
the form of radiation.

Let $\tau_{W}$ be the lifetime of the remnant, namely the time lapsed from the
black-to-white transition of the hole to the moment the white hole has
completely dissipated. By causality, at this final time the radiation is
spread non-uniformly over a sphere of radius $L=\tau$ (remember that we have
$c=1$). The radiation is emitted radially, hence its propagation direction is
everywhere in a single direction. Developing an idea hinted at in
preskill_black_1992 , we model the radiation along a single radius of this
sphere as a uniform one-dimensional gas of photons in equilibrium. This is a
crude approximation, but presumably sufficient to give us orders of magnitude.

The total energy $E$ available for this gas of photons is only that of the
mass of the remnant, which is of the order of the Planck mass, namely

$E\sim 1$ (2)

in natural units. Its total entropy, needed to purify the Hawking radiation,
is (1), distributed over the length $L$.

A standard derivation, recalled in the appendix for completeness, shows that
the entropy $S$ and energy $E$ of a one dimensional photon gas of temperature
$T$ in a space of length $L$ are schwabl_statistical_2006 ;
skobelev_entropy_2013

$S=\frac{2\pi}{3}LT,\ \ \ \ E=\frac{1}{6}LT^{2}.$ (3)

Inverting these two relations we find

$L=\frac{3\text{S}^{2}}{8\pi^{2}\text{E}}=6m^{4},\ \ \ \
T=\frac{4\pi\text{E}}{\text{S}}=\frac{1}{m^{2}}.$ (4)

From the definition of $L$ we have that the lifetime of the white hole is
therefore

$\tau_{W}\sim 6m^{4}$ (5)

which matches previous estimates of the time needed to release the information
contained inside the remnant preskill_black_1992 ; Marolf2017 ;
rovelli_small_2018 ; bianchi_white_2018 . The estimate indicates that the
temperature of the white hole is much lower than the Hawking temperature of
the parent black hole which evolves as $1/m(t)$ ($m(t)$ is the mass of the
black hole at time $t$, as opposed with $m$ being the mass of the black hole
when it first forms), thus increasing during the evaporation. Notice that the
model corrects the infinite explosion predicted by Hawking111The title of his
1974 Nature paper HAWKING1974  was _‘Black hole explosion?’_ : the black hole
temperature grows only up to the Planckian value, then drops abruptly to the
low value (4).

At this point, we can use the crude model we’ve constructed describing a gas
of photons in equilibrium to furthermore estimate the total number of photons
emitted by the remnant. At equilibrium, the peak frequency of the (Planckian)
distribution of photons is

$\nu=\alpha\;T=\frac{\alpha}{m^{2}}$ (6)

where $\alpha\sim 2.82$. The derivation of the constant $\alpha$ is reported
in the appendix for completeness.222If the frequency is low enough, then free-
free absorption makes interstellar gas opaque King1980 . The Milky Way shows a
spectral turnover at 3 MHz attributed to free-free opacity Lacki2012 . We
thank an anonymous referee for pointing this out..

In natural units, the relation between the energy $\epsilon$ of a single
photon and its frequency $\nu$ is of course $\epsilon=\nu$, hence we can
derive the total number of photons emitted by the remnant of a black hole of
initial mass $m$ to be

$N_{\gamma}=\frac{E}{\epsilon}=\frac{m^{2}}{\alpha}.$ (7)

Since the emission is radial, the gas of photons at the end of the dissipation
process would span a ball of volume $V=\frac{4}{3}\pi L^{3}$. At this time,
the average energy density in this volume is therefore

$\rho_{o}=\frac{E}{\frac{4}{3}\pi L^{3}}=\frac{1}{288\,\pi}m^{-12}$ (8)

and the average photon density

$n_{\gamma}=\frac{N_{\gamma}}{\frac{4}{3}\pi
L^{3}}=\frac{1}{288\,\pi\alpha}m^{-10}.$ (9)

Let us now consider a population of remnants uniformly distributed in space
with number density $\Omega$, resulting from the evaporation of black holes
which all formed at the same time and with the same mass. Neither the energy
emitted by a single remnant nor the density of emitted photons are uniformly
distributed in space, but the ensemble of all photons emitted by a uniform
population of remnants is indeed a uniformly distributed bath of radiation on
cosmic scales.

The total energy emitted at the end of the process must be equal to their
total initial mass. Since they have unit mass (in natural units), this is
equal to their total number. Hence, the energy density of the radiation
$\rho_{tot}$ at the end of the process is equal to the initial number density
of remnants,

$\rho_{tot}=\Omega$ (10)

and the total photon number density (due to the entire population) is

$n=\Omega\,N_{\gamma}=\Omega\,\frac{m^{2}}{\alpha}\,.$ (11)

In a cosmological context, the same result holds, but $\Omega$ and equation
(11) express comoving number densities instead.

###  II.1  Linear emission

Consider a population of black holes formed at a time $t=0$, with mass $m$ and
uniformly distributed in space. Assume that they all evaporate around time
$\tau_{B}\sim m^{3}$ as predicted by Hawking radiation theory, and survive as
white hole remnants for a time $\tau_{W}$ as in (5). Between times $\tau_{B}$
and $\tau_{B}+\tau_{W}$, they emit a steady radiation as described above.
Assuming $m\gg 1$, we can approximate $\tau_{B}+\tau_{W}\sim\tau_{W}$. What is
the radiation observed by an observer at time $t$?

A fist approximate answer can be given assuming that the radiation is emitted
steadily in time. For $t<\tau_{B}$, there is none. For
$\tau_{B}<t<\tau_{B}+\tau_{W}$, the observer will receive only the radiation
emitted by the remnants within a distance $r<(t-\tau_{B})$ because radiation
emitted by more distant remnants has not had enough time to reach the
observer. Radiation emitted at a distance $r$ is diluted by distance by a
factor $1/r^{2}$, but the number of emitters at this distance is proportional
to $r^{2}$, hence the radiation received is proportional to $r<(t-\tau_{B})$.
(Recall that we are working on flat space. The result can be extended to a
cosmological context, but we do not do so here.) For the same reason, when
$t>\tau_{B}+\tau_{W}$ and the remnants have completed their emission, the
radiation received remains constant in time. That is to say, the radiation
density changes in time as

$\rho(t)\left\\{\begin{array}[]{lll}=0&{\rm for}&t<m^{3},\\\
=\big{(}\frac{t-\tau_{B}}{\tau_{W}-\tau_{B}}\big{)}\Omega&{\rm
for}&m^{3}<t<6m^{4},\\\ =\Omega&{\rm for}&t>6m^{4}.\end{array}\right.$ (12)

In other words, the process is a steady (linear in time) transformation of
dust into radiation, on a $m^{4}$ timescale. As we see below, however, quantum
mechanics corrects this result.

Figure 1: Background white hole radiation as a function of time. The solid
black line represents a classical linear emission while the dashed red line
represents a quantum emission.

###  II.2  Quantum description of the emission

The energy of a white hole is related to the area of its horizon. A
continuous, steady emission of energy from each white hole, as the one
described above, would imply a continuous decrease of the white hole horizon
area, even below the Planck area. According to LQG, however, any physical area
is quantized, with the minimum non-zero eigenvalue (the ‘area gap’) being of
the order of the Planck area $A_{Pl}$ Rovelli1993c ; Smolin1994 ; Rovelli1994a
. Therefore a remnant with near-Planckian mass and area cannot emit an amount
of energy smaller than an energy of the order of the Planck energy. That is,
it can only make a single quantum leap into radiation. This is analogous to
conventional nuclear radioactivity, where the steady emission of a macroscopic
sample of radioactive material is realised by individual quantized emissions
of its constituent atoms, governed by a probability distribution
martin_nuclear_2019 .

More precisely, an area gap of the order of the Planck area implies that the
energy of the lowest non-vanishing energy states of the remnants is Planckian.
Therefore, in first-order perturbative formulation, the only allowed
transition with the emission of radiation (which necessarily has energy) is
emission of the entire Planck energy of the remnant.

Let us see, in the language of quantum field theory, what the corresponding
vertex describing the transition could be. The essential point is that (black
and) white holes have many internal degrees of freedom that reflect their
internal structure. A white hole parented by an old black hole evaporated from
a initial mass $m$ has an interior capable of holding information compatible
with (1) even if the area of its horizon is small. A vertex coupling such
remnant to a single or a few photons is therefore forbidden by conservation of
information (unitarity), because a few photons do not have enough degrees of
freedom to match the large number of quantum numbers describing the white hole
interior. Just a few photons cannot carry the entire information that can be
stored in the remnant. Hence the only possible transition is a transition
$remnant\to\gamma_{1}...\gamma_{n}$ to a large number of low energy photons:

This conclusion is interesting in view of an old objection to the remnant
scenario, because of which this scenario was abandoned in the Nineties
preskill_black_1992 . The objection was that the large number of remnant
internal states would entail that their production in particle physics
experiments would be too easy, hence expected be observed already. Here we see
clearly why that conclusion was too quick. The effective vertex responsible
for a remnant production would actually have to be

in order to create a long-living, Planck-size remnant. The number of photons
emitted by a single remnant is given in (7). If the number of photons coming
together is small, they can be highly energetic, but the remnant produced
would correspond to a remnant whose parent is a black hole of Planckian size,
which is short-lived. The process would not be distinguishable by the standard
possibility of collapse predicted by conventional quantum gravity. To produce
an actual long living remnant, on the other hand, we need $m$ to be large, and
hence we would need to focus a large number of low energy photons.

For instance to produce a remnant similar to the one left over from a
primordial black hole formed at reheating (see below) the number of low energy
photons to focus would be staggering:

$5\cdot 10^{38}<N_{\gamma}<5\cdot 10^{48}.$ (13)

Creating such a remnant in the lab is clearly unlikely due to the huge number
of photons required for the process to happen. Therefore not being able to
create remnants by the present experimental settings is not a reason to reject
the theory of black holes turning into Planck size white hole remnants at the
end of their life time cycle.

Let us see how this quantum effect corrects the expected emission of a
population of remnants. If there is a single decay into multiple photons for
each remnant, and the probability of transition is constant in time, then a
_population_ of remnants behaves like a radioactive material: the number of
remnants, and therefore the total energy density $\rho_{rem}$ of the
population of remnants decays exponentially, starting at $t=\tau_{B}$, as

$\rho_{rem}(t)=\Omega\ e^{-\lambda(t-\tau_{B})}.$ (14)

If the lifetime of the white hole is of order $\tau_{W}$ (by Bohr’s
correspondence principle), we expect the decay constant to be

$\lambda\sim(\tau_{W}-\tau_{B})^{-1}.$ (15)

The energy density of the radiation as a function of time is therefore

$\rho_{rad}(t)\left\\{\begin{array}[]{lll}=0&{\rm for}&t<m^{3},\\\
=\Omega\bigg{(}1-e^{-\frac{t-\tau_{B}}{\tau_{W}-\tau_{B}}}\bigg{)}&{\rm
for}&t>m^{3}.\end{array}\right.$ (16)

In Figure 1 we have plotted the energy density of the linear emission (12) as
a black solid line, and the quantum emission (16) as a dashed red line. The
two converge in the two limits $\tau_{B}<t\ll\tau_{W}$ and $t\gg\tau_{W}$.

The energy density $\rho_{rad}$ of the radiation emitted by a population of
remnants with current energy density $\rho_{rem}$, generated by parent black
holes that all formed a time $t$ in the past with mass $m$, is then obtained
by combining equations (14) and (16), which gives

$\rho_{rad}(t,m)=\rho_{rem}\left[\exp\\!\left(\frac{1-tm^{-3}}{1-6m}\right)-1\right].$
(17)

Using (6) we can write the mass in terms of the peak frequency of the
radiation, and give the energy density in radiation as

$\rho_{rad}(t,\nu)=\rho_{rem}\left[\exp\\!\left(\frac{1-t\,(\nu/\alpha)^{3/2}}{1-6\sqrt{\alpha/\nu}}\right)-1\right].$
(18)

##  III Dimensionful estimates for primordial holes

In the cosmological standard model, primordial black holes may have formed at
reheating. To get a sense of the characteristic nature of the diffuse
radiation remnants may emit, we estimate its parameter in this simplest case.
We use here a rough model, which neglects the effect of cosmological
expansion.

For such a population of primordial black holes, we can approximate $t$ in the
above formulas to be the Hubble time $t_{H}$. Notice that this approximation
allows us to deduce the density of an otherwise dark population of remnants
just from the observation of the emitted radiation. (In other cosmological
scenarios, in particular in bouncing models ashtekar_loop_2015 ;
brandenberger_bouncing_2017 ; ijjas_bouncing_2018 , $t$ can be larger.)

Restoring physical units, and denoting the Plank mass, energy, frequency and
time as $m_{Pl}$, $E_{Pl}$, $\nu_{Pl}$, $t_{Pl}$ we find that a population of
primordial black holes of formation mass $m$ and comoving number density
$\Omega$ gives rise to remnants emitting diffuse radiation with density

$\rho_{rad}(m)=\Omega
E_{Pl}\bigg{(}1-e^{\frac{1-t_{H}/t_{Pl}\,(m/m_{Pl})^{-3}}{1-6m/m_{Pl}}}\bigg{)}$
(19)

and frequency given by (6), namely

$\nu=\left(\frac{m}{m_{Pl}}\right)^{\\!\\!-2}\nu_{Pl}.$ (20)

If we are in the era where this radiation forms, we must have
$\tau_{B}<t_{H}<\tau_{W}$, which is to say

$(m/m_{Pl})^{3}<t_{H}/t_{Pl}<6(m/m_{Pl})^{4}.$ (21)

Since $t_{H}\sim 10^{61}t_{Pl}$ this gives the approximate mass range

$10^{15}m_{Pl}<m<10^{20}m_{Pl}.$ (22)

The model is thus entirely determined by a single parameter or order of unity,
that can be taken to be

$x=\log_{10}(m/m_{Pl})\in(15,20].$ (23)

And the relevant quantities are

$\displaystyle m$ $\displaystyle=$ $\displaystyle 10^{x-5}{\rm g},$ (24)
$\displaystyle\nu$ $\displaystyle=$ $\displaystyle 5\cdot 10^{-2x+42}{\rm Hz}$
(25) $\displaystyle\rho_{rad}(x)$ $\displaystyle=$
$\displaystyle\rho_{rem}\left[\exp\\!\left(\frac{10^{61}-10^{3x}}{10^{4x}-10^{3x}}\right)-1\right]$
(26)

This is a mass range

$10^{10}{\rm g}<m<10^{15}{\rm g},$ (27)

and a frequency range

$5\cdot 10^{12}{\rm Hz}>\nu>5\cdot 10^{2}{\rm Hz}.$ (28)

The ratio of the radiation density to the total density of remnants and
radiation, as a function of $x$, is shown in Figure 2. Notice that remnants
originating from parent black holes in the mass range of $10^{10}g<m<2\cdot
10^{10}g$ have emitted most of their energy in radiation as of today, while
remnants originating from more massive black holes, $4\cdot 10^{10}<m<10^{15}$
have emitted close to zero. This is because of the long lifetime of white hole
remnants which is in the order of $\tau_{W}=6m^{4}$.

Figure 2: Ratio of radiation to total mass as function of the single parameter
of the model $x$.

##  IV Conclusions

We have shown that a quantum gravity effect, namely the quantization of the
area, removes a main objection to the possibility that the end of a black
hole’s evaporation leaves behind a Plank scale remnant. The area gap gives a
wide forbidden energy region, because of which the remnant can only couple to
a large number of low energy photons. This makes remnant production hard to
achieve in a lab.

We have estimated that a diffuse radiation bath peaked at frequency $\nu$ with
energy density $\rho_{rad}$ can account for a population of white hole
remnants descending from black holes formed a time $t$ in the past, with mass

$m=10^{x}m_{Pl}$ (29)

where

$x=-\frac{1}{2}\ \log_{10}\frac{\nu}{\nu_{Pl}},$ (30)

and density

$\rho_{rem}=\rho_{rad}\left[\exp\\!\left(\frac{1-t\,(\nu/\alpha)^{3/2}}{1-6\sqrt{\alpha/\nu}}\right)-1\right]^{-1}.$
(31)

For black holes formed a Hubble time in the past, this approximately becomes

$\rho_{rad}(t,\nu)=\rho_{rem}\left[\exp\\!\left(\frac{1-10^{61-3x}}{1-6\cdot
10^{-x}}\right)-1\right],$ (32)

where $x$ can be measured directly from the frequency of the diffused
radiation. These results are in flat spacetime, and can be easily generalised
to an expanding cosmological setting.

To evaluate the hypothesis that remnants might form a component of dark
matter, these results must be corrected by taking into account the
cosmological evolution. For black holes formed after the big bang in the
primordial universe, there are strong constraints barrau_closer_2021 . In
other cosmological scenarios such as big bounce or matter bounce scenarios
white hole remnants might account for a significant portion of dark matter.

$***$

Acknowledgments We thank Emily Adlam, Patrick Fraser, Pietropaolo Frisoni,
Yichen Luo and Farshid Soltani for participating in many discussions about
this work. This work was supported by the QISS JFT grant 61466. Research in
FV’s group at Western University is supported by the Canada Research Chairs
Program and by the Natural Science and Engineering Council of Canada (NSERC)
through the Discovery Grant ‘Loop Quantum Gravity: from Computation to
Phenomenology’. We acknowledge the Anishinaabek, Haudenosaunee, Lūnaapéewak
and Attawandaron peoples, on whose traditional lands Western University is
located.

##  Appendix A

Here we derive of the entropy and the energy of a one-dimensional photon gas
emitted by a remnant. We also derive the frequency of the photon gas as a
function of temperature.

The Entropy of the photon gas can be written as skobelev_entropy_2013 :

$S=\frac{k_{B}^{2}}{\hbar}\delta^{n}\frac{4(n+1)(n-1)!\zeta(n+1)}{\Gamma(\frac{n}{2})},$
(33)

where we have:

$\delta=\frac{LT}{2\sqrt{\pi}},$ (34)

where T is the temperature and L is the distance traveled by the photon gas.
The Gamma function is equal to:

$\Gamma(n)=(n-1)!,$ (35)

and the Riemann series function is equal to:

$\zeta(n)=\sum_{k=1}^{\inf}\frac{1}{k^{n}}.$ (36)

We assume a one-dimensional photon gas radiating from the remnant and put
$n=1$. The entropy of a one-dimensional photon gas is equal to:

$S=\delta\frac{8\zeta(2)}{\Gamma(\frac{1}{2})}=\frac{2\pi}{3}LT.$ (37)

The energy of a radiating photon gas can be written as:

$E=Tt^{n}\frac{4n!\zeta(n+1)}{\Gamma(\frac{n}{2})},$ (38)

where for n = 1 the energy is equal to:

$E=\frac{1}{6}LT^{2}.$ (39)

To calculate the frequency of the photons emitted by the remnant we assume the
system under study to behave like a black body radiation and write a Planck
distribution for it. The energy density is then equal to
schwabl_statistical_2006 :

$u(\nu,T)=\frac{2h\nu^{3}}{c^{2}}\frac{1}{e^{h\nu/k_{B}T}-1}.$ (40)

By taking the derivative of the energy density with respect to frequency and
putting it equal to zero we find the peak frequency as following:

$\frac{h\nu}{k_{B}T}\frac{e^{h\nu/k_{B}T}}{e^{h\nu/k_{B}T}-1}-3=0.$ (41)

By changing the variable as $h\nu/k_{B}T=\alpha$ we have:

$(\alpha-3)e^{\alpha}+3=0,$ (42)

where we get:

$\alpha\approx 2.82.$ (43)

The peak frequency is therefore equal to:

$\nu_{peak}=\frac{k_{B}}{h}\alpha T.$ (44)

where the maximum number of emissions happen at this frequency.

## References

  * (1) J. H. MacGibbon, “Can Planck-mass relics of evaporating black holes close the Universe?,” Nature 329 (1987), no. 24, 308.
  * (2) S. Giddings, “Black holes and massive remnants,” Physical Review D 46 (1992), no. 4, 1347–1352, 9203059.
  * (3) H. M. Haggard and C. Rovelli, “Black to white hole tunneling: An exact classical solution,” International Journal of Modern Physics A 30 (Oct., 2015) 1545015. Publisher: World Scientific Publishing Co.
  * (4) H. M. Haggard and C. Rovelli, “Quantum-gravity effects outside the horizon spark black to white hole tunneling,” Physical Review D 92 (Nov., 2015) 104020. Publisher: American Physical Society.
  * (5) E. Bianchi, M. Christodoulou, F. D’Ambrosio, H. M. Haggard, and C. Rovelli, “White Holes as Remnants: A Surprising Scenario for the End of a Black Hole,” Classical and Quantum Gravity 35 (Nov., 2018) 225003\. arXiv:1802.04264 [gr-qc].
  * (6) C. Rovelli and F. Vidotto, “Small black/white hole stability and dark matter,” Universe 4 (nov, 2018) 127, 1805.03872.
  * (7) R. J. Adler, P. Chen, and D. I. Santiago, “The generalized uncertainty principle and black hole remnants,” General Relativity and Gravitation 33 (jun, 2001) 2101–2108, 0106080.
  * (8) K. Nozari and S. H. Mehdipour, “Gravitational Uncertainty and Black Hole Remnants,” Modern Physics Letters A 20 (2008), no. 38, 2937–2948, 0809.3144.
  * (9) K. Inomata, E. McDonough, and W. Hu, “Primordial black holes arise when the inflaton falls,” Physical Review D 104 (apr, 2021) 2104.03972.
  * (10) S. Di Gennaro and Y. C. Ong, “Feasibility of primordial black hole Remnants as dark matter in view of Hawking radiation recoil,” Journal of Cosmology and Astroparticle Physics 2021 (jul, 2021) 2104.08919.
  * (11) A. M. Green and B. J. Kavanagh, “Primordial black holes as a dark matter candidate,” feb, 2021.
  * (12) B. Eslam Panah, S. H. Hendi, and Y. C. Ong, “Black hole remnant in massive gravity,” Physics of the Dark Universe 27 (jan, 2020) 1808.07829.
  * (13) P. Chen and R. J. Adler, “Black Hole Remnants and Dark Matter,” Nuclear Physics B - Proceedings Supplements 124 (July, 2003) 103–106. arXiv:gr-qc/0205106.
  * (14) P. Chen, “Inflation Induced Planck-Size Black Hole Remnants As Dark Matter,” 0406514.
  * (15) J. V. Narlikar, K. M. V. Appa Rao, and N. Dadhich, “High energy radiation from white holes,” Nature 251 (Oct., 1974) 590–591. Number: 5476 Publisher: Nature Publishing Group.
  * (16) P. Hájíček and C. Kiefer, “Singularity avoidance by collapsing shells in quantum gravity,” International Journal of Modern Physics D 10 (Dec., 2001) 775–779. Publisher: World Scientific Publishing Co.
  * (17) M. Ambrus and P. Hájíček, “Quantum superposition principle and gravitational collapse: Scattering times for spherical shells,” Physical Review D 72 (Sept., 2005) 064025. Publisher: American Physical Society Author: Hajiçek.
  * (18) J. Olmedo, S. Saini, and P. Singh, “From black holes to white holes: a quantum gravitational, symmetric bounce,” Classical and Quantum Gravity 34 (Oct., 2017) 225011. Publisher: IOP Publishing.
  * (19) M. Christodoulou, C. Rovelli, S. Speziale, and I. Vilensky, “Planck-Star Tunnelling-Time: an Astrophysically Relevant Observable from Background-Free Quantum Gravity,” Physical Review D 94 (Oct., 2016) 084035. arXiv:1605.05268 [astro-ph, physics:gr-qc, physics:hep-th].
  * (20) M. Christodoulou and F. D’Ambrosio, “Characteristic Time Scales for the Geometry Transition of a Black Hole to a White Hole from Spinfoams,” Tech. Rep. arXiv:1801.03027, arXiv, Mar., 2018.  arXiv:1801.03027 [gr-qc] type: article.
  * (21) F. D’Ambrosio, M. Christodoulou, P. Martin-Dussaud, C. Rovelli, and F. Soltani, “End of a black hole’s evaporation,” Physical Review D 103 (may, 2021).
  * (22) F. Soltani, C. Rovelli, and P. Martin-Dussaud, “End of a black hole’s evaporation. II.,” Physical Review D 104 (may, 2021) 106014.
  * (23) C. Rovelli and F. Vidotto, “Small black/white hole stability and dark matter,” Universe 4 (Nov., 2018) 127. arXiv:1805.03872 [gr-qc, physics:physics].
  * (24) T. Banks, A. Dabholkar, M. R. Douglas, and M. O’Loughlin, “Are horned particles the end point of Hawking evaporation?,” Physical Review D 45 (May, 1992) 3607–3616. Publisher: American Physical Society.
  * (25) S. B. Giddings and A. Strominger, “Dynamics of extremal black holes,” Physical Review D 46 (July, 1992) 627–637. Publisher: American Physical Society.
  * (26) T. Banks, M. O’Loughlin, and A. Strominger, “Black hole remnants and the information puzzle,” Physical Review D 47 (May, 1993) 4476–4482. Publisher: American Physical Society.
  * (27) S. B. Giddings, “Constraints on black hole remnants,” Physical Review D 49 (Jan., 1994) 947–957. Publisher: American Physical Society.
  * (28) T. Banks, “Lectures on black holes and information loss,” Nuclear Physics B - Proceedings Supplements 41 (Apr., 1995) 21–65.
  * (29) M. Christodoulou and C. Rovelli, “How big is a black hole?,” Physical Review D 91 (2015) 064046, 1411.2854.
  * (30) T. De Lorenzo and A. Perez, “Improved Black Hole Fireworks: Asymmetric Black-Hole-to-White-Hole Tunneling Scenario,” Physical Review D 93 (June, 2016) 124018. arXiv:1512.04566 [gr-qc, physics:hep-th].
  * (31) C. Rovelli, “Black holes have more states than those giving the Bekenstein-Hawking entropy: a simple argument,” 1710.00218.
  * (32) C. Rovelli, “The subtle unphysical hypothesis of the firewall theorem,” Entropy 21 (2019), no. 9, 1902.03631.
  * (33) J. Maldacena, “Recent progress on the black hole information paradox,” in Strings 2020, Cape Town, South Africa.  2020\.
  * (34) J. Preskill, “Do Black Holes Destroy Information?,” Tech. Rep. arXiv:hep-th/9209058, arXiv, Sept., 1992.  arXiv:hep-th/9209058 type: article.
  * (35) B. J. Carr and S. W. Hawking, “Black Holes in the Early Universe,” Monthly Notices of the Royal Astronomical Society 168 (1974), no. 2, 399–415, arXiv:1011.1669v3.
  * (36) L. Derome, A. Barrau, and G. Boudoul, “An Improved Gamma-Ray Limit on the Density of Primordial Black Holes,”. Conference Name: International Cosmic Ray Conference ADS Bibcode: 2003ICRC….3.1697B.
  * (37) F. Schwabl, Statistical Mechanics.  Springer Science & Business Media, June, 2006.  Google-Books-ID: 7VnKAW284PgC.
  * (38) V. V. Skobelev, “Entropy and energy of a photon gas in an n-dimensional space,” Russian Physics Journal 55 (Feb., 2013) 999–1004.
  * (39) D. Marolf, “The Black Hole information problem: past, present, and future,” Reports on Progress in Physics 80 (2017) 092001, 1703.02143.
  * (40) S. W. Hawking, “Black hole explosions?,” Nature 248 (mar, 1974) 30–31.
  * (41) A. King, “Radiative Processes in Astrophysics,” Physics Bulletin 31 (1980), no. 10, 359–359.
  * (42) B. C. Lacki, “Interpreting the low-frequency radio spectra of starburst galaxies: A pudding of strömgren spheres,” Monthly Notices of the Royal Astronomical Society 431 (jun, 2013) 3003–3024, 1206.7100v3.
  * (43) C. Rovelli, “A generally covariant quantum field theory and a prediction on quantum measurements of geometry,” Nuclear Physics, Section B 405 (sep, 1993) 797–815.
  * (44) L. Smolin, “Finite diffeomorphism-invariant observables in quantum gravity,” Physical Review D 49 (feb, 1994) 4028–4040, 9302011v1.
  * (45) C. Rovelli and L. Smolin, “Discreteness of area and volume in quantum gravity,” Nucl. Phys. B442 (1995) 593–622, 9411005.
  * (46) B. R. Martin and G. Shaw, Nuclear and Particle Physics: An Introduction.  John Wiley & Sons, Apr., 2019.  Google-Books-ID: pWmEDwAAQBAJ.
  * (47) A. Ashtekar and A. Barrau, “Loop quantum cosmology: From pre-inflationary dynamics to observations,” Class. Quant. Grav. 32 (2015), no. 23, 234001. _eprint: 1504.07559.
  * (48) R. Brandenberger and P. Peter, “Bouncing Cosmologies: Progress and Problems,” Found. Phys. 47 (2017), no. 6, 797–850. _eprint: 1603.05834.
  * (49) A. Ijjas and P. J. Steinhardt, “Bouncing Cosmology made simple,” Classical and Quantum Gravity 35 (July, 2018) 135004. arXiv:1803.01961 [astro-ph, physics:gr-qc].
  * (50) A. Barrau, L. Ferdinand, K. Martineau, and C. Renevey, “Closer look at white hole remnants,” Physical Review D 103 (Feb., 2021) 043532.


# Supplementary Materials of ”Momentum-space Spin Anti-vortex and Spin
Transport in Monolayer Pb”

Kaijie Yang  Department of Physics, the Pennsylvania State University,
University Park, PA 16802  Yuanxi Wang  2-Dimensional Crystal Consortium,
Materials Research Institute, The Pennsylvania State University, University
Park, PA 16802  Department of Physics, University of North Texas, Denton, TX
76203  Chao-Xing Liu<EMAIL_ADDRESS>Department of Physics, the Pennsylvania
State University, University Park, PA 16802

##  I Numerical linear response of spin transport with short range disorder

###  I.1 formalism

The formalism for CISP $\chi_{ij}$ and SHE $\sigma^{k}_{ij}$ in orbital basis
are derived from the linear response theory Mahan (2013) and are given in the
main text by

$\chi_{ij}=-\frac{1}{2\pi}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\mathop{\mathrm{Tr}}S_{i}G^{R}({\textbf{{k}}},E_{f})\Gamma_{j}({\textbf{{k}}},E_{f})G^{A}({\textbf{{k}}},E_{f})$
(1)

and

$\sigma^{k}_{ij}=-\frac{1}{2\pi}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\mathop{\mathrm{Tr}}J^{k}_{i}({\textbf{{k}}})G^{R}({\textbf{{k}}},E_{f})\Gamma_{j}({\textbf{{k}}},E_{f})G^{A}({\textbf{{k}}},E_{f})$
(2)

where the Green’s functions are defined as

$G^{R}({\textbf{{k}}},E)=(E-H({\textbf{{k}}})-\Sigma^{R}(E))^{-1}.$ (3)

The integrand

$\chi_{ij}({\textbf{{k}}})=-\frac{1}{2\pi}\mathop{\mathrm{Tr}}S_{i}G^{R}({\textbf{{k}}},E_{f})\Gamma_{j}({\textbf{{k}}},E_{f})G^{A}({\textbf{{k}}},E_{f})$
(4)

is defined as the contribution over the Brillouin zone (BZ) to CISP. The
contribution over BZ to SHE is also defined similarly for the integrand

$\sigma^{k}_{ij}({\textbf{{k}}})=-\frac{1}{2\pi}\mathop{\mathrm{Tr}}J^{k}_{i}({\textbf{{k}}})G^{R}({\textbf{{k}}},E_{f})\Gamma_{j}({\textbf{{k}}},E_{f})G^{A}({\textbf{{k}}},E_{f}).$
(5)

The self energies $\Sigma^{R}({\textbf{{k}}},E)$ in the Green’s functions come
from short range disorder. The disorder has the $\delta$-like potential
$V(r)=\sum_{i}V_{0}\delta({\textbf{{r}}}-{\textbf{{R}}}_{i})$ and the density
$n_{i}$ where $V_{0}$ is the disorder strength and $R_{i}$ is the disorder
position. After disorder average, the self energies in the self-consistent
Born approximation are

$\Sigma^{R}(E)=n_{i}V_{0}^{2}\int\frac{d^{2}{\textbf{{k}}}^{\prime}}{(2\pi)^{2}}G^{R}({\textbf{{k}}}^{\prime},E).$
(6)

The vertex correction
$\Gamma_{j}({\textbf{{k}}},E)=J_{j}({\textbf{{k}}})+\gamma_{j}(E)$ is written
in the ladder approximation as in the main text

$\gamma_{j}(E)=n_{i}V_{0}^{2}\int\frac{d^{2}{\textbf{{l}}}}{(2\pi)^{2}}G^{R}({\textbf{{l}}},E)(J_{j}({\textbf{{l}}})+\gamma_{j}(E))G^{A}({\textbf{{l}}},E).$
(7)

The $\gamma_{j}(E)$ is independent of momenta because the $\delta$-like
disorder potential $V({\bf r})$ becomes a constant
$V({\textbf{{k}}}-{\textbf{{l}}})$ for all transferred momenta after the
transform to the momentum space.

When the disorder is dilute, the formalism transformed to energy eigenspace of
Hamiltonian can be reduced to the Fermi surfaces. The self-consistent Born
approximation for self energies is reduced to first-order Born approximation
and reads

$\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E)=n_{i}V_{0}^{2}\sum_{n^{\prime}}\int\frac{d^{2}{\textbf{{k}}}^{\prime}}{(2\pi)^{2}}(-\pi)\delta(E-\epsilon_{{\textbf{{k}}}^{\prime}n^{\prime}})|\langle{\textbf{{k}}}n|{\textbf{{k}}}^{\prime}n^{\prime}\rangle|^{2}.$
(8)

The $|{\textbf{{k}}}n\rangle$ are the eigen states of Hamiltonian and
$\epsilon_{{\textbf{{k}}}n}$ are the corresponding eigen energies. The
imaginary part of the bare Green’s function in energy eigenspace
$G^{R}_{0}({\textbf{{k}}}n,E)=(E-\epsilon_{{\textbf{{k}}}n}+i0^{+})^{-1}$ is
approximated by $-\pi\delta(E-\epsilon_{{\textbf{{k}}}n})$. The
$\delta$-function shows the contribution to self-energies come from the Fermi
surface. The off-diagonal parts of self energies in eigen-energy space are
dropped. Under the time reversal symmetry $\mathcal{T}$, the diagonal self-
energies become

$\begin{split}\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E)&=n_{i}V_{0}^{2}\sum_{n^{\prime}}\int\frac{d^{2}{\textbf{{k}}}^{\prime}}{(2\pi)^{2}}(-\pi)\delta(E-\epsilon_{{\textbf{{k}}}^{\prime}n^{\prime}})|\langle{\textbf{{k}}}n|\mathcal{T}^{\dagger}\mathcal{T}|{\textbf{{k}}}^{\prime}n^{\prime}\rangle|^{2}.\\\
&=n_{i}V_{0}^{2}\sum_{n^{\prime}}\int\frac{d^{2}{\textbf{{k}}}^{\prime}}{(2\pi)^{2}}(-\pi)\delta(E-\epsilon_{{\textbf{{k}}}^{\prime}n^{\prime}})|\langle-{\textbf{{k}}}^{\prime}n|-{\textbf{{k}}}n^{\prime}\rangle|^{2}\\\
&=\mathop{\mathrm{Im}}\Sigma^{R}(-{\textbf{{k}}}n,E).\end{split}$ (9)

The relaxation time is

$\tau_{{\textbf{{k}}}n}(E)=\frac{-1}{2\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E)}.$
(10)

To evaluate the CISP in Eq. (1), we consider the approximation

$G^{R}({\textbf{{k}}}n,E)G^{A}({\textbf{{k}}}n^{\prime},E)\approx\delta_{nn^{\prime}}\delta(E-\epsilon_{{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E)},$
(11)

in which the off-diagonal part $n\neq n^{\prime}$ is of the order
$\mathop{\mathrm{Im}}\Sigma^{R}/(\epsilon_{{\textbf{{k}}}n}-\epsilon_{{\textbf{{k}}}n^{\prime}})$,
as compared to diagonal part $n=n^{\prime}$, and thus is neglected. Under this
approximation, the response coefficient $\chi_{ij}$ of the CISP in energy
eigenspace reads

$\chi_{ij}=-\frac{1}{2\pi}\sum_{n}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\langle{\textbf{{k}}}n|S_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n\rangle\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E_{f})}.$
(12)

Below we will also consider the integrand

$\chi_{ij}({\textbf{{k}}})=-\frac{1}{2\pi}\langle{\textbf{{k}}}n|S_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n\rangle\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E_{f})},$
(13)

which reveals the momentum-resolved contribution around the Fermi surfaces of
the band n. Under the similar approximation, the vertex correction reads

$\gamma_{j}({\textbf{{k}}}n,E)=n_{i}V_{0}^{2}\sum_{m}\int\frac{d^{2}{\textbf{{l}}}}{(2\pi)^{2}}(J_{j}({\textbf{{l}}}m)+\gamma_{j}({\textbf{{l}}}m,E))\lvert\langle{\textbf{{k}}}n|{\textbf{{l}}}m\rangle\rvert^{2}\delta(E-\epsilon_{{\textbf{{l}}}m})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{l}}}m,E)}$
(14)

where
$\gamma_{j}({\textbf{{k}}}n,E)=\langle{\textbf{{k}}}n|\gamma_{j}(E)|{\textbf{{k}}}n\rangle$
and
$J_{j}({\textbf{{k}}}n)=\langle{\textbf{{k}}}n|J_{j}({\textbf{{k}}})|{\textbf{{k}}}n\rangle$.
After obtaining $\gamma_{j}({\textbf{{k}}}n,E)$, the off-diagonal part of the
vertex correction is

$\langle{\textbf{{k}}}n|\gamma_{j}(E)|{\textbf{{k}}}n^{\prime}\rangle=n_{i}V_{0}^{2}\sum_{m}\int\frac{d^{2}{\textbf{{l}}}}{(2\pi)^{2}}(J_{j}({\textbf{{l}}}m)+\gamma_{j}({\textbf{{l}}}m,E))\langle{\textbf{{k}}}n|{\textbf{{l}}}m\rangle\langle
lm|{\textbf{{k}}}n^{\prime}\rangle\delta(E-\epsilon_{{\textbf{{l}}}m})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{l}}}m,E)}.$
(15)

The vertex correction can also be calculated for the spin operator in CISP as

$\langle{\textbf{{k}}}n|\gamma^{S}_{i}(E)|{\textbf{{k}}}n\rangle=n_{i}V_{0}^{2}\sum_{m}\int\frac{d^{2}{\textbf{{l}}}}{(2\pi)^{2}}(\langle{\textbf{{l}}}m|S_{i}|{\textbf{{l}}}m\rangle+\langle{\textbf{{l}}}m|\gamma^{S}_{i}(E)|{\textbf{{l}}}m\rangle)\lvert\langle{\textbf{{k}}}n|{\textbf{{l}}}m\rangle\rvert^{2}\delta(E-\epsilon_{{\textbf{{l}}}m})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{l}}}m,E)}.$
(16)

The spin operator with the vertex correction is

$\langle{\textbf{{k}}}n|\tilde{S}_{i}(E)|{\textbf{{k}}}n\rangle=\langle{\textbf{{k}}}n|S_{i}+\gamma^{S}_{i}(E)|{\textbf{{k}}}n\rangle.$
(17)

For SHE, Eq. (2) is divided into diagonal parts and off diagonal parts in
energy eigenspace as

$\sigma^{k}_{ij}=\sigma^{D}+\sigma^{OD}.$ (18)

The diagonal (intraband) contribution is

$\sigma^{D}=-\frac{1}{2\pi}\sum_{n}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\langle{\textbf{{k}}}n|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}|{\textbf{{k}}}n\rangle\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E_{f})}.$
(19)

Under $\mathcal{T}$, the diagonal part becomes

$\begin{split}\sigma^{D}&=-\frac{1}{2\pi}\sum_{n}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\langle-{\textbf{{k}}}n|(\mathcal{T}J^{k}_{i}\mathcal{T}^{\dagger})|-{\textbf{{k}}}n\rangle^{*}\langle-{\textbf{{k}}}n|(\mathcal{T}\Gamma_{j}\mathcal{T}^{\dagger})|-{\textbf{{k}}}n\rangle^{*}\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E_{f})}\\\
&=-\frac{1}{2\pi}\sum_{n}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\langle{\textbf{{k}}}n|(J^{k}_{i})|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|(-\Gamma_{j})|{\textbf{{k}}}n\rangle\delta(E_{f}-\epsilon_{-{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}(-{\textbf{{k}}}n,E_{f})}\\\
&=-\sigma^{D}.\end{split}$ (20)

So the intraband contribution to SHE is zero. The off diagonal part of the
real SHE conductivity reads

$\sigma^{k}_{ij}=\sigma^{OD}=-\frac{1}{2\pi}\sum_{n\neq
n^{\prime}}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}(-)\mathop{\mathrm{Im}}(G^{R}({\textbf{{k}}}n,E_{f})G^{A}({\textbf{{k}}}n^{\prime},E_{f}))\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle)$
(21)

where the term
$\mathop{\mathrm{Re}}(G^{R}({\textbf{{k}}}n,E)G^{A}({\textbf{{k}}}n^{\prime},E))$
is neglected because it is of the order
$\mathop{\mathrm{Im}}(\Sigma^{R})^{2}/(\epsilon_{{\textbf{{k}}}n}-\epsilon_{{\textbf{{k}}}n^{\prime}})^{2}$
compared to the
$\mathop{\mathrm{Im}}(G^{R}({\textbf{{k}}}n,E)G^{A}({\textbf{{k}}}n^{\prime},E))$.
Adopting the approximation

$\mathop{\mathrm{Im}}(G^{R}({\textbf{{k}}}n,E)G^{A}({\textbf{{k}}}n^{\prime},E))\approx(-\pi)(\delta(E-\epsilon_{{\textbf{{k}}}n})+\delta(E-\epsilon_{{\textbf{{k}}}n^{\prime}}))\frac{1}{\epsilon_{{\textbf{{k}}}n}-\epsilon_{{\textbf{{k}}}n^{\prime}}},$
(22)

we find the SHE from the Fermi surface contribution is

$\sigma^{k}_{ij}=-\frac{1}{2\pi}\sum_{n\neq
n^{\prime}}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}2\pi\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{1}{\epsilon_{{\textbf{{k}}}n}-\epsilon_{{\textbf{{k}}}n^{\prime}}}\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle).$
(23)

We also define the momentum-resolved contribution to the SHE as

$\sigma^{k}_{ij}({\textbf{{k}}}n)=-\frac{1}{2\pi}\sum_{n^{\prime}}2\pi\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{1}{\epsilon_{{\textbf{{k}}}n}-\epsilon_{{\textbf{{k}}}n^{\prime}}}\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle)$
(24)

at the Fermi surface of the band n.

In our numerical simulations, we directly evaluate Eq. (1) and Eq. (2) over
the first Brillouin zone. The numerical results can be well understood from
the simplified expressions of Eq. (12) and Eq. (23), which only depends on the
contribution from the Fermi energy.

###  I.2 Constructing the Tight Binding Model from DFT

Density functional theory (DFT) calculations were based on the relaxed
structure of monolayer Pb on SiC, where Pb coverage is 100%, i.e. one Pb per
surface Si atom. We focus on this surface phase here since it is the
thermodynamic ground state in a wide range of Pb chemical potentials compared
with competing surface phases where Pb coverage is partial. Within candidate
structures satisfying 100% coverage, the registry where each Pb is directly on
top of a surface Si site is preferred against the two other high-symmetry
registries (Pb projecting onto the topmost C site or a “hollow” site, as
defined in Ref. [Briggs  _et al._ , 2020]). The relaxed Pb structure features
one out of every three Pb atom buckling out-of-plane. $C_{3v}$ symmetry is
preserved even when this buckling is considered, since our supercell was
chosen to be minimal and includes just three Pb atoms.

For DFT-level structural relaxations, the Perdew-Burke-Ernzerhof (PBE)
parametrization of the exchange-correlation functional at the generalized
gradient approximation level was used Perdew  _et al._ (1996, 1997). We
employed projector augmented wave pseudopotentials Blöchl (1994); Kresse and
Joubert (1999) as implemented in the Vienna Ab-initio Simulation Package
(VASP) Kresse and Furthmüller (1996). All relaxations were performed with the
semi-empirical DFT-D3 van der Waals correction scheme Grimme  _et al._ (2010)
and with a force threshold of 0.01 eV/Å. Brillouin zone sampling densities
were set to be equivalent to $12\times 12\times 1$ for a minimal SiC surface
unit cell with the in-plane lattice constant $a=3.19$Å. The SiC slab thickness
was 6 SiC units, with hydrogen capping included for the other surface
(opposite to the Pb-covered surface).

To construct the tight-binding model used in the main text, we adopted a
structural approximation – a relaxed minimal SiC surface unit cell with one Pb
atom – to perform preparatory self-consistency DFT calculations for Wannier
interpolation. These calculations were performed using the Quantum Espresso
package Giannozzi  _et al._ (2009) with projector augmented wave
pseudopotentials and with k-point grids of $12\times 12\times 1$, all other
convergence parameters being the same as those used in the VASP calculations.
Maximally localized Wannier functions and tight-binding parameters Marzari and
Vanderbilt (1997) were constructed using the Wannier90 code Mostofi  _et al._
(2008). Band disentanglement followed the standard procedure outlined in Ref.
[Souza  _et al._ , 2001]: a disentanglement window was chosen to cover the
entire energy where bands with all four types of orbital characters (see main
text) were present; a frozen window (where Bloch state are forced to be
preserved) were chosen to cover a smaller energy range where bands with
exclusively all four types of orbital characters are present.

###  I.3 Numerical Realization

This section discusses the numerical methods based on Ref. [Wei  _et al._ ,
2021] to evaluate the formula in Sec. I.1.

For the units used in numerical calculation, the length unit is
$k_{0}^{-1}=1\text{\AA}$ and the energy unit is $\epsilon_{0}=1eV$. The units
for other physical quantities are determined from these two units accordingly.
For example, the unit for CISP is $\chi_{0}=k_{0}e\hbar/\epsilon_{0}$. The
unit for $n_{i}V_{0}^{2}$ in Eq. (6) and Eq. (7) is
$k_{0}^{2}\epsilon_{0}^{2}$. With this unit, the $n_{i}V_{0}^{2}$ is taken to
be $0.05$ in the main text.

The first BZ integration in Eq. (1), Eq. (2), Eq. (6) and Eq. (7) is replaced
by summation over a mesh. The first BZ is a hexagon with three $K$ points
located at

$K=(4\pi/3,0)\quad K_{1}=(4\pi/3\times\cos 2\pi/3,4\pi/3\times\sin
2\pi/3)\quad K_{3}=(4\pi/3\times\cos 4\pi/3,4\pi/3\times\sin 4\pi/3)$ (25)

The mesh points are taken to be $k_{1}K+k_{2}K_{1}$, $k_{1}K_{1}+k_{2}K_{2}$
and $k_{1}K_{2}+k_{2}K$, respecting $C_{3}$ symmetry to ensure the current
conserving (See Sec. I.4). The $k_{1}$ and $k_{2}$ are two real numbers in the
range $[0,1]$ with the mesh step $\Delta k=0.005\times 3/4\pi$. The
corresponding energy change of adjacent mesh points in spectrum
$\Delta\epsilon\sim 0.005$ is roughly an order smaller than the self-energies
$\sim n_{i}V_{0}^{2}$. Thus, the BZ integration in Eq. (1), Eq. (2), Eq. (6)
and Eq. (7) could be replaced by

$\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}f({\textbf{{k}}})\rightarrow\sum_{{\textbf{{k}}}_{i}}\frac{(\Delta
k)^{2}\sin 2\pi/3}{(2\pi)^{2}}f({\textbf{{k}}}_{i})$ (26)

where $sin(2\pi/3)$ comes from the Jacobin of a parallelogram mesh with the
angle $2\pi/3$. In total, there are $n_{k}\sim 10^{6}$ points for the mesh.

The self-consistent self-energies and Green’s functions in Eq. (3) and Eq. (6)
are calculated in an iterative manner. The initial values for self-energies is
taken to be $-in_{i}V_{0}^{2}/2$, same as self-energies for 2DEG. After 12
iteration, the difference in self-energies is $\sim 10^{-5}\ll\Delta\epsilon$.
Thus, the self-energies are considered as converged to the self-consistent
values.

The vertex corrections in Eq. (7) are solved in the following way. We move
$\gamma(E)$ out of the matrix multiplication $G^{R}\gamma G^{A}$ by rewriting
Eq. (7) in the form of the tensor product as

$\gamma_{j,mn}(E)=\sum_{m^{\prime},n^{\prime}}n_{i}V_{0}^{2}\int\frac{d^{2}{\textbf{{l}}}}{(2\pi)^{2}}(G^{R}({\textbf{{l}}},E)\otimes
G^{A}({\textbf{{l}}},E)^{T})_{mn,m^{\prime}n^{\prime}}(J_{j,m^{\prime}n^{\prime}}({\textbf{{l}}})+\gamma_{j,m^{\prime}n^{\prime}}(E))$
(27)

where $n,m$ are the basis vector indices, $(G^{R}\otimes G^{A\
T})_{mn,m^{\prime}n^{\prime}}=G^{R}_{m,m^{\prime}}G^{A}_{n^{\prime},n}$ and
$J_{j,mn}({\textbf{{l}}})$ is treated as a vector with the $(m-1)\times
n_{b}+n$th element as $\langle{\textbf{{l}}}m|J_{j}|{\textbf{{l}}}n\rangle$.
Here $n_{b}$ is the number of basis. With matrix inversion, we find that the
vertex correction is solved as

$\gamma_{j}(E)=(1-n_{i}V_{0}^{2}\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}(G^{R}({\textbf{{k}}},E)\otimes
G^{A}({\textbf{{k}}},E)^{T}))^{-1}n_{i}V_{0}^{2}\int\frac{d^{2}{\textbf{{l}}}}{(2\pi)^{2}}(G^{R}({\textbf{{l}}},E)\otimes
G^{A}({\textbf{{l}}},E)^{T})J_{j}({\textbf{{l}}}).$ (28)

The dimension of the matrices are of order $\sim O(n_{k}n_{b}^{2})$. The basis
for calculation are the Bloch states transformed from wannier orbitals used in
DFT calculation. In these basis, $\gamma_{j}(E)$ is independent of momenta for
short range disorder and only has $n_{b}^{2}$ elements compared to
$n_{k}n_{b}^{2}$ elements in energy eigenbasis.

Then we are going to show how to deal with the the integration with the delta
function $\delta(E-\epsilon_{{\textbf{{k}}}n})$ in Eq. (12), Eq. (23), Eq. (8)
and Eq. (14). The delta function $\delta(E-\epsilon_{{\textbf{{k}}}n})$ is
first integrated out by radial integration as

$\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\delta(E-\epsilon_{{\textbf{{k}}}n})f({\textbf{{k}}})\rightarrow\int\frac{k(\theta)d\theta}{(2\pi)^{2}}\frac{1}{|\partial\epsilon_{{\textbf{{k}}}n}/\partial
k|_{k=k(\theta)}}f(k(\theta),\theta).$ (29)

The angular integration is then replaced by a discrete sum as

$\int\frac{d^{2}{\textbf{{k}}}}{(2\pi)^{2}}\delta(E-\epsilon_{{\textbf{{k}}}n})f({\textbf{{k}}})\rightarrow\sum_{\alpha,\theta_{i}}\frac{k_{n}^{\alpha}(\theta_{i})\Delta\theta}{(2\pi)^{2}}\frac{1}{|\partial\epsilon_{{\textbf{{k}}}n}/\partial
k|_{k=k_{n}^{\alpha}(\theta_{i})}}f(k_{n}^{\alpha}(\theta_{i}),\theta_{i}).$
(30)

The $k_{i}^{\alpha}(\theta)$ is the interpolated wave vector at the Fermi
surface, which is characterized by the parameter $\theta$ (the corresponding
angle of the momentum) for the $\alpha$ branch of the band $n$ Fermi surface.

Rewrite self-energies in Eq. (8) and vertex correction in Eq. (14) with the
interpolated Fermi surfaces as

$\mathop{\mathrm{Im}}\Sigma^{R}(k_{n}^{\alpha}(\theta_{i}),E)=n_{i}V_{0}^{2}\sum_{n^{\prime},\alpha^{\prime},\theta_{i}^{\prime}}\frac{k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})\Delta\theta}{(2\pi)^{2}}\frac{-\pi}{|\partial\epsilon_{{\textbf{{k}}}n}/\partial
k|_{k=k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})}}\lvert\langle
k_{n}^{\alpha}(\theta_{i})|k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})\rangle\rvert^{2}.$
(31)
$\gamma_{j}(k_{n}^{\alpha}(\theta_{i}),E)=n_{i}V_{0}^{2}\sum_{n^{\prime},\alpha^{\prime},\theta_{i}^{\prime}}\frac{k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})\Delta\theta}{(2\pi)^{2}}\frac{-\pi}{|\partial\epsilon_{{\textbf{{k}}}n}/\partial
k|_{k=k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})}}\lvert\langle
k_{n}^{\alpha}(\theta_{i})|k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})\rangle\rvert^{2}\frac{1}{\mathop{\mathrm{Im}}\Sigma^{R}(k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime}),E)}(J_{j}(k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime}))+\gamma_{j}(k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime}),E)).$
(32)

###  I.4 Current conservation

In this section, we focus on solving the vertex correction in Eq. (7) and Eq.
(14). We are going to show that the matrix to be inversed for vertex
correction has zero determinant and its relation to current conservation.

The zero determinant of the matrix in Eq. (28) to be inversed is illustrated
with Eq. (31) and Eq. (32) of Fermi surface contributions. By denoting

$\begin{split}M_{{\textbf{{k}}}^{\prime}n^{\prime}}&=M_{k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})}=\frac{k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})\Delta\theta}{(2\pi)^{2}}\frac{-\pi}{|\partial\epsilon_{{\textbf{{k}}}n}/\partial
k|_{k=k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})}}\\\
V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}&=V(k_{n}^{\alpha}(\theta_{i}),k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime}))=n_{i}V_{0}^{2}\lvert\langle
k_{n}^{\alpha}(\theta_{i})|k_{n^{\prime}}^{\alpha^{\prime}}(\theta_{i}^{\prime})\rangle\rvert^{2},\end{split}$
(33)

we have the self energies and the vertex correction as

$\displaystyle\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E)=\sum_{{\textbf{{k}}}^{\prime}n^{\prime}}V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}^{\prime}n^{\prime}}$
(34)
$\displaystyle\gamma_{j}({\textbf{{k}}}n,E)=\sum_{{\textbf{{k}}}^{\prime}n^{\prime}}V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}^{\prime}n^{\prime}}\frac{1}{\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}^{\prime}n^{\prime},E)}(J_{j}({\textbf{{k}}}^{\prime}n^{\prime})+\gamma_{j}({\textbf{{k}}}^{\prime}n^{\prime},E)).$
(35)

The matrix $V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}$ is
symmetric. The matrix to be inversed for the vertex correction is

$I_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}=\delta_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}-V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}^{\prime}n^{\prime}}\frac{1}{\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}^{\prime}n^{\prime},E)}.$
(36)

For the singular value decomposition (SVD) of the matrix $I$, $I$ has a left
eigen vector $e_{L0}=M_{{\textbf{{k}}}n}$ and a right eigen vector
$e_{R0}=\mathop{\mathrm{Im}}\Sigma({\textbf{{k}}}^{\prime}n^{\prime},E)$ as
shown by

$\displaystyle\sum_{{\textbf{{k}}}n}M_{{\textbf{{k}}}n}I_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}=M_{{\textbf{{k}}}^{\prime}n^{\prime}}-(\sum_{{\textbf{{k}}}n}M_{kn}V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}})M_{{\textbf{{k}}}^{\prime}n^{\prime}}\frac{1}{\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}^{\prime}n^{\prime},E)}=0$
(37)
$\displaystyle\sum_{{\textbf{{k}}}^{\prime}n^{\prime}}I_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}\mathop{\mathrm{Im}}\Sigma({\textbf{{k}}}^{\prime}n^{\prime},E)=\mathop{\mathrm{Im}}\Sigma({\textbf{{k}}}n,E)-\sum_{{\textbf{{k}}}^{\prime}n^{\prime}}V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}^{\prime}n^{\prime}}=0$
(38)

with Eq. (34). Both eigenvectors have the eigenvalue 0. As a result, the
matrix $I_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}$ has zero
determinant.

The existence of a solution to the vertex correction in Eq. (35) implies the
current conservation physically. To see that, we re-write the vertex
correction in the form of the SVD decomposition of the matrix $I$ as

$\begin{pmatrix}e_{R0}&e_{R1}&\cdots&e_{Rn}\end{pmatrix}\begin{pmatrix}0&&&\\\
&\lambda_{1}&&\\\ &&\ddots&\\\
&&&\lambda_{n}\end{pmatrix}\begin{pmatrix}e_{L0}^{\dagger}\\\
e_{L1}^{\dagger}\\\ \vdots\\\
e_{Ln}^{\dagger}\end{pmatrix}\gamma_{j}=\tilde{J_{j}}$ (39)

where $e_{Li}\ e_{Ri}\ \lambda_{i}$ are other left eigenvectors, right
eigenvectors, eigenvales of $I$, statisfying
$e_{Li}^{\dagger}e_{Rj}=\delta_{ij}$. The elements of $\gamma_{j}$ is
$\gamma_{j}({\textbf{{k}}}n)$ in Eq. (14) and the elements of $\tilde{J}_{j}$
is
$\sum_{{\textbf{{k}}}^{\prime}n^{\prime}}V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}^{\prime}n^{\prime}}J_{j}({\textbf{{k}}}^{\prime}n^{\prime})/\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}^{\prime}n^{\prime},E)$.
The condition of the existence of $\gamma_{j}$ is

$e_{L0}^{\dagger}\tilde{J_{j}}=\sum_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}n}V_{{\textbf{{k}}}n,{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}^{\prime}n^{\prime}}J_{j}({\textbf{{k}}}^{\prime}n^{\prime})/\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}^{\prime}n^{\prime},E)=\sum_{{\textbf{{k}}}^{\prime}n^{\prime}}M_{{\textbf{{k}}}^{\prime}n^{\prime}}J_{j}({\textbf{{k}}}^{\prime}n^{\prime})=0.$
(40)

Restoring to the integral form, the condition becomes

$\int\frac{\mathrm{d}^{2}{\textbf{{k}}}}{(2\pi)^{2}}\delta(E-\epsilon_{{\textbf{{k}}}n})J_{j}({\textbf{{k}}}n)=0.$
(41)

Physically, this means the total current of the electrons in equilibrium is
zero. Since the charge $\rho$ is conserved, the current conservation
$\partial_{t}\rho+\nabla\cdot{\textbf{{J}}}=0$ is statisfied.

To solve the vertex correction, the inverse in Eq. (35) shall be replaced by
Moore-Penrose pseudo-inverse Barata and Hussein (2012), given by

$\begin{pmatrix}e_{L0}^{\dagger}\gamma_{j}\\\ e_{L1}^{\dagger}\gamma_{j}\\\
\vdots\\\ e_{Ln}^{\dagger}\gamma_{j}\end{pmatrix}=\begin{pmatrix}0&&&\\\
&1/\lambda_{1}&&\\\ &&\ddots&\\\
&&&1/\lambda_{n}\end{pmatrix}\begin{pmatrix}e_{L0}^{\dagger}\tilde{J_{j}}\\\
e_{L1}^{\dagger}\tilde{J_{j}}\\\ \vdots\\\
e_{Ln}^{\dagger}\tilde{J_{j}}\end{pmatrix}$ (42)

in the SVD decomposition. Another way to solve for the vertex correction with
this non-iterative method is using the finite frequency response
$\chi_{ij}(\omega)$ and then decreasing $\omega$ for a convergent result in
the $\omega\rightarrow 0$ limit.

##  II Numerical Results of Spin Transport

###  II.1 Rashba 2DEG

Figure 1:  (a) The spectrum of Rashba 2DEG. (b) The Fermi surface for
$E_{f}=1$ and the spin texture on the Fermi surface. (c) The spectral weight
shown over the Brillouin zone. (d)(e) The contribution to CISP $\chi_{yx}$ on
the Fermi surface of Eq. (13) without and with vertex correction. (e)(f) The
contribution to CISP $\chi_{yx}$ over the Brillouin zone of Eq. (4) without
and with vertex correction. (h) the relaxation time over the Fermi surfaces
$\tau_{{\textbf{{k}}}n}+20$. (i)(j) The contribution to SHE $\sigma^{z}_{yx}$
on the Fermi surface of Eq. (24) without and with vertex correction. (k)(l)
The contribution to SHE $\sigma^{z}_{yx}$ over the Brillouin zone of Eq. (5)
without and with vertex correction. (m) The y component of the spin directions
on the Fermi surface. The solid (dashed) lines are the spin vertex without
(with) the vertex correction. The spin components are shown in the polar angle
$\theta$ of the Brillouin zone.

To validate our numerical methods, we will first evaluate the spin transport
of 2-dimensional electron gas (2DEG) with Rashba spin-orbit coupling, which
can be directly compared with the analytical solutions in Ref. [Inoue  _et
al._ , 2003].

The Hamiltonian of the Rashba 2DEG is

$H(k)=k^{2}/2+\alpha_{R}(k_{y}\sigma_{x}-k_{x}\sigma_{y}).$ (43)

Fig. 1(a) shows the spectrum of two spin splitting bands, with the parameters
$\alpha_{R}=0.5$ and $E_{f}=1$. Fig. 1(b) shows the Fermi surfaces and the
spin texture on the Fermi surfaces. The Fermi surfaces are two concentric
circles. The spins on the Fermi surfaces of band 1 and 2 are opposite. Fig.
1(c) shows The spectral weight
$-2\mathop{\mathrm{Im}}TrG^{R}({\textbf{{k}}},E_{f})$ with
$n_{i}V_{0}^{2}=0.05$. The two bands are clearly resolved in the dilute
disorder limit.

We directly compare our numerical results for CISP with theoretical values
obtained from Ref. [Inoue  _et al._ , 2003] for a typical set of parameters
listed above. The CISP from the full numerical calculations of the formula Eq.
(1) is $-0.79878$ without vertex correction and $-1.5865$ with vertex
correction. The simplified formula Eq. (12) only with the Fermi surface
contribution gives $-0.79572$ without vertex corrections and $-1.5914$ with
vertex correction. Thus, the two results from different expressions give
nearly identical the CISP response coefficients that also agree with the
theoretical values $-\alpha_{R}/4\pi n_{i}V_{0}^{2}=-0.79578$ without and
$-\alpha_{R}/2\pi n_{i}V_{0}^{2}=-1.5916$ with vertex corrections from Ref.
[Inoue  _et al._ , 2003] within $1\%$ . The error can be reduced by taking
denser meshes.

Next we provide an intuitive physical picture for the general understanding of
the origin of the CISP $\chi_{yx}$. Under the electrical field along the $+x$
direction, the electrons are in the non-equilibrium states and the
corresponding Fermi sphere shifts to the $-x$ direction. For the Fermi
surfaces of the band 1 in Fig. 1 (b) , less electrons occupy the states with
$k_{x}>0$ and spin directions $s_{y}>0$ while more electrons occupy the states
with $k_{x}<0$ and spin directions $s_{y}<0$. As a result, a charge current in
$+x$ direction gives rise to a negative contribution to the CISP $\chi_{yx}$
($s_{y}<0$) for the Fermi surface of the band 1, as shown by the blue color in
Fig. 1(d)(e)(f)(g) . Thus, the electrons on this Fermi surface contributes
$\propto(k_{f}+2\alpha_{R})\tau$. $\tau$ is the relaxation time and constant
over the whole Fermi surfaces for both bands as shown in Fig. 1(h) with values
close to the theoretical one $-n_{i}V_{0}^{2}/2$. The electrons on the Fermi
surface of the band 2 gives a positive contribution (red color) to CISP due to
the opposite spin texture as shown in Fig. 1(d)(e)(f)(g), which is
proportional to $k_{f}\tau$. Adding the contributions from two Fermi surfaces,
we find $\chi_{yx}\propto\alpha_{R}\tau$.

For the vertex correction, it is more convenient to consider the correction to
spin operator instead of current operator. Theoretically, the spin operator
with the vertex correction
$\langle{\textbf{{k}}}n|\tilde{S}_{y}(E)|{\textbf{{k}}}n\rangle$ in Eq. (17)
doubles the bare spin operator matrix elements
$\langle{\textbf{{k}}}n|S_{y}|{\textbf{{k}}}n\rangle$. Thus, CISP with the
vertex correction is $\tilde{\chi}_{yx}=2\chi_{yx}$. Fig. 1(m) shows the
numerical spin operator matrix element from Eq. (17) . The spin operator with
the vertex correction doubles the bare spin operator over the whole Fermi
surfaces, consistent with the theoretical results.

Finally, we show that the numerical results for SHE also agree with
theoretical results. The calculated SHE $\sigma^{z}_{yx}$ from the formula Eq.
(2) are $0.0397802$ without vertex correction and $1.26112\times 10^{-4}$ with
vertex correction. The simplified formula Eq. (23) gives $0.039789$ without
vertex correction and $-2.75749\times 10^{-4}$ with vertex corrections. Both
agree with each other and are consistent with the theoretical values $1/8\pi$
without vertex correction and $0$ with vertex correction from Ref. [Inoue  _et
al._ , 2004] within the error $1\%$.

Fig. 1(i)-(l) shows the contribution to the response coefficient
$\sigma^{z}_{yx}$ of SHE defined by Eq. (24) and Eq. (5). Without vertex
correction, the off-diagonal terms contribute a positive value to SHE (red
color), which is maximal on the y axis and nearly zero on the x-axis. The
diagonal contribution is almost zero as
$\langle{\textbf{{k}}}n|J_{y}^{z}|{\textbf{{k}}}n\rangle$ is zero. When
including the vertex correction, we find the SHE coefficient $\sigma^{z}_{yx}$
is almost zero, as predicted by the analytical theory This reveals the
cancellation of the intrinsic and intrinsic-side-jump contributions to SHE
Sinova  _et al._ (2015).

###  II.2 Spin texture of monolayer Pb

Figure 2:  (a) Spin directions along the $\Gamma-K$ line in 3D. (b) The
projection of states along $\Gamma-K$ to different orbitals. The total is the
sum of all orbitals. The grey region is the anti crossing region. (c) The
contribution to y components of spin directions from different orbitals. (d)
The energy spectrum over 2D Brillouin zone around $K$ in 3D. (e) The y
components of spin directions from all orbitals for different $V_{z}$ from
$-0.6eV$ to $2.0eV$ with an interval $0.2eV$ indicated by the colorful arrow.
(f) The energy spectrum with $V_{z}=1eV$.

In this section, we will discuss the physical properties and origins of the
momentum-space spin anti-vortex, as well as the relation to the Lifshitz
transition of Fermi surfaces.

First, we show the matrix element
$\langle{\textbf{{k}}}n|S_{i}|{\textbf{{k}}}n\rangle$ of spin operator at the
location of spin anti-vortex approaches zero. Fig. 2(a) shows the spin
direction along the $\Gamma-K$ line. At the high symmetry points $\Gamma$ and
$K$, the states are degenerate and the spin directions are not well defined.
We find the spin directions along this line mainly lies in the plane and point
to the $y$ axis as a consequence of the $x$-directional mirror and time
reversal symmetries. As the spin directions near $\Gamma$ and $K$ are
opposite, it has to approach zero at certain momentum in the $\Gamma$-$K$
line, which is labelled as $k_{S}$ in Fig. 2(a).

The spin anti-vortex originates from the anti-crossing between the bands with
opposite spin directions. Fig. 2(b) shows the projection of eigen-states to
different orbitals, defined as
$\langle{\textbf{{k}}}n|P_{i}|{\textbf{{k}}}n\rangle$ with the projection
operator $P_{i}=|p_{i}\rangle\langle p_{i}|$ to the orbital $p_{i}$, along
$\Gamma-K$. Around $k_{S}$, the dominant orbital nature of the bands changes
from the $p_{z}$ orbital of SiC (red line) to the $p_{x}$ and $p_{y}$ orbitals
of Pb (blue lines). As the orbital nature varies rapidly across the anti-
crossing regime, we choose the following condition (the gradient of the
probability difference between different orbitals)

$\lvert\nabla_{\textbf{{k}}}\langle{\textbf{{k}}}n|P_{x}+P_{y}-P_{zSiC}|{\textbf{{k}}}n\rangle\rvert>1$
(44)

to determine the anti-crossing regime, as shown by the grey region in Fig.
2(b) and (c). Fig. 2(c) shows the contribution to $s_{y}$ components from
different orbitals
$\langle{\textbf{{k}}}n|P_{i}s_{y}P_{i}|{\textbf{{k}}}n\rangle$. The
$p_{zSiC}$ orbital carries negative $s_{y}$ while the $p_{x}$ and $p_{y}$
orbitals carry positive $s_{y}$. As the orbital nature of the bands varies
from the $p_{zSiC}$ orbital to the $p_{x}$ and $p_{y}$ orbitals, zero $s_{y}$
occurs in the anti-crossing regime and determines the location of the spin
anti-vortex $k_{S}$.

Besides spin anti-vortex, the band anti-crossing also leads to the Lifshitz
transition of Fermi surfaces. The Lifshitz transition happens at the saddle
points of the spectrum. Fig. 2(d) shows the spectrum around $K$ in the 2D BZ
and the saddle point is labeled by $k_{L}$. In the anti-crossing region
between two black lines, the spectrum shows a loop of local maxima in the
radial direction. Along the tangential direction of this loop, the red points
are the local minima and these minima define the saddle points in the energy
spectrum. The Lifshitz transition happens at these saddle points. In summary,
both spin anti-vortex and Lifshitz transition happen in the anti-crossing
region as shown in Fig. 2(d).

In addition, the location of the spin anti-vortex can controlled by an
external gate that tunes the difference in potentials between the Pb layer and
SiC substrate. To model the effect, we include a constant potential $V_{z}$ to
the $p_{z}$ orbital of SiC and the energy spectrum for $V_{z}=1.0$eV is shown
in Fig. 2(f). Compared to the Fig. 2(d) in the main text for $V_{z}=0$eV, the
$j=1/2$ bands from the SiC for $V_{z}=1$eV have higher energies and the band
anti-crossing region moves away from $K$ point, which leads to the change of
the location of both the Lifshitz transition point $k_{L}$ and the spin anti-
vortex point $k_{S}$. To show the movable range of the spin anti-vortex, we
plot the y component of spin directions on the $\Gamma-K$ line from all
orbitals for different $V_{z}$ in 2(e). The zeros of these lines around
$k_{S}(V_{z}=0)$ indicate the locations of spin anti-vortex points for
different $V_{z}$. The spin anti-vortex could be moved towards $\Gamma$ by
positive $V_{z}$ (blue lines) or $K$ by negative $V_{z}$ (red lines).

###  II.3 Numerical CISP for monolayer Pb

Figure 3:  (a) Fermi surfaces and spin texture for $E_{F}=1eV$. (b) Spectral
weight for $E_{F}=1eV$ over the Brillouin zone. (c) Relaxation time over the
Fermi surfaces. (d)(e) The contribution to CISP on Fermi surfaces without and
with vertex correction. (f) The $s_{y}$ over Fermi surfaces for band 1 and 2.
The red line (blue) is band 1 (2). The solid line is the bare spin vertex. The
dashed line is the corrected spin vertex.  Figure 4:  (a) CISP. The solid
(dashed) line is calculated from the Brillouin zone formula Eq. (1) with
(without) vertex correction. The red circles (blue triangles) are calculated
from the Fermi surface formula Eq. (12). (b)(c) The contribution to CISP and
relaxation time at the intersection points $k_{1}k_{2}k_{3}k_{4}$ between the
Fermi surfaces of band 1 and 2 and the high symmetry lines $\Gamma-K$ and
$\Gamma-M$ lines. (d)(e)(f)(g) the contribution to CISP over the Fermi
surfaces Eq. (13) for $E_{F}=1eV,0eV,-0.1eV,-0.2eV$.

In this section, we will focus on understanding the numerical results of CISP
for monolayer Pb by comparing with Rashba 2DEG at high $E_{f}$ and explaining
the behaviors of CISP when varying $E_{f}$. The numerical calculation by Eq.
(12) at $E_{F}=1eV$ for the band 1 and 2 in Fig. 3(a) is consistent with the
physical picture of Fermi surfaces shift in Rashba 2DEG. To see that, Fig.
3(a) shows the spin texture on the Fermi surfaces of bands 1 and 2, similar to
Rashba 2DEG and opposite to each other. All the calculations of CISP are done
within dilute disorder limits, for which the band broadening induced by
disorder is smaller than spin splitting energies and the spin splitting bands
are well resolved in the spectral weight as shown in Fig. 3(b). The
corresponding relaxation time is shown in Fig. 3(c), and we find the
relaxation time to be uniform for each Fermi surface, but its values for the
bands 1 and 2 are slightly different. Another difference from Rashba 2DEG is
the Fermi surfaces of band 1 and 2 are hole pockets, not the electron pockets.
When electrical fields are applied in the $+x$ direction, rather than less
electrons with $k_{x}>0$ and $s_{y}<0$ for electron pockets, there are more
electron with $k_{x}>0$ and $s_{y}<0$ and less electrons with $k_{x}<0$ and
$s_{y}>0$. The contribution from this Fermi pocket to CISP is negative (blue
color) as shown in Fig. 3(d)(e). The red Fermi surface has the opposite spin
texture and contribute positively (red color) as shown in Fig. 3(d)(e). As
Rashba 2DEG, the contribution from the blue Fermi pockets is canceled partly
by the red one and make a total negative contribution as shown in Fig. 4(a) at
$E_{f}=1eV$.

Under $\mathcal{T}$ symmetry, the contribution to CISP in Eq. (13) becomes

$\begin{split}\chi_{ij}({\textbf{{k}}}n)&=-\frac{1}{2\pi}\langle{\textbf{{k}}}n|S_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n\rangle\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E_{f})}\\\
&=-\frac{1}{2\pi}\langle-{\textbf{{k}}}n|(\mathcal{T}S_{i}\mathcal{T}^{\dagger})|-{\textbf{{k}}}n\rangle^{*}\langle-{\textbf{{k}}}n|(\mathcal{T}\Gamma_{j}({\textbf{{k}}},E_{f})\mathcal{T}^{\dagger})|-{\textbf{{k}}}n\rangle^{*}\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}({\textbf{{k}}}n,E_{f})}\\\
&=-\frac{1}{2\pi}\langle-{\textbf{{k}}}n|(-S_{i})|-{\textbf{{k}}}n\rangle^{*}\langle-{\textbf{{k}}}n|(-\Gamma_{j}({\textbf{{k}}},E_{f}))|-{\textbf{{k}}}n\rangle^{*}\delta(E_{f}-\epsilon_{-{\textbf{{k}}}n})\frac{\pi}{-\mathop{\mathrm{Im}}\Sigma^{R}(-{\textbf{{k}}}n,E_{f})}\\\
&=\chi_{ij}(-{\textbf{{k}}}n).\end{split}$ (45)

This means the contribution to CISP should be the same at k and
$-{\textbf{{k}}}$, which is clearly shown by the numerical calculations in
Fig. 3(d)(e) and Fig. 4(d)-(g).

We also find that the vertex correction increases the matrix elements of the
spin operator, similar to the case of Rashba 2DEG. Fig. 3(f) shows the
corrected spin operator in Eq. (17) is proportional to and larger than the
bare spin operator and share similar profiles over $\theta$.

The numerical calculations of the CISP from Eq. (12) and Eq. (1) are shown by
solid (dashed) lines and circles (triangles) in Fig. 4(a) for the cases with
(without) vertex corrections for different $E_{F}$. Interestingly, we notice a
sign change of the CISP when lowering $E_{F}$ from $1eV$ to $0eV$ and a rapid
change around $E_{F}=-0.1eV$. The sign change originates from the different
relaxation time for the bands 1 and 2. Fig. 4(c) shows the relaxation time on
four typical points $k_{1}$ to $k_{4}$ for different $E_{F}$. The four points
are intersection points between Fermi surfaces of band 1 and 2 with the high
symmetry lines $\Gamma-K$ and $\Gamma-M$. Fig. 4(b) shows the contribution to
CISP at these points. Around $E_{F}=1eV$, the relaxation time for all the
momenta in Fig. 4(c) has a similar magnitude, which indicates a similar
relaxation time for the bands 1 and 2. Correspondingly, the sign of the CISP
is determined by the Fermi surface size of the bands 1 and 2, and since the
band 2 has a large Fermi surface size, the CISP coefficient $\chi_{yx}$ has a
negative sign. When lowering $E_{F}$ to $0eV$, the relaxation time of the band
2 (blue lines) decreases more rapidly than that of the band 1 (red lines) in
Fig. 4(c), which leads to a dominant contribution from the band 1 for the CISP
and thus, $\chi_{yx}$ changes its sign to be positive around $E_{F}=0eV$.

Besides the sign change, we also notice a rapid change in CISP coefficient
$\chi_{yx}$ between $E_{F}=0eV$ and $E_{F}=-0.2eV$, which is closely related
to the spin anti-vortex around $E_{F}=-0.1eV$. When lowering $E_{F}=0eV$ to
$E_{F}=-0.2eV$, the Fermi pocket $\alpha_{2}$ around the $K$ point comes
across spin anti-vortex at $k_{S}$. The spin texture at two sides of the spin
anti-vortex are opposite along the $\Gamma-K$ line as shown in Fig. 2(a).
Furthermore, the $\alpha_{2}$ Fermi pockets in Fig. Fig. 4(g) have a
relatively large relaxation time. This means the contribution of CISP from the
Fermi pockets $\alpha_{2}$ of the band 2 has the same sign as that from the
band 1, thus leading to a rapid increasing of the CISP around $E_{F}=-0.1eV$
in Fig. 4(a).

The rapid change in CISP shows a special feature of spin anti-vortex that it
is located at non-high-symmetry momenta when compared to the spin vortex
texture located at time reversal invariant momenta in the Rashba model. This
is explained as follows. First, we can consider a time reversal invariant
momentum. Such a momentum is always enclosed by a Fermi pocket since opposite
momenta around a time reversal invariant momentum have to be both occupied or
unoccupied; e.g. $\Gamma$ is always enclosed by the Fermi surface $\beta$
independent of $E_{\text{F}}$, unless the whole Fermi surface completely
disappears. Next, we consider a spin anti-vortex that is located at a non-
high-symmetry momentum. Consequently, it allows the Fermi surface to pass
through the anti-vortex core as $E_{\text{F}}$ varies: e.g. as $E_{\text{F}}$
lowers in Fig. 4(e), the Fermi surface $\alpha$ (red) passes through $k_{S}$
and evolves into the Fermi surface $\alpha_{2}$ (green). The Fermi surfaces
before and after passing through the spin anti-vortex core have opposite spin
textures and thus the contribution from the spin anti-vortex to the CISP can
change its sign rapidly. Fig. 4(e) and (g), where $E_{\text{F}}$ is lowered
and the CISP contribution from the Fermi surface near $k_{S}$ changes from
negative (blue, from the Fermi surface $\alpha$) to positive (red, from the
Fermi surface $\alpha_{2}$). In comparison, the spin polarization at the Fermi
pocket $\beta$ around $\Gamma$ in Fig. 4(e) and (g) keeps the positive sign
(red) and can only vary relatively slowly.

###  II.4 Numerical SHE for monolayer Pb

Figure 5:  (a) SHE. The solid (dashed) line is calculated from the Brillouin
zone formula Eq. (2) with (without) vertex correction. The red circles (blue
triangles) are calculated from the Fermi surface formula Eq. (23). (b)(c) The
diagonal contribution to SHE over Fermi surfaces at $E_{F}=2eV$ without and
with vertex corrections. (d) The diagonal contribution to SHE for different
$E_{F}$. (e)(f) The off diagonal contribution to SHE over Fermi surfaces at
$E_{F}=2eV$ without and with vertex corrections. (g)
$\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|J_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle)$
over Fermi surfaces of $E_{F}=2eV$. (h)
$\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle)$
including vertex corrections. (i)(j)(k)(l) the off diagonal contribution to
SHE over the Fermi surfaces for $E_{F}=1eV,0eV,-0.1eV,-0.2eV$ with vertex
corrections.

In this section, we discuss the numerical calculations of the SHE coefficient
of monolayer Pb shown in Fig. 5(a). For $\sigma^{z}_{yx}$ at high $E_{F}$, we
check our calculation and compare with Rashba 2DEG. When lowering $E_{F}$ from
$1eV$ to $0eV$, there is a sign change of $\sigma^{z}_{yx}$ and for $E_{F}$
between $0.1eV$ and $-0.2eV$, $\sigma^{z}_{yx}$ are changing rapidly, which
are explained next.

To validate our results, we first examine $\mathcal{T}$ symmetry property of
the numerical results. Fig. 5(b) and (c) show the diagonal contribution to the
SHE (Eq. (20)) at Fermi surfaces, which is opposite at k and
$-{\textbf{{k}}}$, as required by $\mathcal{T}$. When summed over Fermi
surfaces, the diagonal contribution shown in Fig. 5(d) is negligible as Eq.
(20) expected. It can be made closer to zero by taking more points on Fermi
surfaces. For off-diagonal contribution ($n^{\prime}\neq n$ in the equation
below), $\mathcal{T}$ symmetry gives

$\begin{split}\sigma^{k}_{ij}({\textbf{{k}}}n)&=-\frac{1}{2\pi}\sum_{n^{\prime}\neq
n}2\pi\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{1}{\epsilon_{{\textbf{{k}}}n}-\epsilon_{{\textbf{{k}}}n^{\prime}}}\mathop{\mathrm{Im}}(\langle-{\textbf{{k}}}n^{\prime}|(\mathcal{T}J^{k}_{i}\mathcal{T}^{\dagger})|-{\textbf{{k}}}n\rangle^{*}\langle-{\textbf{{k}}}n|(\mathcal{T}\Gamma_{j}({\textbf{{k}}},E_{f})\mathcal{T}^{\dagger})|-{\textbf{{k}}}n^{\prime}\rangle^{*}).\\\
&=-\frac{1}{2\pi}\sum_{n^{\prime}\neq
n}2\pi\delta(E_{f}-\epsilon_{{\textbf{{k}}}n})\frac{1}{\epsilon_{{\textbf{{k}}}n}-\epsilon_{{\textbf{{k}}}n^{\prime}}}(-)\mathop{\mathrm{Im}}(\langle-{\textbf{{k}}}n^{\prime}|(J^{k}_{i})|-{\textbf{{k}}}n\rangle\langle-{\textbf{{k}}}n|(-\Gamma_{j}({\textbf{{k}}},E_{f}))|-{\textbf{{k}}}n^{\prime}\rangle)\\\
&=-\frac{1}{2\pi}\sum_{n^{\prime}\neq
n}2\pi\delta(E_{f}-\epsilon_{-{\textbf{{k}}}n})\frac{1}{\epsilon_{-{\textbf{{k}}}n}-\epsilon_{-{\textbf{{k}}}n^{\prime}}}(-)\mathop{\mathrm{Im}}(\langle-{\textbf{{k}}}n^{\prime}|(J^{k}_{i})|-{\textbf{{k}}}n\rangle\langle-{\textbf{{k}}}n|(-\Gamma_{j}({\textbf{{k}}},E_{f}))|-{\textbf{{k}}}n^{\prime}\rangle)\\\
&=\sigma^{k}_{ij}(-{\textbf{{k}}}n).\end{split}$ (46)

So the off diagonal contribution should share the same color at opposite
momenta, consistent with the calculated one shown in Fig. 5(e)(f).

The cancellationg between intrinsic and intrinsic-side-jump contribution in
Rashba 2DEG Sinova  _et al._ (2015) does not happen for monolayer Pb because
the band 1 and 2 at large $E_{F}$ is quantatively different from the Rashba
model. Fig. 5(g) and (h) show
$\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle)$
in Eq. (23) without and with vertex correction. In contrast, for Rashba 2DEG,
it is known that

$\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|J_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle)=\frac{vk_{y}^{2}}{2k}$
(47)

without vertex corrections and $0$ including vertex corrections. Thus, it is
clear that this quantity
$\mathop{\mathrm{Im}}(\langle{\textbf{{k}}}n^{\prime}|J^{k}_{i}|{\textbf{{k}}}n\rangle\langle{\textbf{{k}}}n|\Gamma_{j}({\textbf{{k}}},E_{f})|{\textbf{{k}}}n^{\prime}\rangle)$
shows different behaviors for monolayer Pb and Rashba model.

We next consider the sign change in Fig. 5(a). It can be understood from the
disappearance of the two electron Fermi pockets with positive contribution
(red color) in Fig. 5(i). At $E_{F}=1eV$, the electron Fermi pockets dominate,
giving positive contribution. At $E_{F}=0eV$, the electron Fermi pockets
disappear and the hole Fermi pocktets of band 1 and 2 gives negative
contribution. So when lowering $E_{F}$ from $1eV$ to $0eV$, there is a sign
change of $\sigma_{yx}^{z}$.

The rapid change is related to the Lifshitz transition of Fermi surfaces of
band 2 happened at $k_{L}$ around $E_{F}=-0.1eV$. At $E_{F}=-0.1eV$, the
contribution on Fermi surfaces around $k_{L}$ switches from negative (blue) to
positive (red) values as shown in Fig. 5(j)(k)(l). The contribution remains
negative around spin anti-vortex. Because contribution around $k_{L}$ change
its sign in a narrow energy range, the total contribution has a large
variation.

###  II.5 Estimate of the experimental observables

Figure 6:  (a)(b)(c)(d) Density of states, carrier density, average of spin
splitting energy, average of Fermi velocity with band 1 and 2 for different
$E_{F}$ (e)(f)(g)(h) comparison between numerical results of monolayer Pb and
Rashba 2DEG in longitudinal conductivity, CISP, SHE, spin hall angle over
$E_{F}$. The solid line is the numerical response functions with vertex
corrections. The black circles are calculated with the Rashba 2DEG model.

In this section, we will give a numerical estimate of the magnitude of
$\chi_{yx}$ and $\sigma^{z}_{yx}$ and compare the values with those estimated
from Rashba 2DEG. In the estimate below, we choose the parameter combination
$n_{i}V_{0}^{2}=0.05eV^{2}\text{\AA}^{-2}$ to characterize the impurity
density and strength.

To compare with Rashba 2DEG, the mobility $\mu$ is extracted from the
longitudinal conductivity $\sigma_{xx}=ne\mu$.
$n=\sum_{n{\textbf{{k}}}}\theta(E_{F}-\epsilon_{{\textbf{{k}}}n})$ is the
carrier density. The effective mass is evaluated from the density of states
$D=m/\pi\hbar^{2}$ and the Rashba SOC strength is evaluated from spin
splitting energy $\Delta\epsilon/2k_{F}$ averaged over Fermi surfaces of band
1 and 2.
$\Delta\epsilon=\epsilon_{{\textbf{{k}}}_{F},4}-\epsilon_{{\textbf{{k}}}_{F},3}$
and $k_{F}$ is combined with $m$ for the averaged Fermi velocity
$v_{F}=k_{F}/m$.

The comparisons between monolayer Pb and Rashbs 2DEG are shown in Fig. 6 and
they share same order of magnitude. The mobility is fitted to be $\sim
50cm^{2}/Vs$ from $\sigma_{xx}$. Other response quantities for Rashba 2DEG are
calculated in the following way

$\chi_{yx}=-\frac{1}{2\pi}\tau\alpha_{R}\frac{me}{\hbar}=-\frac{1}{4}D\mu\frac{\Delta\epsilon}{v_{F}}$
(48) $\sigma^{z}_{yx}=\frac{e}{8\pi}$ (49)

The spin hall angle is

$\frac{e}{\hbar}\frac{\sigma^{z}_{yx}}{\sigma_{xx}}=\frac{e}{\hbar}\frac{1}{8\pi
n\mu}$ (50)

## References

  * Mahan (2013) G. D. Mahan,  _Many-particle physics_ (Springer Science & Business Media, 2013).
  * Briggs  _et al._ (2020) N. Briggs, B. Bersch, Y. Wang, J. Jiang, R. J. Koch, N. Nayir, K. Wang, M. Kolmer, W. Ko, A. D. L. F. Duran, _et al._ , Nature materials 19, 637 (2020).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Perdew  _et al._ (1997) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 78, 1396 (1997).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999), arXiv:0927-0256(96)00008 [10.1016] .
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Grimme  _et al._ (2010) S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J. Chem. Phys. 132, 154104 (2010).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. Fabris, G. Fratesi, S. de Gironcoli, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys. Condens. Matter 21, 395502 (2009), arXiv:0906.2569 .
  * Marzari and Vanderbilt (1997) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
  * Mostofi  _et al._ (2008) A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
  * Souza  _et al._ (2001) I. Souza, N. Marzari, and D. Vanderbilt, Phys. Rev. B 65, 035109 (2001).
  * Wei  _et al._ (2021) Y.-W. Wei, C.-K. Li, Y. Cao, and J. Feng, Computer Physics Communications 258, 107551 (2021).
  * Barata and Hussein (2012) J. C. A. Barata and M. S. Hussein, Brazilian Journal of Physics 42, 146 (2012).
  * Inoue  _et al._ (2003) J.-i. Inoue, G. E. Bauer, and L. W. Molenkamp, Physical Review B 67, 033104 (2003).
  * Inoue  _et al._ (2004) J.-i. Inoue, G. E. Bauer, and L. W. Molenkamp, Physical Review B 70, 041303 (2004).
  * Sinova  _et al._ (2015) J. Sinova, S. O. Valenzuela, J. Wunderlich, C. Back, and T. Jungwirth, Reviews of Modern Physics 87, 1213 (2015).