# Proposal for realizing and probing topological crystalline insulators in
optical lattices

Jing-Xin Liu<EMAIL_ADDRESS>National Laboratory of Solid State
Microstructures and School of Physics, Nanjing University, Nanjing 210093,
China  Collaborative Innovation Center of Advanced Microstructures, Nanjing
210093, China  Jian-Te Wang  National Laboratory of Solid State
Microstructures and School of Physics, Nanjing University, Nanjing 210093,
China  Collaborative Innovation Center of Advanced Microstructures, Nanjing
210093, China  Shi-Liang Zhu<EMAIL_ADDRESS>Key Laboratory of Atomic and
Subatomic Structure and Quantum Control (Ministry of Education), Guangdong
Basic Research Center of Excellence for Structure and Fundamental Interactions
of Matter, School of Physics, South China Normal University, Guangzhou 510006,
China  Guangdong Provincial Key Laboratory of Quantum Engineering and Quantum
Materials, Guangdong-Hong Kong Joint Laboratory of Quantum Matter, Frontier
Research Institute for Physics, South China Normal University, Guangzhou
510006, China  Quantum Science Center of Guangdong-Hong Kong-Macao Greater Bay
Area, Shenzhen, China

###### Abstract

We develop a lattice model which exhibits topological transitions from $Z_{2}$
topological insulators to mirror symmetry-protected topological crystalline
insulators by introducing additional spin-orbit coupling terms. The
topological phase is characterized by the mirror winding number, defined
within the mirror symmetry invariant subspace, which ensures the protection of
gapless edge states and zero-energy corner states under specific boundary
conditions. Additionally, we propose a feasible scheme using ultracold atoms
confined in a stacked hexagonal optical lattice with Raman fields to realize
the two-dimensional topological crystalline insulators. Detection of the
mirror winding number in these systems can be achieved by implementing a
simple quench sequence and observing the evolution of the time-of-flight
patterns.

††preprint: APS/123-QED

##  I Introduction

Topological insulators (TIs) exhibit the bulk-boundary correspondence behavior
and generate symmetry-protected edge statesHasan and Kane (2010); Qi and Zhang
(2011).The bulk is insulating but supports gapless edge states, and these
local boundary states are robust against perturbations as long as the
corresponding symmetry is not broken. Beyond the bulk-boundary correspondence,
Benalcazar et. al. proposed the concept of electric multi-pole
insulatorsBenalcazar  _et al._ (2017a, b), which exhibit bulk-corner
correspondence. In 2D, these insulators do not have edge states but instead
have isolated corner states corresponding to quantized higher electric multi-
pole moments. This extends the notion of correspondences between
$d$-dimensional bulk and $(d-n)$-dimensional boundaries to the concept of
higher-order TIs. Consequently, there is a new possibility that TIs, which
were originally considered to be in the trivial phase, may actually hide
richer higher-order phasesLiu and Wakabayashi (2017); Langbehn  _et al._
(2017); Li  _et al._ (2017); Ezawa (2018a); Geier  _et al._ (2018); Sheng  _et
al._ (2019); Trifunovic and Brouwer (2019); Călugăru  _et al._ (2019); Park
_et al._ (2019); Pan  _et al._ (2019); Zhang  _et al._ (2020); Lee  _et al._
(2020); Trifunovic and Brouwer (2021).

The appearance of corner modes is typically protected by crystalline
symmetries or a combination of other symmetries. Therefore, we can use a broad
class known as topological crystalline insulators (TCIs) to define systems
with corner, edge, and hinge modes Fu (2011); Hsieh  _et al._ (2012).
Describing their topological properties involves creating topological
classification tables that consider crystalline symmetry as well as all
symmetries within the Altland-Zirnbauer tenfold classification scheme Chiu
_et al._ (2016, 2013); Langbehn  _et al._ (2017); Cornfeld and Chapman (2019).
Furthermore, different mechanisms exist to distinguish and define the
topological invariants of TCIs. One class can be straightforwardly classified
by considering the bulk and boundary gap. However, for multi-pole higher-order
topological insulators, the phase transition from a trivial insulator is
characterized by a Wannier band gapGhorashi  _et al._ (2020); Kheirkhah  _et
al._ (2020); Li  _et al._ (2020); Wu  _et al._ (2020); Ezawa (2020); Asaga and
Fukui (2020); Claes and Hughes (2020); Khalaf  _et al._ (2021); Yang  _et al._
(2020); Ezawa (2018b). For more detailed classification of higher-order
topological insulatorJia  _et al._ (2023); Yang  _et al._ (2020), multi-pole
moments Benalcazar  _et al._ (2017a); Schindler  _et al._ (2018a); Zhang  _et
al._ (2013); Benalcazar  _et al._ (2017b); Song  _et al._ (2017); Peng  _et
al._ (2017) and bulk polarization Liu and Wakabayashi (2017) are used to
identify Type-I higher-order TIs, while nested Wilson loops and Wannier band
polarizations Khalaf  _et al._ (2021) are applied to Type-II higher-order TIs.
Therefore, the existence of topologically protected corner modes arises from
various physical mechanisms. Experimental realizations of TCIs states have
been widely reported in materials like bismuth Schindler  _et al._ (2018b) and
various synthetic systems such as phononic metamaterialsSerra-Garcia  _et al._
(2018); Ni  _et al._ (2020); Du  _et al._ (2022); Chen  _et al._ (2022),
photonic metamaterials Mittal  _et al._ (2019), and electronic circuits Imhof
_et al._ (2018); Peterson  _et al._ (2018).

Remarkably, ultracold atoms confined in optical lattices provide an excellent
platform for emulating a wide range of systems in condensed matter physics
Zhang  _et al._ (2018); Bloch  _et al._ (2008); Lewenstein  _et al._ (2007);
Bloch  _et al._ (2012); Schäfer  _et al._ (2020). Synthetic gauge fields and
spin-orbit coupling (SOC) can be realized using various techniques, including
trap rotation Fetter (2009), microrotation Sørensen  _et al._ (2005); Lim  _et
al._ (2008, 2010); Mei  _et al._ (2012), and Raman laser-induced transitions
Dalibard  _et al._ (2011); Lin  _et al._ (2009a, 2011); Zhu  _et al._ (2011);
Lin  _et al._ (2009b); Aidelsburger  _et al._ (2011); Osterloh  _et al._
(2005); Li  _et al._ (2022); Jaksch and Zoller (2003); Gerbier and Dalibard
(2010); Goldman  _et al._ (2009); Huang  _et al._ (2016) By combining laser-
induced tunneling and superlattice techniques, strong Abelian Jaksch and
Zoller (2003) and non-Abelian Osterloh  _et al._ (2005) gauge fields can be
achieved, enabling the simulation of topological insulators and various
topological matters Jotzu  _et al._ (2014); Gerbier and Dalibard (2010);
Goldman  _et al._ (2009); Zhu  _et al._ (2006); Goldman, N. and Gaspard, P.
(2007); Shao  _et al._ (2008); Wu (2008); Liang  _et al._ (2023); Stanescu
_et al._ (2009); Zhu  _et al._ (2006); Liu  _et al._ (2007); Beeler  _et al._
(2013); Alba  _et al._ (2011); Huang  _et al._ (2016); Béri and Cooper (2011);
Goldman  _et al._ (2010). Especially, these technologies have been adopted for
realizing atomic spin Hall effects Beeler  _et al._ (2013); Zhu  _et al._
(2006); Liu  _et al._ (2007), quantum anomalous Hall effects Jotzu  _et al._
(2014); Shao  _et al._ (2008); Wu (2008) and 3D spin-orbit couplingWu  _et
al._ (2016); Wang  _et al._ (2018). Based on this, there are proposals to
realize second-order topological superconductors in cold atom systemsZeng  _et
al._ (2019); Wu  _et al._ (2023); Huang and Liu (2020).

In this article, we present an experimentally feasible scheme to achieve
topological crystalline insulators with ultracold atoms in stacked hexagonal
optical lattices possessing mirror symmetry. By incorporating additional
gradient potential and Raman fields, we can establish a model in the quantum
spin Hall phase with $s_{x}$ spin symmetry. Adjustment of nearest-neighbor
spin-flip hopping via Raman processes allows for the breaking of time-reversal
$\mathcal{T}$ and $s_{x}$ symmetry while preserving mirror symmetry
$\mathcal{M}_{x}$. In the TCIs phase, two $\mathcal{M}_{x}$ symmetry-protected
corner states emerge, characterized by a $\mathbb{Z}_{2}$ class mirror winding
number. This mirror winding number can be determined from the dynamics of
time-of-flight (ToF) images Alba  _et al._ (2011); Deng  _et al._ (2014),
utilizing a straightforward quench sequence Hauke  _et al._ (2014); Fläschner
_et al._ (2016); Yi  _et al._ (2023).

The paper is structured as follows: In Section II, we introduce a model that
demonstrates non-trivial topology, edge states in a specific direction, and
zero-energy corner states under certain boundary conditions. We explore the
origin of corner states, the effective edge Hamiltonian, and the bulk
topological invariant known as the mirror winding number. Section III
discusses the experimental setup involving optical lattices and Raman fields
to achieve the mirror symmetry-protected TCIs phase. In Section IV, we outline
a method to extract the mirror winding number from time-of-flight (ToF) images
by manipulating the quench phase. Finally, Section V concludes the paper with
a summary of key findings and conclusions.

##  II Model

We consider two-component fermionic atoms loaded in a tight-binding honeycomb
lattice with spin-flip hopping terms, which can be viewed as a four-band
model. The corresponding coupling terms are shown in Fig.1(a), and the tight-
binding Hamiltonian can be expressed in three parts as follows:

$\displaystyle H_{0}$
$\displaystyle=-\sum_{\mathbf{r}\in\mathbf{R}_{A}}\sum_{i}t_{i}c^{\dagger}_{\mathbf{r}+\mathbf{e}_{i}}c_{\mathbf{r}}+\mathrm{h.c.},$
(1) $\displaystyle H_{1}$
$\displaystyle=\sum_{\mathbf{r}\in\mathbf{R}_{A}}\left(c^{\dagger}_{\mathbf{r}+\mathbf{e}_{0}-\mathbf{e}_{1}}c_{\mathbf{r}}+c^{\dagger}_{\mathbf{r}+\mathbf{e}_{2}-\mathbf{e}_{0}}c_{\mathbf{r}}\right)(J_{1}s^{-}+J_{2}s^{+})-\sum_{\mathbf{r}\in\mathbf{R}_{B}}\left(c^{\dagger}_{\mathbf{r}+\mathbf{e}_{0}-\mathbf{e}_{1}}c_{\mathbf{r}}+c^{\dagger}_{\mathbf{r}+\mathbf{e}_{2}-\mathbf{e}_{0}}c_{\mathbf{r}}\right)(J_{1}s^{-}+J_{2}s^{+})+\mathrm{h.c.},$
(2) $\displaystyle H_{2}$
$\displaystyle=\sum_{\mathbf{r}\in\mathbf{R}_{A}}\left(g_{1}c^{\dagger}_{\mathbf{r}+\mathbf{e}_{1}}c_{\mathbf{r}}+g_{2}^{\ast}c^{\dagger}_{\mathbf{r}+\mathbf{e}_{2}}c_{\mathbf{r}}\right)s^{-}+\sum_{\mathbf{r}\in\mathbf{R}_{A}}\left(g_{2}c^{\dagger}_{\mathbf{r}+\mathbf{e}_{1}}c_{\mathbf{r}}+g_{1}^{\ast}c^{\dagger}_{\mathbf{r}+\mathbf{e}_{2}}c_{\mathbf{r}}\right)s^{+}+\mathrm{h.c.}$
(3)

where $c^{\dagger}_{\mathbf{r}}$ is the creation operator for a fermion
located at site $\mathbf{r}$, and $s^{+}=(s^{x}+is^{y})/2$ represents the spin
raising operator. The unit direction $\mathbf{e}_{j}\ (j=0,1,2)$ is defined in
Fig. 1a . To create a model with non-trivial topology, we set
$J_{1}=J_{2}=J_{R}+iJ_{I}$ and $g_{1}=-g_{2}=ig$, where the parameters
$J_{R,I}$ and $g$ describe the strengths of next-nearest neighbor and nearest
neighbor spin-flip hopping, respectively.

By performing a Fourier transform, the corresponding Bloch Hamiltonian is
given by

$H(\mathbf{k})=\sum_{\mu,\nu=0}^{3}\bm{d}_{\mu\nu}(\mathbf{k})\cdot\left(\bm{s}_{\mu}\otimes\bm{\sigma}_{\nu}\right),$
(4)

where

$\begin{split}\bm{d}_{01}&=-\sum_{j}t_{j}\cos{(\mathbf{k}\cdot\mathbf{e}_{j})},~{}~{}~{}\bm{d}_{02}=-\sum_{j}t_{j}\sin{(\mathbf{k}\cdot\mathbf{e}_{j})},\cr\bm{d}_{13}&=4\cos{\frac{3k_{x}}{2}}\left(J_{R}\cos{\frac{\sqrt{3}k_{y}}{2}}+J_{I}\sin{\frac{\sqrt{3}k_{y}}{2}}\right),\cr\bm{d}_{21}&=2g\cos{\frac{k_{x}}{2}}\cos{\frac{\sqrt{3}k_{y}}{2}},\cr\bm{d}_{22}&=-2g\sin{\frac{k_{x}}{2}}\cos{\frac{\sqrt{3}k_{y}}{2}}.\end{split}$
(5)

Figure 1:  (a) Hopping structure of the model. (b) Energy spectrum for a
lattice with rhombus boundary.(c) Phase diagram: spin Hall phase $g=0$ (black
solid line), topological crystalline insulator phase boundary (black dashed
line). (d) Density distribution of zero-energy states (DoS). (e) Energy
spectrum and DoS of a lattice with a circular boundary.

###  II.1 Non-trivial bulk band topology with a mirror winding number

In the simplest case where $g=0$ and $J_{R}=0$, this model is in the quantum
spin Hall phase with conserved $s_{1}$. By introducing the nearest-neighbor
spin-flip terms, the phase diagram can be calculated from the bulk gap with
eigenenergies $E=\pm\sqrt{A\pm 2\sqrt{B}}$ where
$A=\bm{d}_{01}^{2}+\bm{d}_{02}^{2}+\bm{d}_{13}^{2}+\bm{d}_{21}^{2}+\bm{d}_{22}^{2}$
and
$B=(\bm{d}_{01}^{2}+\bm{d}_{13}^{2})\bm{d}_{21}^{2}+(\bm{d}_{02}^{2}+\bm{d}_{13}^{2})\bm{d}_{22}^{2}+2\bm{d}_{01}\bm{d}_{02}\bm{d}_{21}\bm{d}_{22}$.
The gap results are shown in Fig.1(c). With the increase of $|g|$, in the TCIs
region, we observe two isolated zero-energy corner states for diamond-shaped
boundary in Fig.1(d). To understand this behavior, we first list all
symmetries in this model. In this case, the model retains time-reversal
symmetry $\mathcal{T}=s_{3}\hat{K}$, particle-hole symmetry
$\mathcal{C}=s_{1}\sigma_{3}\hat{K}$ and some crystalline symmetry, i.e.,
mirror symmetry
$\mathcal{M}_{\mu}H(k_{\mu})\mathcal{M}^{-1}_{\mu}=H(-k_{\mu})$, where
$M_{x}=s_{2}\sigma_{1}$ and $M_{y}=s_{2}$.

To demonstrate the unique topology arising from crystalline symmetry, we
introduce a non-zero $J_{R}$ that breaks parity symmetry $\mathcal{P}$, time-
reversal symmetry $\mathcal{T}$, particle-hole symmetry $\mathcal{C}$, and
mirror symmetry $\mathcal{M}_{y}$. Only the combined symmetries
$\mathcal{PT}=s_{3}\sigma_{1}\hat{K}$, $\mathcal{CP}=is_{1}\sigma_{2}\hat{K}$,
chiral symmetry $\mathcal{S}=s_{2}\sigma_{3}$, and mirror symmetry
$\mathcal{M}_{x}=s_{2}\sigma_{1}$ are retained. According to the
classification table of the $\mathcal{P}+\mathrm{AZ}$ classBzdušek and Sigrist
(2017), this model falls into the $\mathrm{CI}$ class and possesses a
$\mathbb{Z}_{2}$ topological invariant.

However, the presence of mirror symmetry $\mathcal{M}_{x}$ leads to symmetry-
invariant subspaces, and the topological invariant can be described by the
mirror-winding number (also known as the Zak phase). There exist two mirror-
invariant lines, namely the $k_{x}=0$ and $k_{x}=2\pi$ planes. Along each
line, we can define the subspace as a decomposition of $H(k)$, where the
Hamiltonian can be represented as $H(0,k_{y})=H^{+}(0,k_{y})\oplus
H^{-}(0,k_{y})$. Here, $H^{\pm}(0,k_{y})$ are spanned by the eigenstates
${\ket{\psi^{\pm}_{n}}}$ such that
$\mathcal{M}_{x}\ket{\psi^{\pm}_{n}}=\pm\ket{\psi^{\pm}_{n}}$:

$\begin{split}H^{\pm}(0,k_{y})&=\pm\bm{d}_{21}\sigma_{0}+\bm{d}_{01}\sigma_{1}\mp\bm{d}_{13}\sigma_{2}\cr\end{split}$
(6)

and this gives an effective 1D model. As $J_{R}$ increases, the bulk gap at
$k_{x}=0$ will close and reopen. For each sector, this effective model is
gapped when $|J_{R}|<\sqrt{3}|J_{I}|$.

The presence of corner modes can be interpreted as resulting from a non-
trivial winding in the mirror-invariant subspace, which can be quantified by
the winding numbersRen  _et al._ (2020):

$\nu_{\pm}=\frac{1}{2\pi}\int_{0}^{4\pi/\sqrt{3}}\mathrm{d}k_{y}~{}\mathrm{arg}\left[\bm{d}_{01}(k_{y})\mp
i\bm{d}_{13}(k_{y})\right].$ (7)

Here it holds that $\nu_{+}+\nu_{-}=0$ which relates to $\mathcal{PT}$
symmetry. For the case of half occupation, the mirror invariant can be defined
as $\nu=(\nu_{+}-\nu_{-})/2$, where $\nu_{+}$ and $\nu_{-}$ are the winding
numbers for two sectors. The wave functions of the two valence bands can then
be expressed as

$\ket{\psi_{1}}=\left[1,-e^{-i\theta},-ie^{-i\theta},i\right]^{T}/2,$ (8)

and the wave function in the other sector is related as
$\ket{\psi_{2}}=s_{1}\ket{\psi_{1}}$.

We calculate the energy spectrum using two different open boundary conditions
(OBC) and observed that the zig-zag case is always gapless, whereas the
armchair case is gapped by a non-zero $g$. The system in the armchair boundary
remains gapped in the topological crystalline insulator phase. However, for
bulk semi-metal phases, the results are more complex. The OBC spectrum and
edge gaps are presented in Appendix A.

###  II.2 Edge theory and corner modes

Since edge states only exist along one direction, topologically protected
corner states emerge at the intersection of two boundaries. To enhance
clarity, we analyze our model using the edge theory Lee  _et al._ (2007);
Sitte  _et al._ (2012); Yan  _et al._ (2018); Ezawa (2018b); Tan  _et al._
(2022), which is analogous to the Jackiw-Rebbi solitons Jackiw and Rebbi
(1976), to provide a more intuitive understanding .

To illustrate the origin of corner states, we perform a Taylor expansion at
the $\Gamma$ point for different boundaries. Considering that the model we are
discussing involves a hexagonal lattice, we focus on the four boundaries of a
rhombus, as depicted in Fig.1(d). As an example, by replacing
$\sin{k_{\perp}}\to-i\partial_{r_{\perp}}$, $\cos{k_{\perp}}\to
1-\frac{1}{2}\partial^{2}_{r_{\perp}}$, and $\sin{k_{\parallel}}\to
k_{\parallel}$ while omitting insignificant $k_{\parallel}^{2}$ terms, the
Hamiltonian at boundary I, shown in Fig.1(b) with
$r^{(I)}_{\parallel}=(\sqrt{3}/2,1/2)$ and
$r^{(I)}_{\perp}=(-1/2,\sqrt{3}/2)$, is given by

$\displaystyle H(-i\partial_{r_{\perp}},k_{\parallel})$
$\displaystyle=H_{0}(-i\partial_{r_{\perp}},k_{\parallel})+H_{1}(-i\partial_{r_{\perp}},k_{\parallel}),$
(9) $\displaystyle H_{0}(-i\partial_{r_{\perp}},k_{\parallel})$
$\displaystyle=\left(A_{0}+A_{2}\partial^{2}_{r_{\perp}}\right)\sigma_{x}+A_{1}\partial_{r_{\perp}}s_{1}\sigma_{3},$
(10)

where $A_{0}=-(t_{0}+2t_{1})$, $A_{1}=i3J_{I}/2$, and
$A_{2}=-(t_{0}+5t_{1})/8$. Here, we set $t_{2}=t_{1}$. $H_{1}$ collects the
remaining terms, and further details are provided in Appendix B.

By solving the eigenvalue equation with zero energy,
$H_{0}\psi_{\alpha}(r_{\perp})=0$, under the boundary condition
$\psi_{\alpha}(0)=\psi_{\alpha}(+\infty)=0$, the zero-mode solution takes the
form
$\psi_{\alpha}(r_{\perp})=e^{ik_{\parallel}r_{\parallel}}f(r_{\perp})\chi_{\alpha}$,
which satisfies the equation:

$\begin{split}A_{2}\partial^{2}_{r_{\perp}}f(r_{\perp})+\alpha
A_{1}\partial_{r_{\perp}}f(r_{\perp})+A_{0}f(r_{\perp})=0.\end{split}$ (11)

Here we have introduced
$s_{x}\sigma_{z}\chi_{\alpha}=\alpha\sigma_{y}\chi_{\alpha}$. The equation can
be solved by the solutions
$f(r_{\perp})=\mathcal{N}_{r_{\perp}}\sin{(\kappa_{1}r_{\perp})}e^{-\kappa_{2}r_{\perp}}$
with $|\kappa_{1}|=\left|\sqrt{4A_{0}A_{2}-A_{1}^{2}}/2A_{2}\right|$,
$\kappa_{2}=A_{1}/2A_{2}$ and normalization factor
$\mathcal{N}_{r_{\perp}}=2\sqrt{\kappa_{2}\left(1+\kappa_{2}^{2}/\kappa_{1}^{2}\right)}$.
The condition $\alpha=+i$ corresponds to $\kappa_{2}>0$ and $f(+\infty)=0$.
Thus, the solution for the $\alpha=+i$ case represents the zero edge mode at
boundary $\mathrm{I}$. The equations for the spinor part are given by
$s_{x}\sigma_{z}\chi_{\alpha}=+i\sigma_{x}\chi_{\alpha}$, which yields

$\displaystyle\chi^{(\mathrm{I})}_{1}=\ket{s_{1}=+1}\ket{\sigma_{1}=-1},$ (12)
$\displaystyle\chi^{(\mathrm{I})}_{2}=\ket{s_{1}=-1}\ket{\sigma_{1}=+1}.$ (13)

Similarly, at the opposite boundary $\mathrm{III}$ with $\alpha=-i$, the
spinor components are

$\displaystyle\chi^{(\mathrm{III})}_{1}=\ket{s_{1}=+1}\ket{\sigma_{1}=+1},$
(14)
$\displaystyle\chi^{(\mathrm{III})}_{2}=\ket{s_{1}=-1}\ket{\sigma_{1}=-1}.$
(15)

Figure 2:  (a) Experimental setup of a stacked optical lattice with blue and
red laser fields. Additional laser fields propagate along the
$\hat{z}$-direction for the Raman process. (b) Example of $~{}^{40}\mathrm{K}$
atom levels chosen to satisfy $\lambda_{R}/\lambda_{B}\approx 2$. A gradient
potential generates a linear energy offset in the optical lattice along one
specific direction. (c) Distribution of optical potential and Raman fields for
blue and red lasers.

The boundary Hamiltonian $H_{\mathrm{I}}$ with the zero boundary mode basis is
given by

$H_{\mathrm{I}}=v_{z}k_{\parallel}\tau_{3}+m_{x}\tau_{1}.$ (16)

For the opposite boundary, we have
$H_{\mathrm{III}}=-H_{\mathrm{I}}=-v_{z}k_{\parallel}\tau_{3}-m_{x}\tau_{1}$.
Similarly, we calculate the boundary Hamiltonians at boundaries $\mathrm{II}$
and $\mathrm{IV}$: $H_{\mathrm{II}}=-v_{z}k_{\parallel}\tau_{3}+m_{x}\tau_{1}$
and $H_{\mathrm{IV}}=v_{z}k_{\parallel}\tau_{3}-m_{x}\tau_{1}$.

The presence of $g$ introduces an effective mass term at each boundary. The
difference in mass terms for the boundary massive Dirac fermion results in the
Jackiw-Rebbi model Jackiw and Rebbi (1976), where the domain wall leads to
corner states. This leads to corner states located at the intersection of
boundaries $\mathrm{I},\mathrm{IV}$ and $\mathrm{II},\mathrm{III}$ as shown in
Fig.1(b,d).

Besides the above example of the rhombus boundary, we also depict the curved
states with circular boundaries, and the results are shown in Fig.1(e). Due to
the effective edge Hamiltonian having no mass term in specific direction, edge
states preferentially form along the $\hat{y}$ direction. This explains the
observed curved states demonstrated in Fig.1(e).

##  III Physical Realization

In this section, we present the details of deriving effective Hamiltonians
based on realistic cold atom candidates and analyze the corresponding tight-
binding model.

The formation of a hexagonal optical lattice Shao  _et al._ (2008); Zhu  _et
al._ (2007); Duan  _et al._ (2003); Soltan-Panahi  _et al._ (2011); Lee  _et
al._ (2009) can be achieved through the intersection of three laser beams at a
$120^{\circ}$ angle. We generate the hexagonal lattice and spin-orbit coupling
(SOC) by stacking optical lattices with wavelengths $\lambda_{R}$ and
$\lambda_{B}$, along with additional laser fields propagating along the $z$
direction as shown in Fig.2(a). To simultaneously match the spatial periods of
the two optical lattices, the wavelengths of the red and blue lasers should
satisfy $\lambda_{R}/\lambda_{B}=2$. The total electric field of optical
lattice is denoted as

$\begin{split}E_{OL}(\mathbf{r})&=\sum_{n=1}^{2}E_{n}e^{-i\omega
t/n}\sum_{j=1}^{3}\cos{\left(n\mathbf{k}_{j}\cdot\mathbf{r}+\frac{n\pi}{2}\right)}\hat{\mathbf{z}},\cr\end{split}$
(17)

where the wavevector $\mathbf{k}_{1}=(1,0)k_{0}$,
$\mathbf{k}_{2}=(-\frac{1}{2},\frac{\sqrt{3}}{2})k_{0}$,
$\mathbf{k}_{3}=(-\frac{1}{2},-\frac{\sqrt{3}}{2})k_{0}$ with the unit
$k_{0}=2\pi/\lambda_{B}$.

The optical potential for atoms in the ground state can generally be expressed
as
$V(\mathbf{r})=u_{s}|\mathbf{E}|^{2}+iu_{v}\left(\mathbf{E}^{\ast}\times\mathbf{E}\right)\cdot\mathbf{S}$,
where the first term represents the scalar potential, and the second term
denotes a vector light shift. The choice of laser wavelength can be addressed
by selecting the $D_{1}$ line and $6^{2}S_{1/2}\to 7^{2}P_{3/2}$ transition of
${}^{133}\mathrm{Cs}$ atoms with a nuclear spin $I=7/2$, resulting in
$\lambda_{D_{1}}^{(b)}/\lambda_{1}^{(b)}=1.963$, where $\lambda_{1}^{(b)}$
represents the wavelength of the $6^{2}S_{1/2}\to 7^{2}P_{3/2}$ transition
Steck (2010). For the fermionic case, we can choose the $D_{1}$ line and
$4^{2}S_{1/2}\to 5^{2}P_{3/2}$ transition of ${}^{40}\mathrm{K}$ atoms with a
nuclear spin $I=4$, giving $\lambda_{D_{1}}/\lambda_{1}^{(f)}=1.904$, where
$\lambda_{1}^{(f)}$ represents the wavelength of the $4^{2}S_{1/2}\to
5^{2}P_{3/2}$ transition Tiecke (2019).

Before discussing TCIs phases, we should first realize four-band topological
insulators (TIs). The realization of a $\mathcal{T}$-symmetric model has been
achieved in an optical lattice through the incorporating of a gradient
magnetic field and Raman laser fields Kennedy  _et al._ (2013); Aidelsburger
_et al._ (2013); Goldman  _et al._ (2010); Zhang  _et al._ (2016). An approach
we can consider involves introducing an additional gradient potential. This
potential can facilitate the independent coupling of spins at different sites.
By applying an artificial electric field via accelerated motion of the optical
lattice, we can derive the 2D effective Hamiltonian as follows:

$H_{\mathrm{eff}}(\mathbf{r})=\frac{\mathbf{p}^{2}}{2m}+V_{OL}(\mathbf{r})+\mathcal{F}\cdot\mathbf{r}+\frac{\mu_{B}g_{F}}{\hbar}\mathbf{B}\cdot\mathbf{F}.$
(18)

By adding laser fields $\mathbf{E}_{R1}$ and $\mathbf{E}_{R2}$ propagating
along the $\hat{z}$ direction, the electric fields in the $z=0$ plane are:

$\displaystyle\mathbf{E}_{R1}$
$\displaystyle=E_{R1}\left(e^{-i(\omega_{R1}+\delta_{1})t}+e^{-i(\omega_{R1}+\delta_{2})t}\right)(\hat{\mathbf{x}}+i\hat{\mathbf{y}}),$
(19) $\displaystyle\mathbf{E}_{R2}$
$\displaystyle=E_{R2}\left(e^{-i(\omega_{R2}+\delta_{1})t}+e^{-i(\omega_{R2}+\delta_{2})t}\right)(\hat{\mathbf{x}}+i\hat{\mathbf{y}}),$
(20)

where $\delta_{1}$ and $\delta_{2}$ compensate for the detuning caused by the
gradient potential during different spin-flip processes. These laser fields
can drive the $\sigma_{+}$ and $\pi$ transitions from the ground states
$\ket{g_{\uparrow}}=\ket{9/2,+9/2}$ and $\ket{g_{\downarrow}}=\ket{9/2,+7/2}$
to all possible excited states that satisfy the selection rules. For the
$D_{1}$ and $D_{2}$ lines in ${}^{40}\mathrm{K}$ atoms, the state
$\ket{9/2,+7/2}$ couples to the excited states $\ket{F,+9/2}$ ($\sigma_{+}$).
The red (blue) Raman fields, as depicted in Fig.2(b), correspond to the
strengths $E_{1}$ ($E_{2}$) and are derived from the $D_{1}$ line
($4^{2}S_{1/2}\to 5^{2}P_{1/2}$) and $D_{2}$ line ($4^{2}S_{1/2}\to
5^{2}P_{3/2}$) transitions.

The strengths of the blue and red lasers’ optical potentials,
$V^{(1,2)}(\mathbf{r})$, and the corresponding Raman fields,
$M_{R1,R2}(\mathbf{r})$, are illustrated in Fig.2(c). Details explanations of
the optical lattice potential, the construction of Raman processes, and the
relevant parameters can be found in Appendix C. The overall optical potential
is a result of contributions from both the blue and red lasers. Fig.2(c)
depicts the potential distributions of the two lasers. By configuring the
system such that $E_{1}\gg E_{2}$, the blue laser’s contribution to the
potential becomes negligible and we don’t need to consider the third
sublattice degree. This setup ensures that the energy minimum points are
located at the positions marked by black solid points, with a lattice constant
of $2\lambda_{B}/3$.

To obtain the effective tight-binding model, we need to calculate the overlap
integral of two Wannier-Stark states. First, let’s consider a rough analysis:
suppose there is no gradient potential. In such a honeycomb lattice, the
Wannier function maintains $C_{3}$ symmetry (since the energy minima regions
exhibit $C_{3}$ symmetry). In this case, the center of the Wannier function
will coincide with the energy minimum point of the potential, as shown by the
black point in Fig.2(c). The corresponding strengths of next-nearest-neighbor
and nearest-neighbor spin-flip couplings can be calculated by overlapping two
Wannier functions, taking into account the Raman field background. For the
contributions of the red lasers, the nearest-neighbor spin-flip terms can be
eliminated due to the $C_{6}$ anti-symmetry and $C_{3}$ symmetry of the
Wannier function in Fig.2(c). By controlling the strength of $M_{R1}$, the
next-nearest-neighbor spin-flip term can become sufficient large. The nearest-
neighbor spin-flip hopping is controlled by the blue Raman field $M_{R2}$. The
overlap of the next-nearest-neighbor Wannier function is small enough that we
can disregard the next-nearest-neighbor terms contributed by the blue Raman
field. The phase and strength of the Raman fields $M_{R1}$ and $M_{R2}$ can
also be controlled independently by adjusting $E_{R1}$ and $E_{R2}$. By tuning
the Rabi frequencies of the laser fields to satisfy $J_{1}=J_{2}$ and
$g_{1}=-g_{2}=-ig$, we can realize the effective lattice model described in
Section II.

However, the presence of a gradient potential in Eq. 18 breaks the original
$C_{6}$ symmetry of the potential while retaining mirror symmetry
$\mathcal{M}_{x}$. Suppose the gradient potential along the $\hat{y}$
direction is not large enough, so the centers of the Wannier-Stark functions
are only slightly displaced. In this situation, the overlap of two next-
nearest-neighbor Wannier-Stark functions also results in opposite values for
different sublattices. Additionally, there is an extra contribution from the
nearest-neighbor spin-flip term caused by the blue Raman field. This
introduces staggered coupling, which is also influenced by the $C_{6}$ anti-
symmetry of the blue Raman field (as detailed in Appendix D). The effect of
these terms exhibits a similar form to the contribution from the red Raman
field and will not break mirror symmetry $\mathcal{M}_{x}$.

##  IV Probe Mirror winding number

In the $\mathcal{M}_{x}$-invariant subspace, the system can be viewed as two
isolated 1D models, and its topology can be simply described by a single
parameter $\theta(k_{y})$ mentioned in Sec II. In the cold atom platform, this
parameter of Bloch states can be extracted from interference patterns of ToF
imagesHauke  _et al._ (2014); Fläschner  _et al._ (2016); Zhang  _et al._
(2018); Hauke  _et al._ (2014); Fläschner  _et al._ (2016). In this section,
we develop a method to probe the mirror winding number.

For this four-band model, the overall density distribution of the ToF image in
momentum space can be calculated using the following expression:

$\begin{split}n(\mathbf{k})&=f(\mathbf{k})\sum_{m,n}e^{-i\mathbf{k}\cdot(\mathbf{r}_{m}-\mathbf{r}_{n})}e^{i(\mu_{m,\sigma}-\mu_{n,\sigma})t}\cr&\times\bra{G}C_{m}^{\dagger}C_{n}\ket{G},\cr\end{split}$
(21)

where
$C_{n}=[a_{n,\uparrow},b_{n,\uparrow},a_{n,\downarrow},b_{n,\downarrow}]^{T}$
and $a_{n,\sigma}$($b_{n,\sigma}$) represents annihilation operators for site
$\mathbf{r}_{n}$ in sublattice A(B) with spin $\sigma$. The factor
$f(\mathbf{k})$ denotes a broad envelope function determined by the momentum
distribution of the Wannier function, and $\mu_{n,\sigma}$ indicates the
global trap potential strength at site $n$. In the absence of a global
trapping difference, we have $\mu_{n,\sigma}=\mu_{s}$. The state $\ket{G}$
represents the many-particle states with Fermi energy $E_{f}$ in an ideal
periodic system without trapping differences. The ToF image can be expressed
as
$n(\mathbf{k})=f(\mathbf{k})\bra{G}s_{0}\sigma_{0}+s_{0}\sigma_{1}+s_{1}\sigma_{0}+s_{1}\sigma_{1}\ket{G}$,
which contains all contributions from interference. In the case of periodic
boundary condition (PBC), we sum the contribution of the two valence bands
using Eq.(8). The contributions from the $s_{1}\sigma_{0}$ and
$s_{1}\sigma_{1}$ terms are zero, yielding $n(k_{y})=2(1-\cos{\theta(k_{y})})$
theoretically. In this scenario, a quench involving the Pauli matrix
$\sigma_{3}$ generates the evolution operator $\exp{(-i\sigma_{3}\varphi/2)}$,
which can extract $\theta(k_{y})$ independently. This quench sequence leads to
the ToF image into

$n(0,k_{y},\varphi)=2(1-\cos{\varphi}\cos{\theta(k_{y})}).$ (22)

and the parameter $\theta(k_{y})$ can be extracted via adjusting the quench
phases.

However, in a finite-size system, the ground states differ from Bloch waves in
PBC. Consequently, the numerical results, including contributions from corner
states, don’t precisely match Eq. 22. The actual results can be expressed
similarly but with an additional background term due to OBC:
$n_{1}(k_{y},\varphi)=f(k)(1-\cos{\varphi}\cos{\theta})+b(k)$ where $b(k)$ is
background caused by OBC. To account for the broadening caused by Wannier
functions $f(k)$ and the background contributions from corner modes, we
introduce another quench term $\exp{(-i\sigma_{2}\xi/2)}$, where
$n_{2}(k_{y},\xi)=f(k)(1-\cos{(\theta-\xi)})+b(k)$. Consequently, we have
$f(k)=(\mathrm{max}\\{n_{2}(k)\\}-\mathrm{min}\\{n_{2}(k)\\})$ and
$b(k)=\mathrm{min}\\{n_{2}(k)\\}$.

Figure 3:  (a) Bulk gap versus $J_{R}$ and $J_{I}$. The red solid line
indicates the boundary of topological phase transition. The black dashed line
marks the boundary of the gapless region. (b-e) Extracted $\theta_{\pm}$ for
two $\mathcal{M}_{x}$-invariant subspaces. (b) Topological phase with
$J/t=0.4i$; (c) $J/t=0.2+0.4i$; (d) gapless region with $J/t=0.8+0.4i$; (e)
trivial phase with $J/t=1.2+0.4i$ in open boundary condition lattice with
$8\times 8\times 2$ sites and $g/t=0.2$. Subgraphs show modified ToF image
evolution for different momentum $k_{y}$ by varying quench phases.

In Eq. 6, we see that the mirror winding number only relates to the ratio
$|J_{R}/J_{I}|$. However, the additional spin-orbit coupling $g$ changes the
bulk band structure from a gapped phase to gapless one. This implies that in
some regions, we cannot extract the topological invariant because the valance
bands are not well-defined. Therefore, even though the mirror winding numer
remains non-trivial in some regions, we cannot extract it. We give the bulk
gap for a system with constant nearest neighbour coupling $g$ in Fig.3(a). The
red line $|J_{R}/J_{I}|=\sqrt{3}$ marks the boundary between topological and
trivial phases. For regions with finite gap, we can extract information
related to the topology of the valence bands.

From the ToF image, we can extract $\theta(k_{y})$ as shown in Fig.3(b-e). In
a system in the TCIs phase with a non-trivial mirror winding number and a
finite bulk gap, $\cos{\theta(k_{y})}$ falls within the range of $[-1,1]$. In
other gapless regiona, the ToF image contains contributions from both valence
and conduction bands, and the winding number is not well-defined. For a system
in the trivial phase with a finite bulk gap, $\cos{\theta}$ cannot cover the
entire value range and the mirror winding number must be zero.

##  V Conclusion

We have developed a lattice model that exhibits $Z_{2}$ TIs. By introducing
additional SOC terms, we observe a phase transition from $Z_{2}$ TIs to mirror
symmetry-protected TCIs. The topological invariant is defined in the mirror-
symmetry-invariant subspace as the $\mathbb{Z}_{2}$ mirror winding number.
This mirror symmetry protects gapless edge states at specific boundaries, and
zero-energy corner states can be observed in particular open boundary
conditions. We have proposed a feasible scheme using cold atoms with optical
lattices and Raman fields to realize the TCIs, where all parameters in the
Hamiltonian can be independently adjusted. The corresponding mirror winding
number can also be probed by measuring the ToF image. This work provides a
method to control and simulate various topological states of matter in cold
atom systems.

## Acknowledgment

We thank Peng He and Y.Q. Zhu for helpful discussions. This work was supported
by National Key R&D Program of China (Grant No. 2022YFA1405300) and the
National Natural Science Foundation of China (Grant No. 12074180).

##  Appendix A Edge gap

Figure 4:  Edge energy spectrum with zig-zag boundary in the x-direction OBC
(a) and armchair boundary of y-direction OBC (b) with parameters
$g/t=0.2,J_{I}/t=0.4$. (c) The edge gap of armchair boundary exhibits the same
behavior as the bulk gap in the TCIs phase region.

Here, we present the edge band spectrum for zig-zag and armchair open boundary
cases. With the introduction $g$, a gap opens at specific directions. From the
edge gap diagram, we can clearly observe the vanishing of edge mode regions,
which exhibits similar behaviors to the bulk gap diagram.

##  Appendix B Boundary Hamiltonian

To demonstrate the origin of corner modes, we perform a Taylor expansion at
the $\Gamma$ point for different boundaries. Given that the model we are
discussing involves a hexagonal lattice, we can consider the boundaries with
wave vectors $k_{1}=k_{x},~{}k_{2}=k_{x}/2-\sqrt{3}k_{y}/2$ and
$k_{3}=k_{x}/2+\sqrt{3}k_{y}/2$. We replace $\sin{k_{\perp}}$ with
$\to-i\partial_{r_{\perp}}$ and $\cos{k_{\perp}}$ with $1$, while omitting
insignificant $k_{\parallel}^{2}$ terms.

The Hamiltonian at boundary I, with
$r^{(I)}_{\parallel}=\sqrt{3}\hat{x}/2+\hat{y}/2$ and
$r^{(I)}_{\perp}=-\hat{x}/2+\sqrt{3}\hat{y}/2$, can be decomposed into three
parts:

$\begin{split}H_{0}&=-\left[t_{0}+2t_{1}+\frac{1}{8}(t_{0}+5t_{1})\partial^{2}_{r_{\perp}}\right]\sigma_{1}+\left(i\frac{3}{2}J_{I}\partial_{r_{\perp}}\right)s_{1}\sigma_{3}\cr
H_{1}&=i\frac{\sqrt{3}}{4}k_{\parallel}(t_{0}-t_{1})\partial_{r_{\perp}}\sigma_{1}+J_{R}\left(2+\frac{9}{8}\partial^{2}_{r_{\perp}}\right)s_{1}\sigma_{3}\cr\cr&-\left[\frac{\sqrt{3}}{2}J_{I}\left(1-\frac{9}{8}\partial^{2}_{r_{\perp}}\right)k_{\parallel}+i\frac{3\sqrt{3}}{4}J_{R}k_{\parallel}\partial_{r_{\perp}}\right]s_{1}\sigma_{3}\cr&-\frac{1}{16}g\left[\sqrt{3}\left(8+\partial^{2}_{r_{\perp}}\right)k_{\parallel}+8i\partial_{r_{\perp}}\right]s_{2}\sigma_{2}\cr
H_{2}&=\frac{1}{16}(t_{0}-t_{1})\left[\sqrt{3}\left(\partial^{2}_{r_{\perp}}+8\right)k_{\parallel}+8i\partial_{r_{\perp}}\right]\sigma_{2}\cr&+\frac{1}{8}g\left[5\partial^{2}_{r_{\perp}}+2\sqrt{3}i\partial_{r_{\perp}}k_{\parallel}+16\right]s_{2}\sigma_{1}\cr\end{split}$
(23)

where the first part $H_{0}$ has a zero-energy solution. The second and third
parts, $H_{1}$ and $H_{2}$, include all other terms in a complex form. The
zero-energy solutions are given in Sec.II. Here, we express the Hamiltonian in
the subspace consisting of two zero-energy edges states:

$H_{\mathrm{edge}}=\left(\begin{matrix}\bra{\psi_{1}}H\ket{\psi_{1}}&\bra{\psi_{1}}H\ket{\psi_{2}}\\\
\bra{\psi_{2}}H\ket{\psi_{1}}&\bra{\psi_{2}}H\ket{\psi_{2}}\\\
\end{matrix}\right)$ (24)

At boundary I, only the terms involving $\sigma_{2}$ and $s_{2}\sigma_{1}$
have non-zero values. Consequently, the contributions from $H_{0}$ and $H_{1}$
are zero in this subspace. In this scenario, the effective continuum
Hamiltonian at the boundary can be represented as:

$\begin{split}H_{I}&=v_{z}k_{\parallel}\sigma_{z}+m_{x}\sigma_{x}.\end{split}$
(25)

and this is essentially a massive Weyl fermion. The velocity and mass term can
be calculated using these integrations:

$\begin{split}\int_{0}^{+\infty}\mathrm{d}r_{\perp}f^{\ast}(r_{\perp})\partial_{r_{\perp}}f(r_{\perp})&=0\cr\int_{0}^{+\infty}\mathrm{d}r_{\perp}f^{\ast}(r_{\perp})\partial^{2}_{r_{\perp}}f(r_{\perp})&=-\mathcal{N}_{y}^{2}\frac{\kappa_{1}^{2}}{4\kappa_{2}}\cr\end{split}$
(26)

##  Appendix C Optical lattice and Raman field setting

The red Raman fields, related to strength $E_{1}$, are contributed by the
$D_{1}$ line and $D_{2}$ line transitions, where
$E_{+1}^{(1)}(\mathbf{r})=E_{R1}$ and
$E_{0}^{(1)}(\mathbf{r})=E_{1}\sum_{j=1}^{3}\cos{2\mathbf{k}_{j}\cdot\mathbf{r}}$.
The optical potential formed by red lasers in Fig.2(c) for spin $\tau$ is

$\begin{split}V^{(1)}_{\tau}(\mathbf{r})&=\sum_{n=1}^{2}\sum_{F}\left(\frac{|\Omega^{D_{n}}_{\tau
F,+1}|^{2}}{\Delta_{D_{n}}}+\frac{|\Omega^{D_{n}}_{\tau
F,0}|^{2}}{\Delta_{D_{n}}}\right).\cr\end{split}$ (27)

where $\Omega_{\tau F}$ represents the Rabi frequency of all possible
transitions from the state $\ket{4^{2}S_{1/2},9/2,+7/2}$. The Raman field
caused in Fig.2(c) is

$\begin{split}M_{R1}(\mathbf{r})&=\sum_{F}\left(\frac{\Omega^{D_{2}\ast}_{\uparrow
F,+1}\Omega^{D_{2}}_{\uparrow
F,0}}{\Delta_{D_{2}}}+\frac{\Omega^{D_{1}\ast}_{\uparrow
F,+1}\Omega^{D_{1}}_{\uparrow
F,0}}{\Delta_{D_{1}}}\right)\cr&\left(\sum_{j=1}^{2}\delta(\Delta+\omega-\omega_{R1}-\delta_{j}(\mathbf{r}))\right)\ket{g_{\downarrow}}\bra{g_{\uparrow}}+h.c.\cr\end{split}$
(28)

where $\Delta$ is the level splitting caused by the external Zeeman field and
$\delta$ is the delta function. For example, the Rabi frequencies of the
$D_{1}$ line mentioned in main text are given by $|\Omega_{\uparrow
F,\sigma}^{D_{1}}|=|\mu_{\uparrow F}^{D_{1}}||\vec{E}_{\sigma}|/\hbar$, where
$\mu^{D_{1}}_{\uparrow,F}$ is the electric dipole moment:

$\begin{split}\mu^{D_{1}}_{\uparrow,F}&=e\bra{4^{2}S_{1/2},9/2,+9/2}\mathbf{r}\cdot\hat{\epsilon}_{\sigma}\ket{4^{2}P_{1/2},F,+9/2+\sigma}\cr\end{split}$
(29)

The $\hat{\epsilon}_{\sigma}$ represents the polarization direction of the
eletric field. Because the strength of $E_{+1}^{(1)}$ is constant, the Raman
field $M_{R1}$ can be represented as

$M_{R1}(\mathbf{r})=\mathrm{Re}[M_{1}(\mathbf{r})]\sigma_{x}+\mathrm{Im}[M_{1}(\mathbf{r})]\sigma_{y},$
(30)

where $M_{1}(\mathbf{r})$ is proportional to
$\sum_{j=1}^{3}\cos{2\mathbf{k}_{j}\cdot\mathbf{r}}$.

The blue Raman fields shown in Fig.2(b) are related to the strength $E_{2}$
and arise from the $4^{2}S_{1/2}\to 5^{2}P_{1/2}$ and $4^{2}S_{1/2}\to
5^{2}P_{3/2}$ transitions. Here, $E_{+1}^{(2)}(\mathbf{r})=E_{R2}$ and
$E_{0}^{(2)}=E_{2}\sum_{j=1}^{3}\sin{\mathbf{k}_{j}\cdot\mathbf{r}}$. The
optical potential for spin $\tau$ can be calculated using the A.C. Stark
effect and is denoted as

$\begin{split}V^{(2)}_{\tau}(\mathbf{r})&=\sum_{F}\left(\frac{|\Omega^{5^{2}P_{3/2}}_{\tau
F,+1}|^{2}}{\Delta_{5^{2}P_{3/2}}}+\frac{|\Omega^{5^{2}P_{3/2}}_{\tau
F,0}|^{2}}{\Delta_{5^{2}P_{3/2}}}\right)\cr&+\sum_{F}\left(\frac{|\Omega^{5^{2}P_{1/2}}_{\tau
F,+1}|^{2}}{\Delta_{5^{2}P_{1/2}}}+\frac{|\Omega^{5^{2}P_{1/2}}_{\tau
F,0}|^{2}}{\Delta_{5^{2}P_{1/2}}}\right).\cr\end{split}$ (31)

The Raman field caused in Fig.2(b) is

$\begin{split}M_{R2}(\mathbf{r})&=\sum_{F}\left(\frac{\Omega^{5^{2}P_{3/2}\ast}_{\uparrow
F,+1}\Omega^{5^{2}P_{3/2}}_{\uparrow
F,0}}{\Delta_{5^{2}P_{3/2}}}+\frac{\Omega^{5^{2}P_{1/2}\ast}_{\uparrow
F,+1}\Omega^{5^{2}P_{1/2}}_{\uparrow
F,0}}{\Delta_{5^{2}P_{1/2}}}\right)\cr&\left(\sum_{j=1}^{2}\delta(\Delta+\omega-\omega_{R2}-\delta_{j}(\mathbf{r}))\right)\ket{g_{\downarrow}}\bra{g_{\uparrow}}+h.c.,\cr\end{split}$
(32)

The Raman field $M_{R2}$ can be represented as

$M_{R2}(\mathbf{r})=\mathrm{Re}[M_{2}(\mathbf{r})]\sigma_{x}+\mathrm{Im}[M_{2}(\mathbf{r})]\sigma_{y},$
(33)

where $M_{2}(\mathbf{r})$ is proportional to
$\sum_{j=1}^{3}\sin{\mathbf{k}_{j}\cdot\mathbf{r}}$.

##  Appendix D Tight-binding parameters

Figure 5:  Illustration of nearest-neighbor spin-flip terms caused by blue
Raman fields. The imbalanced distribution of Wannier-Stark functions leads to
staggered coupling.

The parameters of the tight-binding model can be calculated from the Wannier-
Stark functions. We define the nearest-neighbor hopping with the Wannier-Stark
state $w(\mathbf{r})$:

$\begin{split}t_{\mu}&=\int\mathrm{d}^{2}\mathbf{r}~{}w_{B\sigma}^{\ast}(\mathbf{r}-\mathbf{R}_{A}-\mathbf{e}_{\mu})\left(\frac{\mathbf{\mathbf{p}}^{2}}{2m}+V_{OL}(\mathbf{r})\right)\cr&\times
w_{A\sigma}(\mathbf{r}-\mathbf{R}_{A}).\cr\end{split}$ (34)

Mirror symmetry enforces $t_{1}=t_{2}\neq t_{0}$. The blue Raman fields
$M_{1}(\mathbf{r})$ generates nearest-neighbor spin-flip terms

$\begin{split}g_{A\mu}&=\int\mathrm{d}^{2}\mathbf{r}~{}w_{B\downarrow}^{\ast}(\mathbf{r}-\mathbf{R}_{A}-\mathbf{e}_{\mu})\left[M_{R2,x}(\mathbf{r})\right.\cr&\left.+iM_{R2,y}(\mathbf{r})\right]w_{A\uparrow}(\mathbf{r}-\mathbf{R}_{A}),\cr
g_{B\mu}&=\int\mathrm{d}^{2}\mathbf{r}~{}w_{A\downarrow}^{\ast}(\mathbf{r}-\mathbf{R}_{B}-\mathbf{e}_{\mu})\left[M_{R2,x}(\mathbf{r})\right.\cr&\left.+iM_{R2,y}(\mathbf{r})+\right]w_{B\uparrow}(\mathbf{r}-\mathbf{R}_{B}),\end{split}$
(35)

and for next-nearest-neighbor spin-flip hopping terms caused by red Raman
fields $M_{2}(\mathbf{r})$:

$\begin{split}J_{A\mu\nu}&=\int\mathrm{d}^{2}\mathbf{r}~{}w_{A\downarrow}^{\ast}(\mathbf{r}-\mathbf{R}_{A}-\mathbf{e}_{\mu}+\mathbf{e}_{\nu})\left[M_{R1,x}(\mathbf{r})\right.\\\
&\left.+iM_{R1,y}(\mathbf{r})\right]w_{A\uparrow}(\mathbf{r}-\mathbf{R}_{A})\\\
J_{B\mu\nu}&=\int\mathrm{d}^{2}\mathbf{r}~{}w_{B\downarrow}^{\ast}(\mathbf{r}-\mathbf{R}_{B}-\mathbf{e}_{\mu}+\mathbf{e}_{\nu})\left[M_{R1,x}(\mathbf{r})\right.\\\
&\left.+iM_{R1,y}(\mathbf{r})\right]w_{B\uparrow}(\mathbf{r}-\mathbf{R}_{B})\end{split}$
(36)

where $J_{A\mu\nu}=-J_{B\mu\nu}$ since $C_{6}$ anti-symmetry of red Raman
field $M_{1}(\mathbf{r})$.

Considering the actual case where the overlap center of two Wannier-Stark
functions is not strictly located on the $C_{6}$ anti-symmetric line, we still
need to account for the influence of these effects on the tight-binding model.
This imply the red Raman fields also introduce nearest-neighbor spin-flip
terms, which exhibit staggered coupling, as shown in Fig. 5. The red and blue
arrows indicate the main contribution regions of nearest-neighbor spin-flip
hopping. Following the setting in the main text, by setting
$J_{1}=J_{2}=iJ_{I}$ and $g_{1}=-g_{2}=-ig$, the new terms cause by the
properties of the Wannier-Stark functions can be considered as a modification
of $g$, and this will not influence the form of the tight-binding model.

## References

  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Benalcazar  _et al._ (2017a) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Science 357, 61 (2017a).
  * Benalcazar  _et al._ (2017b) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B 96, 245115 (2017b).
  * Liu and Wakabayashi (2017) F. Liu and K. Wakabayashi, Phys. Rev. Lett. 118, 076803 (2017).
  * Langbehn  _et al._ (2017) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017).
  * Li  _et al._ (2017) L. Li, H. H. Yap, M. A. N. Araújo, and J. Gong, Phys. Rev. B 96, 235424 (2017).
  * Ezawa (2018a) M. Ezawa, Phys. Rev. Lett. 120, 026801 (2018a).
  * Geier  _et al._ (2018) M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Phys. Rev. B 97, 205135 (2018).
  * Sheng  _et al._ (2019) X.-L. Sheng, C. Chen, H. Liu, Z. Chen, Z.-M. Yu, Y. X. Zhao, and S. A. Yang, Phys. Rev. Lett. 123, 256402 (2019).
  * Trifunovic and Brouwer (2019) L. Trifunovic and P. W. Brouwer, Phys. Rev. X 9, 011012 (2019).
  * Călugăru  _et al._ (2019) D. Călugăru, V. Juričić, and B. Roy, Phys. Rev. B 99, 041301 (2019).
  * Park  _et al._ (2019) M. J. Park, Y. Kim, G. Y. Cho, and S. Lee, Phys. Rev. Lett. 123, 216803 (2019).
  * Pan  _et al._ (2019) X.-H. Pan, K.-J. Yang, L. Chen, G. Xu, C.-X. Liu, and X. Liu, Phys. Rev. Lett. 123, 156801 (2019).
  * Zhang  _et al._ (2020) R.-X. Zhang, F. Wu, and S. Das Sarma, Phys. Rev. Lett. 124, 136407 (2020).
  * Lee  _et al._ (2020) E. Lee, R. Kim, J. Ahn, and B.-J. Yang, npj Quantum Materials 5, 1 (2020).
  * Trifunovic and Brouwer (2021) L. Trifunovic and P. W. Brouwer, physica status solidi (b) 258, 2000090 (2021).
  * Fu (2011) L. Fu, Phys. Rev. Lett. 106, 106802 (2011).
  * Hsieh  _et al._ (2012) T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Nature communications 3, 982 (2012).
  * Chiu  _et al._ (2016) C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * Chiu  _et al._ (2013) C.-K. Chiu, H. Yao, and S. Ryu, Phys. Rev. B 88, 075142 (2013).
  * Cornfeld and Chapman (2019) E. Cornfeld and A. Chapman, Phys. Rev. B 99, 075105 (2019).
  * Ghorashi  _et al._ (2020) S. A. A. Ghorashi, T. L. Hughes, and E. Rossi, Phys. Rev. Lett. 125, 037001 (2020).
  * Kheirkhah  _et al._ (2020) M. Kheirkhah, Z. Yan, Y. Nagai, and F. Marsiglio, Phys. Rev. Lett. 125, 017001 (2020).
  * Li  _et al._ (2020) C.-A. Li, B. Fu, Z.-A. Hu, J. Li, and S.-Q. Shen, Phys. Rev. Lett. 125, 166801 (2020).
  * Wu  _et al._ (2020) X. Wu, W. A. Benalcazar, Y. Li, R. Thomale, C.-X. Liu, and J. Hu, Phys. Rev. X 10, 041014 (2020).
  * Ezawa (2020) M. Ezawa, Phys. Rev. B 102, 121405 (2020).
  * Asaga and Fukui (2020) K. Asaga and T. Fukui, Phys. Rev. B 102, 155102 (2020).
  * Claes and Hughes (2020) J. Claes and T. L. Hughes, Phys. Rev. B 102, 100203 (2020).
  * Khalaf  _et al._ (2021) E. Khalaf, W. A. Benalcazar, T. L. Hughes, and R. Queiroz, Phys. Rev. Res. 3, 013239 (2021).
  * Yang  _et al._ (2020) Y.-B. Yang, K. Li, L.-M. Duan, and Y. Xu, Phys. Rev. Res. 2, 033029 (2020).
  * Ezawa (2018b) M. Ezawa, Phys. Rev. B 98, 045125 (2018b).
  * Jia  _et al._ (2023) W. Jia, X.-C. Zhou, L. Zhang, L. Zhang, and X.-J. Liu, Phys. Rev. Res. 5, L022032 (2023).
  * Schindler  _et al._ (2018a) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Science Advances 4, eaat0346 (2018a).
  * Zhang  _et al._ (2013) F. Zhang, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 110, 046404 (2013).
  * Song  _et al._ (2017) Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * Peng  _et al._ (2017) Y. Peng, Y. Bao, and F. von Oppen, Phys. Rev. B 95, 235143 (2017).
  * Schindler  _et al._ (2018b) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, _et al._ , Nature physics 14, 918 (2018b).
  * Serra-Garcia  _et al._ (2018) M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature 555, 342 (2018).
  * Ni  _et al._ (2020) X. Ni, M. Li, M. Weiner, A. Alù, and A. B. Khanikaev, Nature communications 11, 2108 (2020).
  * Du  _et al._ (2022) J. Du, T. Li, X. Fan, Q. Zhang, and C. Qiu, Phys. Rev. Lett. 128, 224301 (2022).
  * Chen  _et al._ (2022) J.-H. Chen, Z.-Z. Yang, W.-J. Yang, A.-Y. Guan, X.-Y. Zou, and J.-C. Cheng, Applied Physics Letters 120, 253508 (2022).
  * Mittal  _et al._ (2019) S. Mittal, V. V. Orre, G. Zhu, M. A. Gorlach, A. Poddubny, and M. Hafezi, Nature Photonics 13, 692 (2019).
  * Imhof  _et al._ (2018) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Nature Physics 14, 925 (2018).
  * Peterson  _et al._ (2018) C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. Bahl, Nature 555, 346 (2018).
  * Zhang  _et al._ (2018) D.-W. Zhang, Y.-Q. Zhu, Y. X. Zhao, H. Yan, and S.-L. Zhu, Advances in Physics 67, 253 (2018).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
  * Lewenstein  _et al._ (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, and U. Sen, Advances in Physics 56, 243 (2007).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbene, Nature Physics 8, 267 (2012).
  * Schäfer  _et al._ (2020) F. Schäfer, T. Fukuhara, S. Sugawa, Y. Takasu, and Y. Takahashi, Nature Reviews Physics 2, 411 (2020).
  * Fetter (2009) A. L. Fetter, Rev. Mod. Phys. 81, 647 (2009).
  * Sørensen  _et al._ (2005) A. S. Sørensen, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 94, 086803 (2005).
  * Lim  _et al._ (2008) L.-K. Lim, C. M. Smith, and A. Hemmerich, Phys. Rev. Lett. 100, 130402 (2008).
  * Lim  _et al._ (2010) L.-K. Lim, A. Lazarides, A. Hemmerich, and C. Morais Smith, Phys. Rev. A 82, 013616 (2010).
  * Mei  _et al._ (2012) F. Mei, S.-L. Zhu, Z.-M. Zhang, C. H. Oh, and N. Goldman, Phy. Rev. A 85, 013638 (2012).
  * Dalibard  _et al._ (2011) J. Dalibard, F. Gerbier, G. Juzeliūnas, and P. Öhberg, Rev. Mod. Phys. 83, 1523 (2011).
  * Lin  _et al._ (2009a) Y.-J. Lin, R. L. Compton, K. Jiménez-García, J. V. Porto, and I. B. Spielman, Nature 462, 628 (2009a).
  * Lin  _et al._ (2011) Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83 (2011).
  * Zhu  _et al._ (2011) S.-L. Zhu, L.-B. Shao, Z. D. Wang, and L.-M. Duan, Phy. Rev. Lett. 106, 100404 (2011).
  * Lin  _et al._ (2009b) Y.-J. Lin, R. L. Compton, A. R. Perry, W. D. Phillips, J. V. Porto, and I. B. Spielman, Phys. Rev. Lett. 102, 130401 (2009b).
  * Aidelsburger  _et al._ (2011) M. Aidelsburger, M. Atala, S. Nascimbène, S. Trotzky, Y.-A. Chen, and I. Bloch, Phys. Rev. Lett. 107, 255301 (2011).
  * Osterloh  _et al._ (2005) K. Osterloh, M. Baig, L. Santos, P. Zoller, and M. Lewenstein, Phys. Rev. Lett. 95, 010403 (2005).
  * Li  _et al._ (2022) J.-Z. Li, C.-J. Zou, Y.-X. Du, Q.-X. Lv, W. Huang, Z.-T. Liang, D.-W. Zhang, H. Yan, S. Zhang, and S.-L. Zhu, Phys. Rev. Lett. 129, 220402 (2022).
  * Jaksch and Zoller (2003) D. Jaksch and P. Zoller, New Journal of Physics 5, 56 (2003).
  * Gerbier and Dalibard (2010) F. Gerbier and J. Dalibard, New Journal of Physics 12, 033007 (2010).
  * Goldman  _et al._ (2009) N. Goldman, A. Kubasiak, A. Bermudez, P. Gaspard, M. Lewenstein, and M. A. Martin-Delgado, Phys. Rev. Lett. 103, 035301 (2009).
  * Huang  _et al._ (2016) L. Huang, Z. Meng, P. Wang, P. Peng, S.-L. Zhang, L. Chen, D. Li, Q. Zhou, and J. Zhang, Nature Physics 12, 540 (2016).
  * Jotzu  _et al._ (2014) G. Jotzu, M. Messer, R. Desbuquois, M. Lebrat, T. Uehlinger, D. Greif, and T. Esslinger, Nature 515, 237 (2014).
  * Zhu  _et al._ (2006) S.-L. Zhu, H. Fu, C.-J. Wu, S.-C. Zhang, and L.-M. Duan, Phys. Rev. Lett. 97, 240401 (2006).
  * Goldman, N. and Gaspard, P. (2007) Goldman, N. and Gaspard, P., EPL 78, 60001 (2007).
  * Shao  _et al._ (2008) L. B. Shao, S.-L. Zhu, L. Sheng, D. Y. Xing, and Z. D. Wang, Phys. Rev. Lett. 101, 246810 (2008).
  * Wu (2008) C. Wu, Phys. Rev. Lett. 101, 186807 (2008).
  * Liang  _et al._ (2023) M.-C. Liang, Y.-D. Wei, L. Zhang, X.-J. Wang, H. Zhang, W.-W. Wang, W. Qi, X.-J. Liu, and X. Zhang, Phys. Rev. Res. 5, L012006 (2023).
  * Stanescu  _et al._ (2009) T. D. Stanescu, V. Galitski, J. Y. Vaishnav, C. W. Clark, and S. Das Sarma, Phys. Rev. A 79, 053639 (2009).
  * Liu  _et al._ (2007) X.-J. Liu, X. Liu, L. C. Kwek, and C. H. Oh, Phys. Rev. Lett. 98, 026602 (2007).
  * Beeler  _et al._ (2013) M. C. Beeler, R. A. Williams, K. Jiménez-García, L. J. LeBlanc, A. R. Perry, and I. B. Spielman, Nature 498, 201 (2013).
  * Alba  _et al._ (2011) E. Alba, X. Fernandez-Gonzalvo, J. Mur-Petit, J. K. Pachos, and J. J. Garcia-Ripoll, Phys. Rev. Lett. 107, 235301 (2011).
  * Béri and Cooper (2011) B. Béri and N. R. Cooper, Phys. Rev. Lett. 107, 145301 (2011).
  * Goldman  _et al._ (2010) N. Goldman, I. Satija, P. Nikolic, A. Bermudez, M. A. Martin-Delgado, M. Lewenstein, and I. B. Spielman, Phys. Rev. Lett. 105, 255302 (2010).
  * Wu  _et al._ (2016) Z. Wu, L. Zhang, W. Sun, X.-T. Xu, B.-Z. Wang, S.-C. Ji, Y. Deng, S. Chen, X.-J. Liu, and J.-W. Pan, Science 354, 83 (2016).
  * Wang  _et al._ (2018) B.-Z. Wang, Y.-H. Lu, W. Sun, S. Chen, Y. Deng, and X.-J. Liu, Phys. Rev. A 97, 011605 (2018).
  * Zeng  _et al._ (2019) C. Zeng, T. D. Stanescu, C. Zhang, V. W. Scarola, and S. Tewari, Phys. Rev. Lett. 123, 060402 (2019).
  * Wu  _et al._ (2023) Y.-B. Wu, Z. Zheng, X.-G. Qiu, L. Zhuang, G.-C. Guo, X.-B. Zou, and W.-M. Liu, Phys. Rev. A 107, 043304 (2023).
  * Huang and Liu (2020) B. Huang and W. V. Liu, Phys. Rev. Lett. 124, 216601 (2020).
  * Deng  _et al._ (2014) D.-L. Deng, S.-T. Wang, and L.-M. Duan, Phys. Rev. A 90, 041601 (2014).
  * Hauke  _et al._ (2014) P. Hauke, M. Lewenstein, and A. Eckardt, Phys. Rev. Lett. 113, 045303 (2014).
  * Fläschner  _et al._ (2016) N. Fläschner, B. S. Rem, M. Tarnowski, D. Vogel, D.-S. Lühmann, K. Sengstock, and C. Weitenberg, Science 352, 1091 (2016).
  * Yi  _et al._ (2023) C.-R. Yi, J. Yu, H. Yuan, R.-H. Jiao, Y.-M. Yang, X. Jiang, J.-Y. Zhang, S. Chen, and J.-W. Pan, Phys. Rev. Res. 5, L032016 (2023).
  * Bzdušek and Sigrist (2017) T. c. v. Bzdušek and M. Sigrist, Phys. Rev. B 96, 155105 (2017).
  * Ren  _et al._ (2020) Y. Ren, Z. Qiao, and Q. Niu, Phys. Rev. Lett. 124, 166804 (2020).
  * Lee  _et al._ (2007) D.-H. Lee, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. 99, 196805 (2007).
  * Sitte  _et al._ (2012) M. Sitte, A. Rosch, E. Altman, and L. Fritz, Phys. Rev. Lett. 108, 126807 (2012).
  * Yan  _et al._ (2018) Z. Yan, F. Song, and Z. Wang, Phys. Rev. Lett. 121, 096803 (2018).
  * Tan  _et al._ (2022) Y. Tan, Z.-H. Huang, and X.-J. Liu, Phys. Rev. B 105, L041105 (2022).
  * Jackiw and Rebbi (1976) R. Jackiw and C. Rebbi, Phys. Rev. D 13, 3398 (1976).
  * Zhu  _et al._ (2007) S.-L. Zhu, B. Wang, and L.-M. Duan, Phys. Rev. Lett. 98, 260402 (2007).
  * Duan  _et al._ (2003) L.-M. Duan, E. Demler, and M. D. Lukin, Phys. Rev. Lett. 91, 090402 (2003).
  * Soltan-Panahi  _et al._ (2011) P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick, W. Plenkers, G. Meineke, C. Becker, P. Windpassinger, M. Lewenstein, and K. Sengstock, Nature Physics 7, 434 (2011).
  * Lee  _et al._ (2009) K. L. Lee, B. Grémaud, R. Han, B.-G. Englert, and C. Miniatura, Phys. Rev. A 80, 043411 (2009).
  * Steck (2010) D. A. Steck, “Cesium d line data(v2.3.2),”  (2010).
  * Tiecke (2019) T. G. Tiecke, “Properties of potassium(v1.03),”  (2019).
  * Kennedy  _et al._ (2013) C. J. Kennedy, G. A. Siviloglou, H. Miyake, W. C. Burton, and W. Ketterle, Phys. Rev. Lett. 111, 225301 (2013).
  * Aidelsburger  _et al._ (2013) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Phys. Rev. Lett. 111, 185301 (2013).
  * Zhang  _et al._ (2016) D.-W. Zhang, Y. X. Zhao, R.-B. Liu, Z.-Y. Xue, S.-L. Zhu, and Z. D. Wang, Phy. Rev. A 93, 043617 (2016).