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

# SUPPLEMENTARY INFORMATION
for
”Two-dimensional non-Hermitian skin effect in an ultracold Fermi gas”

Entong Zhao  Department of Physics, The Hong Kong University of Science and
Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Zhiyuan Wang  International Center
for Quantum Materials, School of Physics, Peking University, Beijing 100871,
China  Chengdong He  Department of Physics, The Hong Kong University of
Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Ting Fung Jeffrey Poon
International Center for Quantum Materials, School of Physics, Peking
University, Beijing 100871, China  Ka Kwan Pak  Department of Physics, The
Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Yu-Jun Liu  Department of Physics,
The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Peng Ren  Department of Physics,
The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Xiong-Jun Liu
<EMAIL_ADDRESS>International Center for Quantum Materials, School of
Physics, Peking University, Beijing 100871, China  Hefei National Laboratory,
Hefei 230088, China  International Quantum Academy, Shenzhen 518048, China
Gyu-Boong Jo<EMAIL_ADDRESS>Department of Physics, The Hong Kong University of
Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  IAS Center for Quantum
Technologies, The Hong Kong University of Science and Technology, Hong Kong,
China

##  I Experimental setup and procedure

###  I.1 Experimental setup

The spin-dependent 1D optical lattice potential along $y$ direction is
generated by retro-reflecting a Gaussian laser beam with the frequency
$\omega$ blue detuned by 1GHz from the
${}^{1}S_{0}(F=\frac{5}{2})\rightarrow{{}^{3}}P_{1}(F^{\prime}=\frac{7}{2})$
transition of ${}^{173}\rm{Yb}$ atoms and the polarization linearly polarized
along $x$ direction (Figure S1). The electromagnetic field of the lattice beam
can be written as $\mathbf{E_{y}}\propto\hat{e}_{x}E_{y}\cos(k_{0}y+\omega
t+\phi)=\mathbf{E_{+}+E_{-}}$, with
$\mathbf{E}_{\pm}=\frac{1}{\sqrt{2}}\hat{e}_{\pm}E_{y}\cos(k_{0}y+\omega
t+\phi)$ and $\hat{e}_{\pm}=\frac{1}{\sqrt{2}}(\hat{e}_{x}\pm i\hat{e}_{y})$,
forming the lattice potential

$\displaystyle
V^{\text{latt}}_{\sigma}(y)=\sum_{F^{\prime}}\left(\frac{\hbar\Omega^{2}_{y,\sigma,+,F^{\prime}}}{4\Delta_{y,F^{\prime}}}+\frac{\hbar\Omega^{2}_{y,\sigma,-,F^{\prime}}}{4\Delta_{y,F^{\prime}}}\right)=V_{\sigma}\cos^{2}(k_{0}y),$
(1)

where the effective Rabi frequency $\Omega_{y,\sigma,\pm,F^{\prime}}=\langle
F,\sigma|e\frac{x+y}{\sqrt{2}}|F^{\prime},m_{F}\pm 1\rangle E_{\pm}$ and
single-photon detunings $\Delta_{y,F^{\prime}}$ are determined from all
relevant transitions to the excited states
$(F^{\prime}=\frac{7}{2},\frac{5}{2},\frac{3}{2})$ in the ${}^{3}P_{1}$
manifold.

Another linearly ($\pi$) polarized light with polarization direction $z$ and
frequency $\omega^{\prime}=\omega+\delta\omega$, denoted as
$\mathbf{E_{x}}\propto\hat{e}_{z}E_{x}e^{-ik_{1}y+ik_{2}x+\omega^{\prime}t+\phi^{\prime}}$,
is applied along the direction with angle $76^{o}$ tilted from the lattice
incident beam and induces the Raman coupling between
$|{\uparrow}\rangle=|m_{F}=\frac{5}{2}\rangle$ and
$|{\downarrow}\rangle=|m_{F}=\frac{3}{2}\rangle$, where
$k_{1}=k_{0}\cos\theta$ and $k_{2}=k_{0}\sin\theta$. The frequency difference
$\delta\omega$ between the lattice lights and Raman light can be tuned by the
acousto-optical modulators. The Raman coupling potential $M_{R}$ has the form

$\displaystyle
M_{R}(x,y)=\sum_{F^{\prime}}\frac{\hbar\Omega_{x,\sigma,+,F^{\prime}}\Omega_{R,\sigma,F^{\prime}}}{4\Delta_{y,F^{\prime}}}=M_{0}\cos(k_{0}y)e^{-ik_{1}y+ik_{2}x},$
(2)

with effective Rabi frequency $\Omega_{R,\sigma,F^{\prime}}=\langle
F,\sigma|ez|F^{\prime},m_{F}\rangle E_{x}$.

Figure S 1: a, Schematic energy level diagram with relevant transitions. Both
the lattice and the Raman beams are blue-detuned from
${}^{1}S_{0}(F=\frac{5}{2})\rightarrow{{}^{3}}P_{1}(F^{\prime}=\frac{7}{2})$
intercombination transition and induce the Raman transition between the
$|{\uparrow}\rangle=|{m_{F}=5/2}\rangle$ and
$|{\downarrow}\rangle=|{m_{F}=3/2}\rangle$ hyperfine states of the
${}^{1}S_{0}$ ground manifold. b, Our optical Raman lattice system, which
consists of a 1D optical lattice along the y direction and another Raman beam
along the direction with angle $76^{o}$ tilted from the lattice incident beam,
generating a periodic Raman potential. A loss beam with $\sigma^{-}$
polarization is further applied along -z direction to realize the non-
Hermiticity. c, Absorption image of $m_{F}$ = 5/2 atoms without quenching the
two photon detuning ($t_{f}$=0ms) after 15ms time-of-flight. d, When the two
photon detuning is quenched from $\delta_{i}$ to $\delta_{f}$, two $m_{F}$ =
3/2 atom clouds are coupled out from $m_{F}$ = 5/2 clouds after holding
$t_{f}$=0.15ms, which shows the resonant spin flipping on particular quasi-
momentum subspace.

To separate out spin-1/2 space from other hyperfine states of the ground
manifold, we use an additional $\sigma^{-}$ polarized 556nm light along the z
direction called lift beam to induce spin-dependent ac stark shift. Notably,
the spin-dependent lattice potential
$V_{\uparrow(\downarrow)}=3.2{E_{r}}(2.2{E_{r}})$ will induce an additional
on-site energy difference $\delta_{0}$ to the Bloch bands of two spin states,
which can be compensated by the two-photon detuning. The quantized axis is
fixed by the bias magnetic field of 13.6G along the z-direction, which will
induce an 8.1MHz Zeeman splitting between adjacent $m_{F}$ states. The loss
beam is $\sigma^{-}$ polarized and red detuned by $\Delta_{\downarrow}=1.0$ or
1.6MHz and $\Delta_{\uparrow}=9.1$ or 9.7MHz with respect to
${}^{1}S_{0}(F=\frac{5}{2},m_{F}=\frac{3}{2})\rightarrow{{}^{3}P_{1}}(F^{\prime}=\frac{7}{2},m_{F^{\prime}}=\frac{1}{2})$
and
${}^{1}S_{0}(F=\frac{5}{2},m_{F}=\frac{5}{2})\rightarrow{{}^{3}P_{1}}(F^{\prime}=\frac{7}{2},m_{F^{\prime}}=\frac{3}{2})$.
Based on CG-coefficients, the estimated ratio of loss rate between two spin
states $\gamma_{\downarrow}/\gamma_{\uparrow}$ is around 110 for
$\Delta_{\downarrow}=-1.0$MHz and 248 for $\Delta_{\downarrow}=-1.6$MHz, which
means the $\gamma_{\uparrow}$ term can be ignored. Remarkably, the loss beam
will also induce an additional energy shift in two photon detuning, which we
can compensate by tuning the frequency difference between Raman beam and
lattice beam in our experiment.

###  I.2 Experimental procedure

Figure S 2: Experimental sequence for momentum resolved Rabi spectroscopy (a)
and spin polarization texture measurement in equilibrium (b).

The experimental sequence for the momentum resolved Rabi spectroscopy can be
briefly described in Figure S2a. A single-component degenerate
${}^{173}\text{Yb}$ Fermi gas in the $|{\uparrow}\rangle$ state is prepared
with the optical pumping process during (with 399nm light) and after (with
556nm light) the evaporation. The atoms are then loaded into a 1D optical
lattice potential by ramping up both the lattice and Raman laser intensity in
3ms with the initial two photon detuning $\delta_{i}=-27E_{r}$ far enough to
exclude the Raman coupling. The final lattice depth and Raman coupling
strength are set to $V_{\uparrow(\downarrow)}=3.2(2.2)E_{r}$ and
$M_{0}=2.0E_{r}$ ($\text{Er}=\hbar^{2}k_{0}^{2}/2m=h\times 3.735\text{kHz}$),
which can be calibrated from Kapitza-Dirac diffraction and the two-photon Rabi
oscillation, respectively. After the initial state preparation, we further
hold the system for another 1ms and then quench the two-photon detuning from
$\delta_{i}$ to $\delta_{f}$ at time $t=0$ by controlling the RF signal of the
acousto-optical modulator and the loss pulse is switched on at the same time.
The atoms initially prepared in $|{\uparrow}\rangle$ state start the Raman-
Rabi oscillation between the lowest two bands in our system. Holding a
variable time $t_{f}$, we record spin-sensitive momentum distribution after
10ms (or 15ms) time of flight expansion by applying the 556nm blast pulses
resonant to
${}^{1}S_{0}(F^{\prime}=5/2)\rightarrow{}^{3}P_{1}(F^{\prime}=7/2)$
transition.

In our Rabi spectroscopy, we found that there is no obvious difference between
the lattice loading in 3ms and 9ms, probably because the final lattice depth
is not quite deep. From the estimation of Kapitza-Dirac diffraction, even
though the lattice is suddenly switched on, the maximum population for higher
order is smaller than 30% and the oscillation period is around $60\mu s$, much
smaller than the hold time (1ms). Therefore, the shorter loading time will not
affect the Rabi spectroscopy and will maximumly suppresses the heating and
other detrimental effects. Besides, it is worth noting that the optical dipole
trap is switched off at the same time as the quench dynamics starts. This is
because, in previous experiments song2018-observation ; ren2022-chiral , we
found the intraband transition induced by optical dipole trap will affect the
time evolution. Since the time evolution is smaller than 0.4ms, the atoms fall
less than 1$\mu m$ under gravity, which is much smaller than the beam waist of
the lattice and Raman beam and can be ignored in our experiment.

###  I.3 Spin polarization of equilibrium state in Hermitian regime

Except for the momentum-resolved Rabi spectroscopy, we also measured the spin
polarization texture in equilibrium with experimental sequence in Figure S2b.
An optical ac Stark shift is added by a lift beam 2ms before the optical
lattice potential, which separates out an effective spin-1/2 subspace from
other hyperfine levels. The Raman beam and lattice beams are adiabatically
switched on with an 8-ms exponential ramp to the final value. The 10ms spin-
resolved time-of-light image is taken after a spin sensitive blast sequence to
reconstruct the spin polarization texture. Dependent on the value of two
photon detuning (equivalent to $k_{x}$), our system can form either a
topological band or trivial band along $q_{y}$ direction, as shown in Figure
S3a. The experimentally measured spin polarization $P(q_{y})$ can be
reconstructed from the integrated momentum distribution along $k_{x}$
direction within $k_{x}/k_{0}=-0.2$ to $0.2$, as shown in Figure S3b, similar
to the detection technique in song2019-observation . When the two photon
detuning is in the trivial regime, the spin polarization is either
$|{\uparrow}\rangle$(red)-dominated or $|{\downarrow}\rangle$(blue)-dominated.
However, when the Raman transition resonantly couples the lowest two bands,
spin polarization changes from blue to red and blue again within the first
Brillouin zone, which indicates the nonzero winding number and topological
band song2018-observation .(Figure S3c) Notably, since the atom number at
higher $k_{x}$ is not enough for a good signal-to-noise ratio, we choose to
fix $k_{x}$ and scan the value of two photon detuning instead of fixing the
two photon detuning and changing $k_{x}$ to obtain the whole spin polarization
texture in equilibrium case.

Figure S 3: Spin polarization of equilibrium state. a Numerical simulation of
spin polarization as a function of two photon detuning $\delta$. b
Experimental measurement of spin polarization as a function of two photon
detuning $\delta$. c Theoretical and experimental spin polarization texture of
the lowest band for $\delta=-1.6,-0.3,1.0E_{r}$, respectively.

###  I.4 Calibration of the atomic loss

In our system, the spin dependent loss is induced by an optical pumping beam
at a small detuning from the narrow intercombination transition
${}^{1}S_{0}(F=\frac{5}{2})\rightarrow{}^{3}P_{1}(F=\frac{7}{2})$, as shown in
Figure S4a. The optical pumping beam is red detuned by
$\Delta_{\downarrow}=1.0$ or 1.6MHz and $\Delta_{\uparrow}=9.1$ or 9.7MHz with
respect to
$|{\downarrow}\rangle\rightarrow|{m_{F^{\prime}}=\frac{1}{2}}\rangle$ and
$|{\uparrow}\rangle\rightarrow|{m_{F^{\prime}}=\frac{3}{2}}\rangle$, together
with the relative optical transition strength or CG-coefficients (Figure S4b),
resulting in the spin dependence of atomic loss. In our experiment, we first
prepare the atoms on specific spin state $|{\sigma}\rangle=|{\uparrow}\rangle$
or $|{\downarrow}\rangle$ in optical lattice and then suddenly switching on
the loss pulse. The two photon detuning $m_{z}(\delta)$ is kept constant and
large enough ($\delta\sim-134E_{r}$) to avoid the Raman coupling during the
pulse, and the relative atom number on specific state can then be resolved by
optical Stern Gerlach measurment. The experimental sequence for other optical
beams are same as that in momentum resolved Rabi spectroscopy measurement to
avoid the effect of other types of dissipation. We further calibrate the loss
rate by fitting the relative atom number of spin states
$|{\sigma=\uparrow,\downarrow}\rangle$ with
$f(t)=Ae^{-\frac{\gamma_{\sigma}}{\hbar}t}+B$ and the loss rate is tunable by
controlling the power of the loss beam (Figure S4c). From the calibration, the
loss rate between two spin states $\gamma_{\downarrow}/\gamma_{\uparrow}$ is
larger than 20 for $\Delta_{\downarrow}=-1.6$MHz.

Figure S 4: Calibration of the atomic loss. a Schematic energy level for spin
dependent loss. b Relative optical transition strength of different $m_{F}$
states for the narrow intercombination transition
${}^{1}S_{0}(F=\frac{5}{2})\rightarrow{}^{3}P_{1}(F=\frac{7}{2})$. c Loss rate
calibration for spin-down and spin-up state at same optical intensity. d Loss
rate calibration for spin-down state at variable intensity of the loss beam.
Figure S 5: 2D momentum dependent time evolution for spin polarization texture
$P(k_{x},q_{y},t)$ in a, normalized spin up texture
$\mathcal{D}_{\uparrow}(k_{x},q_{y},t)/\mathcal{D}_{\uparrow}(k_{x},q_{y},t=0)$
in b, normalized spin down texture
$\mathcal{D}_{\downarrow}(k_{x},q_{y},t)/\mathcal{D}_{\uparrow}(k_{x},q_{y},t=0)$
in c with $\gamma_{\downarrow}=0.0,0.4,1.0,2.1E_{r}$, respectively.  Figure S
6: 2D time averaged spin polarization texture obtained by the evolution from 0
ms to 0.30 ms with $\gamma_{\downarrow}=0.0,0.5,1.0,1.2,1.8,2.1E_{r}$,
respectively.

###  I.5 Momentum-dependent time evolution

To obtain the spin information in the quasi-momentum space, we take a series
of spin resolved absorption images $\\{\mathcal{I}_{1}(\text{no atoms
killed}),\mathcal{I}_{2}(|{\uparrow}\rangle\text{atoms
killed}),\mathcal{I}_{3}(|{\uparrow,\downarrow}\rangle\text{atoms killed})\\}$
after a 10ms or 15ms time-of-flight expansion. The atomic distribution of spin
up and spin down can be extracted from
$\mathcal{I}_{\uparrow}(k^{\prime}_{x},k_{y})=\mathcal{I}_{1}-\mathcal{I}_{2}$
and
$\mathcal{I}_{\downarrow}(k^{\prime}_{x},k_{y})=\mathcal{I}_{2}-\mathcal{I}_{3}$.
Notably, a fraction of atoms may occupy other spin states due to the imperfect
isolation of the spin-$\frac{1}{2}$ subspace or spontaneous emission of the
atoms on excited states pumped by the loss beam, which can be eliminated by
the background image $\mathcal{I}_{3}$. The atomic distributions
$\mathcal{I}_{\sigma}(k^{\prime}_{x},k_{y})$ are then folded into the first
Brillouin zone by shifting the second and third Brillouin zone with integer
number of $2k_{0}$,
$\mathcal{D}_{\sigma}(k^{\prime}_{x},q^{\prime}_{y})=\sum_{n=-1}^{1}\mathcal{I}_{\sigma}(k^{\prime}_{x},q^{\prime}_{y}+2nk_{0})$.
After that, based on the spin-momentum locking, momentum distributions in
quasi-momentum are defined as
$\mathcal{D}_{\uparrow}(k_{x},q_{y})=\mathcal{D}_{\uparrow}(k^{\prime}_{x},q^{\prime}_{y})$
and
$\mathcal{D}_{\downarrow}(k_{x},q_{y})=\mathcal{D}_{\downarrow}(k^{\prime}_{x}+k_{0}\sin{\theta},\text{mod}{(q^{\prime}_{y}+k_{0}-k_{0}\cos{\theta},2k_{0})})$,
where mod is the modulo operator. Finally, the spin polarization can be
calculated by
$P(k_{x},q_{y})=\frac{\mathcal{D}_{\uparrow}(k_{x},q_{y})-\mathcal{D}_{\downarrow}(k_{x},q_{y})}{\mathcal{D}_{\uparrow}(k_{x},q_{y})+\mathcal{D}_{\downarrow}(k_{x},q_{y})}$.
The time evolution of spin polarization texture, normalized spin up texture
and spin down texture are shown in Figure S5. In the Hermitian regime with
$\gamma_{\downarrow}=0E_{r}$, the spin polarization at different quasi
momentum oscillates with different frequency following Rabi process after
quenching the two photon detuning and a line structure appears at specific
time (0.08$\sim$0.20ms), which reveals a non-uniform band gap. When the spin-
dependent atomic loss is induced, the oscillation frequency at quasi-momentum
near the line structure become slower and the oscillation amplitude becomes
smaller, which indicates the band gap closing behavior. The line structure,
where the Raman coupling is resonant, can also be found in normalized spin up
and spin down texture as shown in Figure S5b-c. In particular, when the
dissipation is comparable with the Raman coupling, since the spin up atoms
near the line structure are constantly transfered to the spin down states by
the two photon Raman process and the spin down atoms experienced a large
atomic loss, there are few atoms return spin up state, which makes the line
structure in spin up texture (Figure S5b) can be observed for longer duration
than Hermitian case. We can further obtain the time averaged spin texture by
averaging the data with evolution time from 0 ms to 0.30 ms after quenching
$m_{z}(\delta)$ as shown in Figure S6, where the zero time-averaged spin
polarizations indicate the band inversion surface in Hermitian regime
zhang2019-dynamical ; yi2019-observing ; song2019-observation .

Figure S 7: a, Time evolution of spin polarization for different loss rate at
$k_{x}=0$ layer. b, The 2D band gap reconstructed from the spin polarization
texture at different dissipation strength. When the loss rate is closed to the
exceptional point, the significantly reduced atom number makes it hard to
obtain reasonable fitting result for some quasi-momentum, which is marked by
the black cross.

###  I.6 Extraction of band gap information

In Figure S7a, the time dependent spin texture at $k_{x}=0$ for different
dissipation strength is shown. From the time-dependent spin polarization
curve, one can roughly extract the band gap information with a empirical
damped sinusoidal function $f(t)=Ae^{-Bt}\cos{Ct}+1-A$, where B corresponds to
the imaginary part of energy gap and C corresponds to the real part of energy
band. The extracted 2D band gap are shown in Figure S7b, which shows similar
ring-like pattern as simulation from plane wave expansion in Fig.2a.

However, the phase diagram obtained from this empirical method deviates from
the theoretical simulation, as shown in the inset of Fig.3a. Even prior to the
exceptional point, a finite damping rate exists, and the oscillation frequency
obtained from fitting cannot be completely zero or may have a large fitting
error after exceptional point. There are three primary reasons why the phase
diagram deviates from the simulation. Firstly, in addition to the uncertainty
of atom loss caused by the power fluctuation of the loss beam (green shaded
region in main Fig.3a), there is also momentum uncertainty induced by the
finite resolution of the imaging system and the fluctuation of two-photon
detuning in the system. This uncertainty causes the oscillation curve to
dephase even without spin-dependent dissipation. In main Fig.3, the effect of
this type of momentum uncertainty on the phase diagram is labeled by the brown
shaded region and explain the reason that the imaginary band gaps obtained by
empirical damped sinusoidal function appear before the exceptional points.
Secondly, when the dissipation is comparable to exceptional points, the atom
number is significantly reduced, and the oscillation curve becomes a monotonic
decay curve dominated by dissipation, making it difficult to obtain the
oscillation term and further determine the position of the exceptional point.
Finally, the theoretical time evolution is not a perfectly damped oscillation,
which also contributes to the deviation.

Therefore, we also quantify the different behaviors before and after PT-
symmetry breaking transition by fitting with both damped oscillation
$f(t)=Ae^{-Bt}\cos{Ct}+1-A$ and exponential decay $f(t)=Ae^{-Bt}+1-A$, and
then compare the
$R^{2}=1-\frac{\sum_{i=1}^{n}(y_{i}-f(x_{i}))^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}}$,
like shown in main Fig.3b. When the loss rate is large enough, the $R^{2}$ of
both fits overlap with each other, which indicates the needless of the
oscillation term when fitting the time evolution at large loss. With this
method, we can estimate where the transition occurs and then fit the time
evolution curve before and after transition region with sinusoidal function
and exponential function respectively to extract the energy band and damping
rate (Fig 3a). For different $k_{x}/k_{0}$ (or equivalently $\delta$), the
loss rate of EPs are different.

###  I.7 Relation between two photon detuning and $k_{x}$

From the Hamiltonian of our Raman lattice, one can see that although the
lattice potential is one-dimensional in the $y$ direction, the running-wave
term $e^{ik_{2}x}$ in the Raman potential couples the spin-up (spin-down)
states with the x-directional kinetic energy
$\frac{\hbar^{2}}{2M}(k_{x}\pm\frac{k_{2}}{2})^{2}$. This coupling gives rise
to a linear SO term $\frac{\hbar^{2}}{2M}(k_{x}k_{2})\sigma_{z}$, which leads
to a 2D energy band. In fact, the linear SO term effectively shifts the Zeeman
energy $m_{z}$, making the energy band for the $k_{x}=0$ layer obtained by
tuning $m_{z}(\delta)$ identical to that for fixed $m_{z}(\delta)=m_{0}$ by
scanning $k_{x}$ following the relation
$H_{q_{y},k_{x}=0,m_{z}=m_{0}}=H_{q_{y},k_{x}=\frac{2Mm_{0}}{\hbar^{2}k_{2}},m_{z}=0}$.

In main Fig. 3a and b, we chose different detuning values and then transferred
them into different $k_{x}$ values due to the limited atom number at high
$k_{x}$. In Figure S8, we compared the reconstructed band gap by shifting the
two-photon detuning values $\delta$ and $k_{x}$ in Hermitian regime, which
confirmed the relationship between the two-photon detuning and $k_{x}$ in our
experiments.

Figure S 8: Comparison between reconstructed band gaps for $\delta=-0.3E_{r}$
(inset on the right) and $\delta=-1.4E_{r}$(inset on the left). For
$\delta=-1.4E_{r}$, the data is further shifted along $k_{x}$ direction during
the analysis to compensate the difference of two photon detuning. The error
bars stand for the fitting error.

###  I.8 Momentum uncertainty and reconstruction of band difference on complex
plane

As mentioned in Section F, the momentum uncertainty makes it hard to
distinguish dephasing issue from the decay induced by atom dissipation when
fitting the oscillation of spin polarization with damped oscillation.
Therefore, to avoid the effect of dephasing issue, when we reconstruct the
complex eigen spectrum in low dissipation regime, the imaginary gap is
determined by fitting the oscillation curve of spin down atoms by following
equation

$f(t)=A(e^{2\Im\left(E_{2}\right)t}+e^{2\Im\left(E_{1}\right)t})-2B{e}^{\Im\left({E_{1}+E}_{2}\right)t}\cos{\Re\left(E_{2}-E_{1}\right)t}+2B-2A$
(3)

This equation is obtained by combining the time-dependent Schrödinger equation
of tight-binding Hamiltonian and an empirical damped sinusoidal function.
Here, we use the parameters obtained from fitting the spin polarization with
damped sinusoidal function as the initial fitting parameters. This fitting
function enables us to distinguish the decay induced by the imaginary band
from the dephasing induced by the momentum uncertainty. By combining the real
band gap extracted from the oscillation curve of the spin polarization, we can
further reconstruct the band difference on the complex plane. Figure S9
displays the time evolution of spin polarization and normalized spin down for
selected points on the complex eigen spectrum with
$\gamma_{\downarrow}=0.0E_{r}$ and $\gamma_{\downarrow}=0.4E_{r}$,
respectively, along with the fitting curve. When there is no spin-dependent
dissipation, the normalized spin down curve with only dephasing induced by
momentum uncertainty has a non-zero offset after a sufficiently long time.
This behavior can be described by the above fitting function with
$\Im\left(E_{1}\right)=\Im\left(E_{2}\right)$ and $A=0$, which is an empirical
damped sinusoidal function. Conversely, when $\gamma_{\downarrow}\neq 0$, the
spin down atoms eventually reach zero, which implies that
$\Im\left(E_{1}\right)\neq\Im\left(E_{2}\right)$ and $A=B$ in the above
equation, which is the exact solution of time-dependent Schrödinger equation
of tight-binding Hamiltonian.

Figure S 9: Reconstructing the band difference on the complex plane. a, Time
evolution of spin polarization and normalized spin down for selected point on
the complex eigen spectrum with $\gamma_{\downarrow}=0.0E_{r}$. b, Time
evolution of spin polarization and normalized spin down for selected point on
the complex eigen spectrum with $\gamma_{\downarrow}=0.4E_{r}$. In both
figures, the dots are experiment data and the solid line is the fitting
result.

###  I.9 Band difference on complex plane

Figure S10 presents the band difference on the complex plane at various
$k_{x}$ layers ($k_{x}/k_{0}=-0.4,-0.1,0.1,0.4$) with dissipation
$\gamma_{\downarrow}=0.4E_{r}$. For all the $k_{x}$ layers depicted in Figure
S10, the band difference exhibits a non-zero spectral winding number on the
complex plane. This non-backstepping curve signifies that the system maintains
a universal non-Hermitian skin effect (NHSE) along the $y$ direction, as
corroborated in section II.B. It is also important to highlight that only the
absolute value of the imaginary band gap can be measured in the experiment.
According to the simulation, the original band difference spectrum on the
complex plane manifests a ring-like pattern, and the imaginary part of the
band gap alters its sign when the quasi-momentum traverses the band inversion
surface. By flipping the portion with a negative imaginary band gap into the
first quadrant, the ring-like pattern transforms into the butterfly shape
depicted in the figures. In the case along the $q_{y}$ direction, as no
optical lattice exists in this direction, the band difference spectrum
intersects at $k_{x}\rightarrow\pm\infty$. Since only a portion of the band
difference spectrum is measured in the experiment, it does not display a non-
zero spectral winding. However, for specific $k_{x}$ layers
($q_{y}/k_{0}=-0.2,0.6$), the currently measured region is sufficiently
expansive to observe the non-backstepping behavior. Consequently, the system
also upholds a universal NHSE along the $x$ direction, and based on these two
points, the entire 2D system must exhibit a non-zero spectral area.

Figure S 10: Band difference at different $k_{x}$ layers on complex plane with
$\gamma_{\downarrow}=0.4E_{r}$. The top row is the experimental data, while
the bottom row is the corresponding simulation result. Figure S 11: Band
difference at different $q_{y}$ layers on complex plane with
$\gamma_{\downarrow}=0.4E_{r}$. The top row is the experimental data, while
the bottom row is the corresponding simulation result.

##  II Spectral topology and non-Hermitian skin effect

In this section, we focus on the complex spectral topology and the NHSE. In
the first subsection we review the theorems for skin effect in lattice systems
and generalize the theory to continuous space. In the second subsection, based
on the band difference obtained in the experiment we get access to non-trivial
spectral topology and further predict the NHSE under open boundary condition
(OBC). In the third subsection we numerically solve a 2D lattice model and
discuss the features of the skin effect from a symmetry point of view.

###  II.1 The criteria for non-Hermitian skin effect in 1D and 2D systems

In this subsection we review the theorems for the occurrence of the skin
effect in (quasi-)1D and 2D lattice system and generalize the theory to the
case of continuous space. The NHSE has a topological origin manifested by the
bulk-boundary correspondence between the spectral topology of the periodic
boundary eigenspectrum and the skin modes in the real space. For (quasi-)1D
lattices, it has been rigorously proved that the emergence of skin effect is
equivalent to the existence of nonzero spectral winding numbers

$w(E_{0})=\int_{\rm{FBZ}}\frac{dk}{2\pi i}\partial_{k}\det(H(k)-E_{0}),$ (4)

where $H(k)$ is the Bloch Hamiltonian with quasi-momentum $k$, and $E_{0}$
could be any gap point sato2020-origin . An equivalent and more intuitive
description is that the periodic boundary spectrum $E_{i}(k)$ ($i$ runs over
all bands) is non-backstepping on the complex plane. Here we make a
generalization that this also applies for the skin effect in the continuous
space. This is because for any given eigenstate of a Hamiltonian in the
continuous space we could treat it as an eigenstate of a lattice model by
applying the finite difference approximation with a proper step length. The
discretized model obtained in this way serves as a good approximation for the
low-energy and long-wavelength spectrum of the continuous model. So if the
spectrum of the continuous model is non-backstepping, that of the lattice
model should also be so and the corresponding open boundary eigenstates of the
lattice model should be skin modes. Therefore the continuous model should also
have skin modes and exhibit NHSE. We try this criterion on two continuous two-
band models which in the momentum space are written as:

$\displaystyle H_{1}(k)$
$\displaystyle=\begin{pmatrix}\hbar^{2}(k-Q)^{2}/2m&\Omega_{R}\\\
\Omega_{R}^{*}&\hbar^{2}(k+Q)^{2}/2m-i\Gamma\end{pmatrix},$ $\displaystyle
H_{2}(k)$
$\displaystyle=\begin{pmatrix}\hbar^{2}(k-Q)^{2}/2m&\Omega_{R}+\Gamma/2\\\
\Omega_{R}^{*}-\Gamma/2&\hbar^{2}(k+Q)^{2}/2m\end{pmatrix},$ (5)

where $Q,\Omega_{R}$ and $\gamma$ are all constants. If $Q$ takes the value
$k_{0}\sin 76^{o}/2$ and $\Gamma=\gamma_{\downarrow}/2$ is the loss rate of
spin down atoms, $H_{1}$ is just the quasi-1D layer of our experimental model.
We could see in Figure SLABEL:FigS911(a1,a2) that $H_{1}$ has non-backstepping
periodic boundary spectrum and exhibits skin effect under the OBC, while
$H_{2}$ has backstepping spectrum and no skin effect.

In a 2D periodic system, a nonzero periodic boundary spectral area on the
complex plane results in the universal skin effect, which means that the skin
effect persists unless the open boundary geometry coincides with some spatial
symmetry of the bulk zhang2022-universal . And more specifically in a 2D
lattice with a rectangular open boundary, the skin effect can be absent only
if the spectral winding numbers of quasi-1D layers perpendicular to the
boundary

$w_{E}(\vec{k}_{\|})=\int_{\rm{FBZ}}\frac{dk_{\perp}}{2\pi
i}\partial_{k_{\perp}}\det(H(k_{\perp},\vec{k}_{\|})-E_{0})$ (6)

are zero for all gap points $E_{0}$ and fixed momenta $\vec{k}_{\|}$ parallel
to the boundary. Such winding numbers indicate the direction of spectral flow
and on which boundary the skin modes accumulate in the real space. If the
system is continuous space in one or two directions, the spectral winding
numbers are no longer available. In that case, we just need to check whether
the spectrum of the quasi-1D layers are non-backstepping and the spectral flow
will be clear.

Figure S 12: Complex spectra and non-Hermitian skin effect of the 1D and 2D
models. a1,a2, open boundary spectrum (blue) and periodic boundary spectrum
(red, calculated in the momentum space) of $H_{1/2}$ in Eq. 5. The insets are
the density profile summed over all open boundary eigenstates with eigenenergy
$E$ satisfying $|\Re E|<2\rm{Er}$ and spin degree of freedom $\sum_{i,s,|\Re
E_{i}|<2\rm{Er}}|\Psi_{i}(x,s)|^{2}$. b1, low energy part of the periodic
boundary (red) and open boundary (blue) spectra of the 2D lattice model in Eq.
(8). b2, the density profile summed over open boundary eigenstates of the
model in Eq. 8 with eigenenergy $E$ satisfying $\Re E<1.5\rm{Er}$. b3, density
profile of the 50 eigenstates with eigenenergies most closest to the base
point $E_{0}=(1.1+0.08\text{I})E_{r}$. b4, density profile of the 50
eigenstates with eigenenergies closest to the base point
$-E_{0}=(-1.1-0.08\text{I})E_{r}$.

###  II.2 A proof of the nontrivial spectral topology and the prediction of
skin effect

In this subsection we deduce the nontrivial spectral based on the experimental
results and predict the occurrance of NHSE under OBC. The band difference
$\Delta E(k_{x},q_{y})=E_{1}(k_{x},q_{y})-E_{2}(k_{x},q_{y})$ between the two
lowest bands extracted from the quench dynamics of spin polarization shows
nonzero spectral area and non-backstepping curves of quasi-1D layers in
$k_{x}$ and $q_{y}$ directions. But to predict skin effect we need the actual
band dispersion rather than their difference. To bridge the gap, here we show
that there exists at least one topological nontrivial quasi-1D subsystem in
both $x$ and $y$ directions. We make the following assertion and provide a
proof: If the band difference $\Delta E(k_{x},q_{y})$ of the quasi-1D layers
with fixed $k_{x}$ or $q_{y}$ is non-backstepping on the complex plane and
does not pass the origin point ($\Delta E\neq 0$), the bands
$E_{1/2}(k_{x},q_{y})$ should also be non-backstepping.

We prove the case of quasi-1D layer along $q_{y}$ direction with fixed
$k_{x0}$. First we define
$h(k_{x},q_{y})=(E_{1}(k_{x},q_{y})+E_{2}(k_{x},q_{y}))/2$. With the band
difference $\Delta E(k_{x},q_{y})=E_{1}(k_{x},q_{y})-E_{2}(k_{x},q_{y})$ the
bands could be written as

$\displaystyle E_{1/2}(k_{x},q_{y})=h(k_{x},q_{y})\pm\Delta E(k_{x},q_{y})/2.$
(7)

The bands $E_{1/2}(k_{x0},q_{y})$ could be backstepping only in two cases: (i)
$\forall q_{y}$, $\exists q_{y}\neq q_{y}^{\prime},q_{y}^{\prime\prime}$ s.t.
$E_{1}(k_{x0},q_{y})=E_{1}(k_{x0},q_{y}^{\prime})$,
$E_{2}(k_{x0},q_{y})=E_{2}(k_{x0},q_{y}^{\prime\prime})$, or (ii) $\forall
q_{y}$, $\exists$ only one $q_{y}^{\prime}\neq q_{y}$ s.t.
$E_{1}(k_{x0},q_{y})=E_{2}(k_{x0},q_{y}^{\prime})$. In case (i) there should
be $q_{y}^{\prime}=q_{y}^{\prime\prime}$ because such band degeneracy should
be caused by a certain symmetry transforming the momentum in a physical sense,
like $q_{y}\rightarrow-q_{y}$, $q_{y}\rightarrow\pi-q_{y}$ or so on. Then we
equivalently have $h(k_{x0},q_{y})\pm\Delta
E(k_{x0},q_{y})=h(k_{x0},q_{y}^{\prime})\pm\Delta E(k_{x0},q_{y}^{\prime})$
and thus $\Delta E(k_{x0},q_{y})=\Delta E(k_{x0},q_{y}^{\prime})$. This
$q_{y}\rightarrow q_{y}^{\prime}$ relation actually gives a local
homeomorphism $q_{y}^{\prime}(q_{y})$ from the first Brillouin zone (FBZ) to
itself. If the bands $E_{1/2}(k_{x0},q_{y})$ are backstepping,
$q_{y}^{\prime}(q_{y})$ and $q_{y}$ should go in the opposite directions
($q_{y}\rightarrow q_{y}^{\prime}$ have a negative mapping degree, rigorously
speaking). Therefore $\Delta E(k_{x0},q_{y})=\Delta E(k_{x0},q_{y}^{\prime})$
indicates that the $\Delta E(k_{x0},q_{y})$ curve is backstepping. This
contradiction excludes the possibility of the first case. In case (ii) we have
$h(k_{x0},q_{y})+\Delta E(k_{x0},q_{y})=h(k_{x0},q_{y}^{\prime})-\Delta
E(k_{x0},q_{y}^{\prime})$. Since $q_{y}$ and $q_{y}^{\prime}$ run in the
opposite directions, there must be some $q_{y1}$ s.t.
$q_{y1}=q_{y1}^{\prime}$. Then we have $\Delta E(k_{x0},q_{y1})=0$, which
contradicts the condition that the curve of $\Delta E$ does not pass the
origin point. By excluding the possibilities of these two cases, we conclude
that the two bands $E_{1/2}(k_{x},q_{y})=$ themselves are non-backstepping. As
a sidenote, our proof applies to the half-integer winding number case that
across the Brillouin zone the two bands exchange eigenstates rather than
returning to the original eigenstates respectively.

Applying this on the non-backstepping $\Delta E(k_{x},q_{y})$ curves as shown
by Figure S10 and Figure S11, we deduce that the corresponding quasi-1D band
spectra are non-backstepping. Combining this with the lack of reciprocity and
spatial symmetry of our system, we predict the appearance of corner skin modes
under the OBC, which is supported by the numerical simulation shown in Figure
S12b2.

With the non-backstepping spectrum of both quasi-1D layers in $k_{x}$ and
$q_{y}$ directions, we could further infer that the periodic boundary spectrum
of the whole 2D system must have a nonzero spectral area. This is because
otherwise the spectrum on the complex plane should be a 1D curve including
both non-backstepping open curves of $k_{x}$ layers going to real positive
infinity and the non-backstepping loops of $q_{y}$ layers, which is not
possible. The nonzero spectral area also supports the occurrence of 2D NHSE
zhang2022-universal .

###  II.3 Simulation of the non-Hermitian skin effect and the nonlocal feature

To see the NHSE of our optical lattice system, here we solve a 2D lattice
model with open boundary cut along $x$ and $y$ directions. The model is
obtained by applying s-band tight-binding approximation in the $y$ direction
and finite-difference approximation in the $x$ direction. It acts on the local
basis $|\vec{j}s\rangle$ centered at the site $\vec{j}=(j_{x},j_{y})$ with
spin $s=(\uparrow,\downarrow)$, which in the $y$ direction are the s-band
Wannier functions. The Hamiltonian is

$\displaystyle H_{lattice}=$
$\displaystyle\sum_{\vec{j}}\Big{\\{}(m_{z}+i\gamma_{\downarrow}/2)\bigr{(}|\vec{j}\uparrow\rangle\langle\vec{j}\uparrow|-|\vec{j}\downarrow\rangle\langle\vec{j}\downarrow|\bigr{)}+\sum_{k=x,y}\Big{[}\bigr{(}-t_{\uparrow}^{k}|\vec{j}+\vec{e}_{k}\uparrow\rangle\langle\vec{j}\uparrow|+t_{\downarrow}^{k}|\vec{j}+\vec{e}_{k}\downarrow\rangle\langle\vec{j}\downarrow|\bigr{)}+h.c.\Big{]}$
$\displaystyle+\bigr{(}t_{so}e^{ik_{0}\cos\theta
j_{y}}|\vec{j}\downarrow\rangle\langle\vec{j}+\vec{e}_{y}\uparrow|-t_{so}e^{-ik_{0}\cos\theta
j_{y}}|\vec{j}\uparrow\rangle\langle\vec{j}+\vec{e}_{y}\downarrow|+h.c.\bigr{)}\Big{\\}},$
(8)

where $m_{z}$ is the Zeeman energy, $\gamma_{\downarrow}$ is the loss rate of
spin-down atoms, $t_{\uparrow,\downarrow}^{k}$ are the nearest-neighbor spin-
conserved hopping coefficients of the optical lattice in the $k=(x,y)$
directions, $t_{so}$ is the nearest-neighbor spin-flipped hopping coefficient
in the optical lattice direction, and $\theta=76^{o},k_{0}=\pi/a$ with $a$
being the lattice constant in the y direction. Here
$t_{\uparrow}^{x}=t_{\downarrow}^{x}=-i\sin\theta\frac{a}{2l_{x}\pi}+\frac{a^{2}}{l_{x}^{2}\pi^{2}}$
where $l_{x}$ is the finite difference step-length. We could gauge out the
phase factor $e^{\pm ik_{0}\cos\theta j_{y}}$ with a transformation
$|\vec{j}\downarrow\rangle\rightarrow e^{-ik_{0}\cos\theta
j_{y}}|\vec{j}\downarrow\rangle$. Here the tight-binding parameters are set to
be $t_{\uparrow(\downarrow)}=0.32(0.38)E_{r}$ and $t_{so}=0.25E_{r}$
corresponding to the experimental parameters
$V_{\uparrow(\downarrow)}=3.2(2.2)E_{r}$ and $M_{R}=2.0E_{r}$, and the Zeeman
splitting is $m_{z}=0.3E_{r}$. While 2D lattice model may fail in the high-
energy and short wavelength sectors for the absence of higher bands, long-
range hoppings and the discretization of the continuous direction, it could
approximate the low-energy and long wavelength sector of the original system
well. Here the step-length in the $x$ direction is chosen to be
$x_{0}=a_{0}/7$. Figure S13a shows that in this case the periodic boundary
solutions with the real part of the eigenenergies satisfying $\Re E<1.5E_{r}$
already converge well. More quantitatively, the convergence could be
demonstrated by the number of eigenstates in this regime as shown in Figure
S13b. Since the low-energy solutions of the lattice model are good
approximation of the real eigenstates, we could study the NHSE by solving
$H_{lattice}$ under the OBC.

Figure S 13: The convergence of the finite difference method. a Low energy
part of the periodic boundary spectrum of the lattice model in Eq. (8) with
different step-lengths in the $x$ direction. b Number of eigenstates with the
real part of eigenenergy satisfying $\Re E<1.5E_{r}$ varying with the step-
length in the $x$ direction.

The lattice model exhibits NHSE when $\gamma_{\downarrow}\neq 0$, with a
volume-law number of eigenstates being localized on the corners of the open
boundary as shown by Figure S12b2. A noteworthy feature is that besides the
most of the skin modes being localized at the right-bottom corner, there are
also skin modes showing up in the opposite left-top corner. We attribute this
to the (approximate) particle-hole(-like) symmetry (PHS). If
$t_{\uparrow}^{y}=t_{\downarrow}^{y}$, the lattice model has a symmetry
$\hat{S}=(-1)^{j_{x}}\sigma_{x}\hat{P}$, where $\sigma_{x}$ is a Pauli matrix
acting on the spin degree of freedom and $\hat{P}$ is the transpose operation
under the basis $|\vec{j}s\rangle$ (which is the particle-hole conjugation
$\Pi_{\vec{j},s}(c^{\dagger}_{\vec{j},s}+c_{\vec{j},s})$ in the second
quantization representation). The measurement of $\Delta E(k_{x},q_{y})$
reveals the nontrivial spectral winding topology of quasi-1D layers, say
$E_{1}(q_{y})$ with a fixed $k_{x0}$. And then the particle-hole pair related
to one band $E_{1}(q_{y})$ by $\hat{S}$, $-E_{1}(-q_{y})$, should wind in the
opposite direction. Therefore the particle and hole bands have opposite
spectral winding numbers, indicating that the skin modes are distributed in
different (opposite) boundaries sato2020-origin . For the same reason, in the
$k_{x}$ direction this counter-propagating spectral flows also exist in the
low-energy regime. Therefore the system hosts 2D symmetric skin effect with
skin modes being localized in both of the opposite corners as shown in the
plot. In fact, this local PHS relates pairs of skin modes like
$|E\rangle,|-E\rangle$ with opposite eigenenergies $\pm E$, which if not being
topological boundaries states should be localized in opposite boundaries
liu2023-PHS . Figure S12b3 show that the 50 eigenstates with eigenenergies
most closest to the base point $E_{0}=(1.1+0.08\text{I})E_{r}$ are being
localized at the top left corner, and correspondingly the 50 eigenstates with
eigenenergies most closest to $-E_{0}$ are being localized at the bottom right
corner as shown in Figure S12b4. It is a novel phenomenon arising in the non-
Hermitian physics that the local symmetry establishes a nonlocal
correspondence.

Figure S 14: Numerical simulations of the spin polarization of the lowest band
at equilibrium. a Spin polarization of the lowest band of the 2D lattice model
in Eq. (8), with $\gamma_{\downarrow}=0$ and
$t_{\uparrow}=t_{\downarrow}=0.35E_{r}$. b Spin polarization of the lowest
band of the 2D lattice model in Eq. (8), with
$t_{\uparrow}=0.32Er,t_{\downarrow}=0.38E_{r}$. c Spin polarization of the
lowest band of the model of the experimental system at equilibrium with
$V_{\uparrow}=V_{\downarrow}=0.27E_{r}$, calculated by the plane wave method.
d Spin polarization of the lowest band of the model of the experimental system
at equilibrium with $V_{\uparrow(\downarrow)}=0.32E_{r}(0.22E_{r})$,
calculated by the plane wave method.

Our real experimental system deviates from this $\hat{S}$ symmetry, mainly
resulting from the high-energy sectors (higher bands beyond the energy scale
of the s-band in the optical lattice) and that $t_{\uparrow}\neq
t_{\downarrow}$. However, we know that the s-band width being about
$w=2t_{\uparrow}+2t_{\downarrow}$ is much larger than the strength of the
perturbation breaking the PHS at about $\delta
w=2|t_{\uparrow}-t_{\downarrow}|$. Since any pair of eigenstates related by
the PHS have opposite eigenenergies, such pair of skins modes in the low-
energy sector whose eigenenergies are larger than $\delta w$ will not be
largely influenced by the perturbation and their localization properties
should be preserved. Therefore we could still predict this nonlocal feature of
the skin effect for the real system. More concretely, the low-energy
properties of the real system are not qualitatively changed compared with the
ideal model in Eq. (8). We visualise the symmetry breaking by simulating the
spin polarization of the lowest band at equilibrium as shown in Figure S14: In
the ideal symmetric system described by Eq. (8) with
$t_{\uparrow}=t_{\downarrow}$, only the physics within the energy scale of
s-band width of the optical lattice is kept. In this case, the $\hat{S}$
symmetry is perfectly respected, reflected in the features that the Dirac
points are located at $q_{y}/k_{0}=0,\pm 1$ and the reflection symmetry about
the $q_{y}=0$ axis of the spin polarization distribution as shown in Figure
S14a. If $t_{\uparrow}\neq t_{\downarrow}$, the Dirac points are slightly
perturbed from the symmetric lines $q_{y}=0,\pm\pi$ and the reflection
symmetric feature is slightly broken as shown in Figure S14b. If the higher-
bands are included and the optical lattice is spin independent, the Dirac
points are moved away from the symmetric lines and the spin polarization image
is not perfectly symmetric about $q_{y}=0$ as shown by Figure S14c. In the
real experimental system with spin-dependent optical lattices, the band
structure is modified as shown by Figure S14d. We could see that the symmetric
features protected by the PHS are broken to a small extent in the real system.
Also, the simulation shown in Figure S12b1 already shows the case of
$t_{\uparrow}\neq t_{\downarrow}$ in which the nearly-symmetric feature is
observed. The robustness of the nonlocal skin effect with PHS significantly
distinguishes PHS from time-reversal symmetries sato2020-origin ;
critical-2020 ; hofmann2020-reciprocal  and unitary spatial symmetries
highorder-2020 ; mirror-2020 , in the case of which the degenerate pairs
related by the symmetry could be scattered to each other by an arbitrarily
small perturbation breaking the symmetry, resulting in the absence of skin
effect. The nonlocal nature of the fundamental symmetries like the PHS
manifested in non-Hermitian physics, and the symmetry protected feature of the
NHSE are the new directions of our future research. As an aside, the
eigenstates with energy much higher than the s-band energy scale should
recover the Hermiticity and become more extended in the real space, simply
because the kinetic energy part dominates the short wavelength physics.

Figure S 15: Simulated 2D momentum dependent time evolution for spin
polarization texture $P(k_{x},q_{y},t)$.

## Supplementary References

  * (1) B. Song, L. Zhang, C. He, T. F. J. Poon, E. Hajiyev, S. Zhang, X.-J. Liu and G.-B. Jo, Science advances 4, eaao4748 (2018).
  * (2) Z. Ren, D. Liu, E. Zhao, C. He, K. K. Pak, J. Li and G.-B. Jo, Nature Physics 18, 385 (2022).
  * (3) B. Song, C. He, S. Niu, L. Zhang, Z. Ren, X.-J. Liu and G.-B. Jo, Nature Physics 15, 911 (2019).
  * (4) L. Zhang, L. Zhang and X.-J. Liu, Physical Review A 99, 053606 (2019).
  * (5) C.-R. Yi, L. Zhang, L. Zhang, R.-H. Jiao, X.-C. Cheng, Z.-Y. Wang, X.-T. Xu, W. Sun, X.-J. Liu, S. Chen et al., Physical review letters 123, 190603 (2019).
  * (6) N. Okuma, K. Kawabata, K. Shiozaki and M. Sato, Physical review letters 124, 086801 (2020).
  * (7) K. Zhang, Z. Yang and C. Fang, Nature communications 13, 1 (2022).
  * (8) X.-J. L. Zhi-Yuan Wang, Jian-Song Hong, Phys. Rev. B 108, L060204 (2023).
  * (9) Z. Y. Chun-Hui Liu, Kai Zhang and S. Chen, Phys. Rev. Research 2, 043167 (2020).
  * (10) T. Hofmann, T. Helbig, F. Schindler, N. Salgo, M. Brzezińska, M. Greiter, T. Kiessling, D. Wolf, A. Vollhardt, A. Kabaši et al., Physical Review Research 2, 023265 (2020).
  * (11) M. S. Kohei Kawabata and K. Shiozaki, Phys. Rev. B 102, 205118 (2020).
  * (12) T. M. Tsuneya Yoshida and Y. Hatsugai, Phys. Rev. Research 2, 022062(R) (2020).