# Higher-order topological insulator with glide reflection symmetry

Shaoliang Zhang<EMAIL_ADDRESS>School of Physics, Institute of Quantum
Science and Engineering, Huazhong University of Science and Technology, Wuhan
430074, China

###### Abstract

The Benalcazar-Bernevig-Hughes(BBH) model is one of the most representative
higher-order topological insulators with different bulk-edge correspondence
comparing conventional topological states. In this work, we propose to realize
a couple-BBH model with glide reflection symmetry by coupling two
topologically distinct BBH models. Such a system gives rise to intriguing
properties beyond the traditional BBH model. A topological phase transition
from higher-order topological phase to a semimetal has been predicted in this
system. Furthermore, the band crossing points and nodal ring protected by
glide reflection symmetry has also been explored. Our work provides a
viewpoint to demonstrate the versatile phenomena induced by the interplay
between higher-order topology and nonsymmorphic symmetries.

##  I introduction

Two-dimensional Benalcazar-Bernevig-Hughes(BBH) modelBernevig1 ; Bernevig2 ,
as one of the most representative model of higher-order topological
insulator(HOTI), has attracted a lot of interests in condensed matter physics
because of its simple structure and beautiful topological properties.
Different with conventional topological statesmzhasan , there is no gapless
edge state when BBH model is in topological non-trivial phase no matter which
direction is chosen as open boundary. But if the open boundary is chosen as a
square, four topologically protected zero-energy corner states emerge on the
corners of the system. This new type of bulk-edge correspondence not only
extends the understanding of topological physicstrifunovic , but also inspires
new viewpoints to the application of topological mattersbanerjee ; cli .
Recently, the research about HOTI developed rapidly and a series of
theoretical models about HOTI have been proposedhotimodel1 ; hotimodel2 ;
hotimodel3 ; hotimodel4 ; hotimodel5 ; hotimodel6 ; hotimodel7 . Some of them
have been constructed experimentally in some synthetic quantum systems and
metamaterialsmetamaterial1 ; metamaterial2 ; metamaterial3 ; metamaterial4 ;
metamaterial5 ; metamaterial6 , and the gapless corner states has been
observed. The HOTI has also been observed in some condensed matter
systemshotiinbismuth1 ; hotiinbismuth2 .

The topology of higher-order topological states always depend on the
symmetriessymmetry1 ; symmetry2 ; symmetry3 ; symmetry4 ; symmetry5 ;
symmetry6 ; symmetry7 . For example, the bulk topological properties of BBH
model can be characterized by the quantized quadrupole moment, where the
quantization is protected by the chiral symmetry and $C_{\\_}2$
symmetryBernevig1 ; Bernevig2 . In the previous works, these symmetries are
always belonged to the point-group symmetries. In the research of symmetry
protected topological states, people also focus on the role of space-group
symmetries, especially the nonsymmorphic symmetries, which are the combination
of point-group symmetric operations with fractional translations of lattice.
The nonsymmorphic symmetry can lead to a lot of novel topological states
because it will induce additional band crossing, such as nodal points or nodal
linesglide1 ; glide2 ; glide3 ; glide4 ; glide5 ; glide6 ; glide7 ; glide8 .
So far, there are comprehensive and thorough researches about nonsymmorphic
symmetry protected topological states, but only a few of them involve the
interplay between higher-order topology and nonsymmorphic symmetriesglidehoti1
; glidehoti2 ; glidehoti3 .

In this work, we propose a coupled-BBH model which can be constructed in
ultracold atomic systems or some other synthetic quantum systems. It is
constituted by two BBH models coupled with each other. We consider the
ultracold atoms with two hyperfine spin states and each spin component is
confined in a spin-dependent two-dimensional lattice respectively, where the
lattice structure of different spin components has a fractional lattice
distance. By fine tuning the form of coupling between different spin
components, we can ensure the system has the glide reflection symmetry, which
is one of the nonsymmorphic symmetries. We found that, although the system is
a HOTI, the topological properties depend on the amplitude of coupling instead
of the ratio of intra and inter-cell tunneling in BBH model. Furthermore, the
glide reflection symmetry has significant effects to the properties of system.
(1) If the coupling term preserves the glide reflection symmetry and time-
reversal symmetry, and they commute with each other, the band structure has
two-fold degeneracy and two degenerate states has opposite eigenvalues of
glide operator respectively. (2) If the coupling term preserves the glide
reflection symmetry but breaks time-reversal symmetry, the two-fold degeneracy
will be lifted. But the glide reflection symmetry leads to additional band
crossing and changes the topological properties of system. With the increase
of coupling amplitude, the system changes from a HOTI to a semimetal. (3) If
the coupling term breaks glide reflection symmetry and time-reversal symmetry
individually but preserves their combination, which we call it preserves the
magnetic glide symmetry, the symmetry protected nodal ring emerges in the
system and the topological phase transition is also very different with above
two cases.

This work is organized as follows. In Section II, we introduce our model with
on-site inter-spin coupling and discuss its topological properties. In Section
III and IV, we consider two different types of inter-spin couplings, one
breaks the time-reversal symmetry but preserves glide reflection symmetry; the
other breaks both but preserves their combination, i.e., the magnetic glide
symmetry. In Section V, the experimental setup to realize the discussed models
in ultracold atomic systems has been proposed. Section VI is devoted to a
discussion.

##  II Model and topological properties

We consider a square lattice with two different spin components. The lattice
structure of each spin component is similar with BBH modelBernevig1 ;
Bernevig2 , as shown in fig.1, where the lattice structure for different spin
components has a displacement along the diagonal line of $x$ and $y$
directions. An on-site coupling between different spin components has been
added and the tight-binding Hamiltonian can be written as

$H=H_{\\_}\uparrow+H_{\\_}\downarrow+H_{\\_}\mathrm{couple}$ (1)

where $H_{\\_}\uparrow$ and $H_{\\_}\downarrow$ are the Hamiltonian of spin-up
and spin-down component respectively as

$\begin{split}H_{\\_}\uparrow=&\lambda\sum_{\\_}{mn}\big{(}a^{\dagger}_{\\_}{mn,\uparrow}b_{\\_}{mn,\uparrow}-c^{\dagger}_{\\_}{mn,\uparrow}d_{\\_}{mn,\uparrow}\big{)}\\\
+&\gamma\sum_{\\_}{mn}\big{(}b^{\dagger}_{\\_}{mn,\uparrow}a_{\\_}{(m+1)n,\uparrow}-d^{\dagger}_{\\_}{mn,\uparrow}c_{\\_}{(m+1)n,\uparrow}\big{)}\\\
+&\lambda\sum_{\\_}{mn}\big{(}a^{\dagger}_{\\_}{mn,\uparrow}c_{\\_}{mn,\uparrow}+b^{\dagger}_{\\_}{mn,\uparrow}d_{\\_}{mn,\uparrow}\big{)}\\\
+&\gamma\sum_{\\_}{mn}\big{(}c^{\dagger}_{\\_}{mn,\uparrow}a_{\\_}{(m+1)n,\uparrow}+d^{\dagger}_{\\_}{mn,\uparrow}b_{\\_}{(m+1)n,\uparrow}\big{)}+h.c.\end{split}$
(2)
$\begin{split}H_{\\_}\downarrow=&\gamma\sum_{\\_}{mn}\big{(}c^{\dagger}_{\\_}{mn,\downarrow}d_{\\_}{mn,\downarrow}-a^{\dagger}_{\\_}{mn,\downarrow}b_{\\_}{mn,\downarrow}\big{)}\\\
+&\lambda\sum_{\\_}{mn}\big{(}d^{\dagger}_{\\_}{mn,\downarrow}c_{\\_}{(m+1)n,\downarrow}-b^{\dagger}_{\\_}{mn,\downarrow}a_{\\_}{(m+1)n,\downarrow}\big{)}\\\
+&\gamma\sum_{\\_}{mn}\big{(}a^{\dagger}_{\\_}{mn,\downarrow}c_{\\_}{mn,\downarrow}+b^{\dagger}_{\\_}{mn,\downarrow}d_{\\_}{mn,\downarrow}\big{)}\\\
+&\lambda\sum_{\\_}{mn}\big{(}c^{\dagger}_{\\_}{mn,\downarrow}a_{\\_}{(m+1)n,\downarrow}+d^{\dagger}_{\\_}{mn,\downarrow}b_{\\_}{(m+1)n,\downarrow}\big{)}+h.c.\end{split}$
(3)

where $\lambda(\gamma)$ and $\gamma(\lambda)$ are amplitudes of intra- and
inter-cell tunneling for spin-up (down) component respectively.
$H_{\\_}\mathrm{couple}$ is the coupling between different spin components as

$\begin{split}H_{\\_}\mathrm{couple}=&2\Omega\sum_{\\_}{mn}\big{(}a^{\dagger}_{\\_}{mn,\uparrow}a_{\\_}{mn,\downarrow}+b^{\dagger}_{\\_}{mn,\uparrow}b_{\\_}{mn,\downarrow}\\\
&+c^{\dagger}_{\\_}{mn,\uparrow}c_{\\_}{mn,\downarrow}+d^{\dagger}_{\\_}{mn,\uparrow}d_{\\_}{mn,\downarrow}\big{)}+h.c.\end{split}$
(4)

where $2\Omega$ corresponding to the coupling amplitude.

Figure 1: The scheme of lattice structure. Red and blue dots correspond to
different spin components. Thick and thin lines correspond to the amplitude of
nearest-neighbor tunneling $\lambda$ and $\gamma$ respectively. Dotted line
corresponds to the tunneling amplitude has an additional minus sign. Black
double arrow lines correspond to the on-site coupling between different spin
components. The inset on right shows the first Brillouin zone and some high
symmetric points.

When $\Omega=0$, different spin components are independent and the topological
properties of the system is clear. If $\lambda<\gamma$, Hamiltonian
$H_{\\_}\uparrow$ is topological non-trivial but Hamiltonian
$H_{\\_}\downarrow$ is topological trivial. If $\lambda>\gamma$, the answer is
opposite. So no matter in which case, the whole system is topological non-
trivial. With a square open boundary, there will always be a zero energy state
localized on each corner. In this work, we will only consider the
$\lambda<\gamma$ case and set the length of square unit cell as $1$ for
simplicity.

When $\Omega\neq 0$, the Hamiltonian (1) is invariant with combination of spin
flip $\uparrow\leftrightarrow\downarrow$ and the translation of the distance
$1/\sqrt{2}$ along $x\pm y$ directions, which means the system has the glide
reflection symmetry. With Fourier transformation, the above Hamiltonian (1)
can be rewritten in momentum space as $H=\sum_{\\_}{\bf
k}\Psi^{\dagger}_{\\_}{\bf k}H_{\\_}{\bf k}\Psi_{\\_}{\bf k}$ where
$\Psi_{\\_}{\bf k}=[a_{\\_}{{\bf k},\uparrow},b_{\\_}{{\bf
k},\uparrow},c_{\\_}{{\bf k},\uparrow},d_{\\_}{{\bf k},\uparrow},a_{\\_}{{\bf
k},\downarrow},b_{\\_}{{\bf k},\downarrow},c_{\\_}{{\bf
k},\downarrow},d_{\\_}{{\bf k},\downarrow}]^{\mathrm{T}}$ and

$\begin{split}H_{\\_}{\bf k}=&-\Gamma(1-\cos
k_{\\_}x)\tau_{\\_}z\sigma_{\\_}x+\Delta(1+\cos
k_{\\_}x)s_{\\_}z\tau_{\\_}z\sigma_{\\_}x\\\ &+\Gamma\sin
k_{\\_}x\tau_{\\_}z\sigma_{\\_}y+\Delta\sin
k_{\\_}xs_{\\_}z\tau_{\\_}z\sigma_{\\_}y\\\ &+\Delta(1+\cos
k_{\\_}y)\tau_{\\_}x-\Gamma(1-\cos k_{\\_}y)s_{\\_}z\tau_{\\_}x\\\
&+\Delta\sin k_{\\_}y\tau_{\\_}y+\Gamma\sin k_{\\_}ys_{\\_}z\tau_{\\_}y+\Omega
s_{\\_}x\end{split}$ (5)

where $\Gamma=(\lambda-\gamma)/2$ and $\Delta=(\lambda+\gamma)/2$.
$\tau_{\\_}i$, $\sigma_{\\_}i$ and $s_{\\_}i$ ($i=x,y,z$) are Pauli matrices
which corresponding to orbital ($\bm{\tau}$, $\bm{\sigma}$) and spin
($\bm{s}$) degrees of freedom respectively. It can be separated into two parts
as $H_{\\_}{\bf k}=H_{\\_}{d,{\bf k}}+H_{\\_}\mathrm{couple}$. The block
diagonalized matrix $H_{\\_}{d,{\bf k}}=\mathrm{diag}\big{[}H_{\\_}{{\bf
k},\uparrow},H_{\\_}{{\bf k},\downarrow}\big{]}$ corresponds to Hamiltonian of
two different spin components where $H_{\\_}{{\bf k},\uparrow}$ and
$H_{\\_}{{\bf k},\downarrow}$ are Hamiltonian of spin-up and spin-down
components respectively, which can be expressed separately as

$\begin{split}H_{\\_}{{\bf k},\uparrow}=&\big{(}\lambda+\gamma\cos
k_{\\_}x\big{)}\tau_{\\_}z\sigma_{\\_}x+\gamma\sin
k_{\\_}x\tau_{\\_}z\sigma_{\\_}y\\\ &+\big{(}\lambda+\gamma\cos
k_{\\_}y\big{)}\tau_{\\_}x+\gamma\sin k_{\\_}y\tau_{\\_}y\end{split}$ (6)
$\begin{split}H_{\\_}{{\bf k},\downarrow}=&-\big{(}\gamma+\lambda\cos
k_{\\_}x\big{)}\tau_{\\_}z\sigma_{\\_}x-\lambda\sin
k_{\\_}x\tau_{\\_}z\sigma_{\\_}y\\\ &+\big{(}\gamma+\lambda\cos
k_{\\_}y\big{)}\tau_{\\_}x+\lambda\sin k_{\\_}y\tau_{\\_}y\end{split}$ (7)

$H_{\\_}\mathrm{couple}=\Omega s_{\\_}x$ corresponds to the coupling between
different spin components. The glide operator can also be written in momentum
space as $G_{\\_}{\bf k}=s_{\\_}x\tilde{G}_{\\_}{\bf k}$ where

$\tilde{G}_{\\_}{\bf k}=\left(\begin{array}[]{cc}0&1\\\
e^{ik_{\\_}y}&0\end{array}\right)\otimes\left(\begin{array}[]{cc}0&1\\\
e^{ik_{\\_}x}&0\end{array}\right)$ (8)

Here we just consider the glide operation along $x+y$ direction and the
corresponding operation along $x-y$ direction also has the similar result. It
can be verified that the Hamiltonian (5) satisfies the relation

$G^{\dagger}_{\\_}{\bf k}H_{\\_}{\bf k}G_{\\_}{\bf k}=H_{\\_}{\bf k}$ (9)

For each block of matrix $H_{\\_}{d,{\bf k}}$, we also have the relation

$\tilde{G}^{\dagger}_{\\_}{\bf k}H_{\\_}{{\bf k},\uparrow}\tilde{G}_{\\_}{\bf
k}=H_{\\_}{{\bf k},\downarrow}$ (10)

Two consecutive glide operations correspond to one lattice translation along
$x+y$ direction as

$G^{2}_{\\_}{\bf k}=e^{i(k_{\\_}x+k_{\\_}y)}\mathcal{I}$ (11)

where $\mathcal{I}$ is identity matrix, which means that the eigenvalues of
glide operator should be $\pm e^{i(k_{\\_}x+k_{\\_}y)/2}$.

Figure 2: The band structure of system with on-site coupling when
$\lambda/\gamma=0.4$ for different coupling amplitude. (a)
$\Omega/\gamma=0.4$, (b) $\Omega/\gamma=0.76$, (c) $\Omega/\gamma=0.9$.

The band structure is shown in fig.2. The system has time-reversal symmetry
and chiral symmetry with symmetric operator as $\mathcal{T}=\mathcal{K}$ and
$\mathcal{C}=s_{\\_}z\tau_{\\_}z\sigma_{\\_}z$, where $\mathcal{K}$
corresponds to the Hermitian conjugation. Similar with BBH modelBernevig1 ;
Bernevig2 , the system (5) has reflection symmetry along $x$ and $y$
directions with symmetric operator as $M_{\\_}x=\sigma_{\\_}x$ and
$M_{\\_}y=\tau_{\\_}x\sigma_{\\_}z$ respectively. It also has $C_{\\_}2$
rotation symmetry with symmetric operator as
$C_{\\_}2=\tau_{\\_}x\sigma_{\\_}y$, which means that
$C^{-1}_{\\_}2H_{\\_}{(-k_{\\_}x,-k_{\\_}y)}C_{\\_}2=H_{\\_}{(k_{\\_}x,k_{\\_}y)}$.
These symmetries protect the quantization of quadrupole moment and our system
is still a second-order topological insulator. But the topological properties
of the Hamiltonian (5) depends on the coupling amplitude between different
spin components. With the increase of $\Omega$, the system becomes from
topological non-trivial phase to topological trivial phase. As shown in
fig.2(b), the critical point of topological phase transition is at

$\Omega_{\\_}\mathrm{cri}=\sqrt{\frac{\lambda^{2}+\gamma^{2}}{2}}$ (12)

Numerical calculation shows that there are four zero-energy corner states with
square open boundaries when $\Omega<\Omega_{\\_}\mathrm{cri}$. In this case,
if we consider the bulk topology, four Wannier bands are all gapped and the
corresponding Wannier sector polarization are all quantized and two of them
are non-trivial and the other two are trivialBernevig1 ; Bernevig2 . After
passing through the critical point as $\Omega>\Omega_{\\_}\mathrm{cri}$, the
Wannier sector polarizations of all four Wannier bands are trivial and there
is no zero-energy corner states with square open boundaries anymore.

In the whole range of the coupling amplitude $\Omega$, there is four-fold
degeneracy at $Y=(0,\pi)$ point, as shown in fig.2. We should point out that
the degeneracy comes from the form of Hamiltonian (5) at $Y$ point only
contains three terms $s_{\\_}z\tau_{\\_}z\sigma_{\\_}x$, $s_{\\_}z\tau_{\\_}x$
and $s_{\\_}x$, whose square are all identity matrix and anti-commute with
each other. It doesn’t relate to the glide reflection symmetry because this
degeneracy can be lifted by adding an additional energy detuning
$\delta\tau_{\\_}z\sigma_{\\_}z$ which preserve the glide reflection symmetry.
Actually, the time-reversal symmetry and glide reflection symmetry commute
with each other as $\big{[}\mathcal{T},G_{\\_}{\bf k}\big{]}=0$ in this
system. We can find an operator $\Theta=i\tau_{\\_}x\sigma_{\\_}y\mathcal{K}$
which is commute with the Hamiltonian (5) and preserves the Kramers degeneracy
in the whole Brillouin zone(BZ) due to $\Theta^{2}=-1$. This operator and
glide operator $G_{\\_}{\bf k}$ also has the relation as

$\Theta^{-1}G_{\\_}{\bf k}\Theta=-e^{-i(k_{\\_}x+k_{\\_}y)}G_{\\_}{\bf k}$
(13)

which means that if state $|\psi\rangle$ is the common eigenstate of
Hamiltonian (5) and glide operator $G_{\\_}{\bf k}$ with eigenvalue $E$ and
$e^{i(k_{\\_}x+k_{\\_}y)/2}$ respectively, then $\Theta|\psi\rangle$ is
degenerate with $|\psi\rangle$ and also the eigenstate of $G_{\\_}{\bf k}$
with opposite eigenvalue as $-e^{i(k_{\\_}x+k_{\\_}y)/2}$. It implies that the
eigenstates of glide operator with opposite eigenvalues always degenerate with
each other in the whole BZ.

##  III The breaking of time-reversal symmetry

Instead of the on-site coupling between different spin components, one can add
nearest-neighbor tunneling with additional phases along clockwise
($a_{\\_}\uparrow$,$c_{\\_}\uparrow$) and anti-clockwise
($b_{\\_}\uparrow$,$d_{\\_}\uparrow$) directions respectively, as shown in
fig.3, which breaking the time-reversal symmetry. The Hamiltonian of two
different spin components $H_{\\_}{d,{\bf k}}$ has the same form with (5) and
the coupling term changes into

$\begin{split}H_{\\_}\mathrm{couple,{\bf k}}&=\Omega(1-\cos
k_{\\_}x)s_{\\_}x\sigma_{\\_}y+\Omega\sin k_{\\_}xs_{\\_}x\sigma_{\\_}x\\\
+&\Omega(1-\cos k_{\\_}y)s_{\\_}y\tau_{\\_}y\sigma_{\\_}z+\Omega\sin
k_{\\_}ys_{\\_}y\tau_{\\_}x\sigma_{\\_}z\end{split}$ (14)

Figure 3: Left: The scheme of lattice structure with nearest-neighbor
coupling. Purple arrow lines correspond to the nearest-neighbor coupling
between different spin components. Right: The tunneling from a spin-up site to
its nearest-neighbor spin-down sites along different directions with different
additional phases.

It is easy to verify that the system still has the glide reflection symmetry
with the same glide operator as the above discussion.The chiral symmetry and
$C_{\\_}2$ symmetry are also preserved with corresponding symmetric operators
as $\mathcal{C}=\tau_{\\_}z\sigma_{\\_}z$ and
$C_{\\_}2=\tau_{\\_}x\sigma_{\\_}y$, which means that the system is still a
second-order topological insulator. The band structure is shown in fig.4. This
type of coupling also preserve the reflection symmetries along $x$ and $y$
directions with the corresponding symmetric operators change into
$M_{\\_}x=s_{\\_}z\sigma_{\\_}x$ and
$M_{\\_}y=s_{\\_}z\tau_{\\_}x\sigma_{\\_}z$ respectively. In most area of the
BZ, two-fold degeneracy has been lifted. But along some high-symmetric lines,
such as along $k_{\\_}x=0$ and $k_{\\_}y=0,\pi$ lines, the band structure
still has two-fold degeneracy. We can define another operator
$\Theta^{\prime}=is_{\\_}z\tau_{\\_}x\sigma_{\\_}y\mathcal{K}$ which is
commute with original Hamiltonian (5) in Section II when $\Omega=0$. It also
commute with those $k_{\\_}x$ dependent terms in the coupling Hamiltonian
(14). So along $k_{\\_}y=0$ line, $\Theta^{\prime}$ is commute with the
Hamiltonian of the whole system and the band has Kramers degeneracy. Along
$k_{\\_}x=0$ line, the operator $\Theta=i\tau_{\\_}x\sigma_{\\_}y\mathcal{K}$
can be used to explain the degeneracy. Along $k_{\\_}y=\pi$ line, the Kramers
degeneracy can not be explained easily because the corresponding symmetric
operator should be ${\bf k}$-dependent.

Figure 4: Upper: The band structure of system with nearest-neighbor coupling
when $\lambda/\gamma=0.4$ for different coupling amplitudes. (a)
$\Omega/\gamma=0.3$, (b) $\Omega/\gamma=0.48$, (c) $\Omega/\gamma=0.6$. Lower:
The band structure with open boundary along $x$ direction with corresponding
system parameters.

The most interesting phenomena is along $k_{\\_}x=\pi$, where the two-fold
degeneracy is lifted. With the increase of coupling amplitude, there are two
band crossing points with quadratic dispersions along this line at critical
point and the system has a topological phase transition, as shown in fig.4(b)
(it only shows the $k_{\\_}y>0$ one). After passing through the critical
point, these two band crossing points separate into four band crossing points
with linear dispersions and the system becomes a semimetal, as shown in
fig.4(c). To explain these transition analytically, one can consider an ${\bf
k}$-dependent unitary transformation as $U_{\\_}{\bf
k}=\mathrm{diag}\big{[}\mathcal{I}_{\\_}{4\times
4},-ie^{-ik_{\\_}y/2}\tilde{G}_{\\_}{\bf k}\big{]}$ where
$H^{\prime}_{\\_}{\bf k}=U^{\dagger}_{\\_}{\bf k}H_{\\_}{k}U_{\\_}{\bf k}$.
With this transformation, from relation (10), two diagonal blocks of
$H^{\prime}_{\\_}{\bf k}$ become the same and the form of whole Hamiltonian
along $k_{\\_}x=\pi$ line can be written as

$\begin{split}H^{\prime}_{\\_}{k_{\\_}x=\pi}&=(\lambda-\gamma)\tau_{\\_}z\sigma_{\\_}x+(\lambda+\gamma\cos
k_{\\_}y)\tau_{\\_}x+\gamma\sin k_{\\_}y\tau_{\\_}y\\\
-&2\Omega\sin\frac{k_{\\_}y}{2}s_{\\_}x\sigma_{\\_}x+2\Omega\cos\frac{k_{\\_}y}{2}s_{\\_}x\tau_{\\_}x+2\Omega\sin\frac{k_{\\_}y}{2}s_{\\_}x\tau_{\\_}y\end{split}$
(15)

which can be diagonalized analytically to get the energy spectrum as

$E_{\\_}{k_{\\_}x=\pi}=\pm
2\Omega\sin\frac{k_{\\_}y}{2}\pm\sqrt{(\lambda-\gamma)^{2}+\Big{|}\lambda+\gamma
e^{ik_{\\_}y}\pm 2\Omega e^{\frac{ik_{\\_}y}{2}}\Big{|}^{2}}$ (16)

Then the critical point can be estimated as

$\Omega_{\\_}\mathrm{cri}=\sqrt{\frac{2\lambda\gamma(\lambda-\gamma)^{2}}{6\lambda\gamma-\lambda^{2}-\gamma^{2}}}$
(17)

We can also choose the open boundary along $x$ direction. As shown in
fig.4(f), the zero-energy edge states emerge in the restricted region between
two band crossing points in the semimetal phase, which is similar with the
bearded and zigzag edges of graphenegraphene . An interesting property is that
the relation (17) is meaningful only when the value under the square root is
positive, which means that the topological phase transition can only happen
when $(3-2\sqrt{2})\gamma<\lambda<\gamma$. If $\lambda\ll\gamma$, there is no
topological phase transition and the system is always a non-trivial HOTI with
the increase of coupling amplitude. Some details of this case is discussed in
Appendix A.

Figure 5: (a) The band structure for $k_{\\_}x=\pi$. Red and blue lines
correspond to the sign as $+$ and $-$ before the eigenvalue of glide operator.
(b) The gap between the second and third bands. The color from dark to light
correspond to the increase of gap energy. The darkest part corresponds to the
band crossing points. Two red dots in (a) and (b) correspond to each other.
The tunneling and coupling amplitude we choose in (a) and (b) are
$\lambda/\gamma=0.4$ and $\Omega/\gamma=0.3$ respectively.

In most area of BZ, Kramers degeneracy is lifted and glide reflection symmetry
induces additional band crossing. As shown in fig.5(a), different bands
correspond to different eigenvalues of the glide operator. Red and blue lines
in fig.5(a) correspond to the sign as $+$ and $-$ before the eigenvalues of
glide operator respectively. For an arbitrary $k_{\\_}y$,
$e^{i(k_{\\_}x+k_{\\_}y)/2}$ will return to the original value after a
transition $k_{\\_}y\rightarrow k_{\\_}y+4\pi$, which passing through twice of
BZ. So there should be at least one band crossing point in the first BZ, as
the red dots in fig.5. Furthermore, the band crossing points protected by the
glide reflection symmetry is robust to the change of system parameters such as
coupling amplitude. As shown in fig.5(b), the band crossing points between the
second and third bands come into a nodal line whose form has no qualitative
changes in distinct topological phases.

##  IV Magnetic glide symmetry

In this section, we consider another type of on-site coupling with the
corresponding coupling Hamiltonian as

$\begin{split}H_{\\_}\mathrm{couple}=&2\Omega\sum_{\\_}{mn}\big{(}e^{i\varphi}a^{\dagger}_{\\_}{mn,\uparrow}a_{\\_}{mn,\downarrow}+e^{-i\varphi}b^{\dagger}_{\\_}{mn,\uparrow}b_{\\_}{mn,\downarrow}\\\
&+e^{-i\varphi}c^{\dagger}_{\\_}{mn,\uparrow}c_{\\_}{mn,\downarrow}+e^{i\varphi}d^{\dagger}_{\\_}{mn,\uparrow}d_{\\_}{mn,\downarrow}\big{)}+h.c.\end{split}$
(18)

where $\varphi$ is an additional phase in the coupling. With the Fourier
transformation, it can be expressed in momentum space as

$H_{\\_}\mathrm{couple}=2\Omega\cos\varphi s_{\\_}x+2\Omega\sin\varphi
s_{\\_}y\tau_{\\_}z\sigma_{\\_}z$ (19)

When $\varphi=0,\pi/2,\pi$, this coupling term is equivalent to the case we
discussed in Section II, which preserve both the glide reflection symmetry and
time-reversal symmetry. But when $0<\varphi<\pi/2$ or $\pi<\varphi<\pi$, both
of these two symmetries has been broken, only their combination $G_{\\_}{\bf
k}\mathcal{T}$ is preserved. This symmetry is one of the nonsymmorphic
magnetic symmetriesglide3  and we call it magnetic glide symmetry. We can also
find that the system still has the chiral symmetry and $C_{\\_}2$ symmetry
with the same symmetric operators as in Section II, which means that the
system is still a HOTI. This type of coupling break the reflection symmetries
along $x$ and $y$ directions, but with unitary operators
$\tilde{M}_{\\_}x=\tau_{\\_}x\sigma_{\\_}z$ and
$\tilde{M}_{\\_}y=\sigma_{\\_}x$, we have the relations
$\tilde{M}^{-1}_{\\_}xH_{\\_}{(-k_{\\_}x,k_{\\_}y)}\tilde{M}_{\\_}x=H^{*}_{\\_}{(k_{\\_}x,k_{\\_}y)}$
and
$\tilde{M}^{-1}_{\\_}yH_{\\_}{(k_{\\_}x,-k_{\\_}y)}\tilde{M}_{\\_}y=H^{*}_{\\_}{(k_{\\_}x,k_{\\_}y)}$
respectively, where $H^{*}_{\\_}{(k_{\\_}x,k_{\\_}y)}$ corresponds to the
Hermitian conjugation of Hamiltonian $H_{\\_}{(k_{\\_}x,k_{\\_}y)}$. So the
band structure is still symmetric along $k_{\\_}x$ and $k_{\\_}y$ directions.

The band structure is shown in fig.6(a)-(c). With the increase of coupling
amplitude, the system has three different topological phases, which is very
different with above two cases. Numerical calculation shows that there are two
critical points in the system as $\Omega_{\\_}{1,\mathrm{cri}}$ and
$\Omega_{\\_}{2,\mathrm{cri}}$ which depend on the value of $\varphi$ and have
no explicit expressions. When the coupling amplitude
$\Omega<\Omega_{\\_}{1,\mathrm{cri}}$, the system is a non-trivial HOTI. When
$\Omega=\Omega_{\\_}{1,\mathrm{cri}}$, the system translates into a semimetal
with two band crossing points with quadrupole dispersion. After passing
through this critical point, these two band crossing points separate into four
band crossing points with linear dispersions and moving in the first BZ with
the increase of $\Omega$. When $\Omega$ approaches to
$\Omega_{\\_}{2,\mathrm{cri}}$, the four band crossing points merge into two
crossing points again. After passing through $\Omega_{\\_}{2,\mathrm{cri}}$,
band gap will be reopened and the system becomes topological trivial. At two
sides when $\Omega<\Omega_{\\_}{1,\mathrm{cri}}$ and
$\Omega>\Omega_{\\_}{2,\mathrm{cri}}$, the topology of bulk and corner are
similar with the case we discussed in Section II.

Figure 6: (a)-(c) correspond to the band structure when $\lambda/\gamma=0.4$
with different coupling amplitude $\Omega$ where $\varphi=\pi/3$. (a)
$\Omega/\gamma=0.3$, (b) $\Omega/\gamma=0.48$, (c) $\Omega/\gamma=0.6$. The
red dots correspond to the band crossing points protected by the magnetic
glide symmetry. (d)-(e) corresponds to the gap between the second and third
bands with different $\varphi$ when the coupling amplitude is fixed as
$\Omega/\gamma=0.3$. The color from dark to light correspond to the increase
of gap energy. The darkest part corresponds to the band crossing points. (d)
$\varphi=\pi/6$, (e) $\varphi=\pi/4$, (f) $\varphi=\pi/3$. The red dots in (f)
correspond to the positions of two band crossing points along $k_{\\_}x=\pi$
and $k_{\\_}y=0$ lines in BZ.

Because the breaking of time-reversal symmetry, two-fold degeneracy has been
lifted in most area of BZ. But the magnetic glide symmetry also lead to
additional band crossing points between the second and third bands, which
forms a nodal ring. As shown in fig.6 (d)-(f), the form of nodal ring has no
qualitative changes in different topological phases with the increase of
$\Omega$. It only depends on the value of the phase of coupling $\varphi$.

##  V Scheme of experimental setup

Exploring the models we discussed above in condensed matter systems directly
requires more in-depth research. But the synthetic quantum systems provide
convenient platforms for simulating novel topological states. For example, we
can try to simulate the coupled-BBH model in ultracold atomic systems with the
following three steps. The first step is simulating BBH model using two-
dimensional superlattice. The second step is constructing a spin-dependent
lattice where lattice structure of different spin component has a displacement
along particular direction. The third step is adding the coupling using
microwave or laser fields.

Step 1. Consider the ultracold atoms confined in two-dimensional superlattice,
which is constructed by two types of counter-propagating lasers with long and
short wavelength $4\pi/k_{\\_}0$ and $2\pi/k_{\\_}0$ respectively ($k_{\\_}0$
is the wave vector of laser with short wavelength) along both $x$ and $y$
directions, as shown in fig.7. The lattice structure can be written as
$V(x,y)=V(x)+V(y)$ where

$\begin{split}V(x)=&V_{\\_}{sx}\sin^{2}(k_{\\_}0x)+V_{\\_}{lx}\sin(k_{\\_}0x+\theta)\\\
V(y)=&V_{\\_}{sy}\sin^{2}(k_{\\_}0y)+V_{\\_}{ly}\sin(k_{\\_}0y)\end{split}$
(20)

where $0<\theta<\pi/2$ is the relative phase between two lasers with different
wavelength propagating along $x$ direction. The unit cell of the two-
dimensional superlattice is a square which including four neighbor sites.
Along $y$ direction, the difference between intra- and inter-cell nearest-
neighbor tunneling lead to a Su-Schrieffer-Heeger modelssh ; zakphase . But as
shown in fig.7(b), the lattice structure along $x$ direction is asymmetric
because of the non-zero relative phase $\theta$, the nearest-neighbor
tunneling along $x$ direction is forbidden by the energy detuning between two
neighbor sites. By adding two additional Raman lasers with the frequency
difference matching this energy detuning exactly, the Raman process leads to a
two-photon assisted tunneling with an additional phase, which corresponding to
adding an effective staggered flux in the superlatticefluxlattice . The
asymmetric of the lattice structure also induce the difference between
intracell and intercell tunneling amplitude. Setting the wavelength of two
Raman lasers equal to the wave length of short laser, the flux can be set as
$\pi$, which corresponding to the BBH model. The amplitude of tunneling along
$x$ direction can be fine tuned by choosing the Rabi frequencies of Raman
lasers. But we should point out that the difference of tunneling amplitudes
along $x$ and $y$ directions will not break the symmetries discussed in our
work.

Figure 7: Scheme of experimental setup. (a) The lattice structure for
different spin components. The spin-up and spin-down are shown by upper and
lower surfaces respectively. The red and blue lines correspond to the lattice
potential along $x$ and $y$ directions for different spin components
respectively. Red and blue arrows are the Raman lasers with wave vectors and
frequencies as ${\bf k}_{\\_}1,\omega_{\\_}1$ and ${\bf
k}_{\\_}2,\omega_{\\_}2$ respectively. (b) The lattice structure along $x$
direction for spin-up (upper) and spin-down (lower) respectively. The photon-
assisted nearest-neighbor tunneling is induced by Raman process. Two dotted
lines with double arrow correspond to the microwave fields coupling different
spin components on different types of sublattices respectively.

Step 2. There are a lot of mature experimental technics to realize the
construction of spin-dependent lattices in ultracold atomic
systemsspinlattice1 ; spinlattice2 ; spinlattice3 . We can also use two types
of lasers with the same wave vector but different frequencies to confine two
different hyperfine spin states of the same type of atoms individually. If the
lattice potential of spin-up component is (20), the lattice potential of spin-
down component can be written as $V^{\prime}(x,y)=V^{\prime}(x)+V^{\prime}(y)$
where
$V^{\prime}(x)=V_{\\_}{sx}\sin^{2}(k_{\\_}0x)-V_{\\_}{lx}\sin(k_{\\_}0x+\theta)$
and $V^{\prime}(y)=V_{\\_}{sy}\sin^{2}(k_{\\_}0y)-V_{\\_}{ly}\sin(k_{\\_}0y)$,
where the lattice potential with long wavelength has an opposite sign. It is
the same with the lattice potential of spin-up but a shift of half lattice
spacing along both $x$ and $y$ direction, which can be realized by a phase
shift between these two types of lasers. Then with the same Raman process on
two different types of superlattices, the lattice structure shown in fig.1 can
be constructed.

Step 3. We have to point out that it is not easy to construct the type of
coupling as discussed in Section III. Because it is not a Rashba-type spin-
orbit couplingsoc2d1 ; soc2d2 ; soc2d3  but a nearest-neighbor coupling with
additional phase along clockwise ($a_{\\_}\uparrow,c_{\\_}\uparrow$) and anti-
clockwise ($d_{\\_}{\uparrow},b_{\\_}\uparrow$) directions respectively. But
the realization of on-site coupling discussed in Section II and IV is not
difficult. As shown in fig.7(b), because the asymmetric lattice potential
along $x$ direction, the energy detuning between different spin components on
$a(c)$ and $b(d)$ sublattices are different. So two microwave fields with
different frequencies need to be added to couple different spin components on
$a(c)$ and $b(d)$ sublattices respectively. The frequencies of two microwave
fields should match the energy detuning of two different types of sublattices
respectively. By locking the relative phase as zero between these two
microwave fields, the on-site coupling in Section II can be realized.

The realization of the on-site coupling in Section IV is more complex. We need
to construct the coupling with the form as
$\Omega_{\\_}1s_{\\_}y\pm\Omega_{\\_}2\cos(k_{\\_}0x)s_{\\_}x$ where
$\varphi\approx\tan^{-1}(\Omega_{\\_}2/\Omega_{\\_}1)$, which can be separated
into two parts. The first constant part with the form $\Omega_{\\_}1s_{\\_}y$
can be realized by adding microwave fields as discussed above. To realize the
spatial-modulated coupling $\Omega_{\\_}2\cos(k_{\\_}0x)s_{\\_}x$, the
microwave fields may not work, but the Raman process can be used to instead.
We can add two counter-propagating Raman lasers with frequency $\omega_{\\_}3$
to constitute the spatial-modulated field as
$\propto\cos(k_{\\_}0x)e^{-i\omega_{\\_}3t}$. If the frequency is fine tuned
as $\omega_{\\_}3-\omega_{\\_}1$ is matching the energy detuning between
different spin components on $a(c)$ sublattices, the
$\Omega_{\\_}2\cos(k_{\\_}0x)s_{\\_}x$ form of coupling can be realized. The
spatial-modulated coupling on $b(d)$ sublattices can be constituted with
another two counter-propagating Raman lasers with frequency $\omega_{\\_}4$
where $\omega_{\\_}4-\omega_{\\_}1$ is matching the energy detuning on $b(d)$
sublattices. A relative phase is necessary to these two Raman lasers to add an
opposite sign before the spatial-modulated coupling on $b(d)$ sublattices
comparing with $a(c)$ sublattices.

##  VI Discussion

In this work, we propose a coupled-BBH model with different types of coupling
and explore the interplay between higher-order topology and glide reflection
symmetry, where interesting band crossings and degeneracies are predicted. We
have also discussed the corresponding experimental implementation methods of
two schemes in ultracold atomic systems. The band crossing points and nodal
rings discussed in our work can be observed with band mapping
technicsbandmapping . In principle, there is no difficulties to construct
these systems and observe the corresponding properties in other synthetic
quantum systems and metamaterials.

In the future, we will extend our exploration to higher-dimensional system.
The influence of nonsymmorphic symmetries to the three-dimensional higher-
order topological systems and their surface states should be more interesting
and have more potential applications than the present work. Another
interesting research in the future is the generalization of the present model
into many-body systems by considering the interaction between different spin
components. We hope that the interplay between higher-order topology,
nonsymmorphic symmetries and interaction will lead to more interesting quantum
phenomena.

## Acknowledgements

We acknowledge the useful discussion with Xi-Wang Luo. The work is supported
by the National Natural Science Foundation of China (Grant No. 12174138).

##  Appendix A The topological properties of system in Section III when
$\lambda\ll\gamma$

When $\lambda\ll\gamma$, the system with coupling Hamiltonian (14) has no
topological phase transition. With the increase of coupling amplitude, the
system is always HOTI with four zero-energy corner states. As shown in fig.8,
with the increase of coupling amplitude, the band width increase and the
energy of middle two bands along $k_{\\_}x=\pi$ will approach to zero. But
they will never touch the zero energy line. From the energy spectrum (16), the
position of gap minimum can be estimated as

$k_{\\_}y=2\cos^{-1}\frac{\Omega(\lambda+\gamma)}{2(\Omega^{2}+\lambda\gamma)}$
(21)

which approaching to $\pi$ when $\Omega$ is far more larger than $\lambda$ and
$\gamma$. The corresponding minimum gap is about $\Delta
E\approx\frac{(\lambda+\gamma)^{2}\Omega^{3}}{4(\Omega^{2}+\lambda\gamma)^{2}}$
which approaching to zero.

With open boundary along $x$ direction, the band structure is different with
the case we discuss in Section III. There are two separated bands which
corresponding to edge states exist in the middle of the bulk energy bands.
With the increase of coupling amplitude, these two edge bands move approach to
zero energy line, but they will also never touch it.

Figure 8: Upper: The band structure of system with nearest-neighbor coupling
when $\lambda/\gamma=0.1$ for different coupling amplitudes. (a)
$\Omega/\gamma=0.3$, (b) $\Omega/\gamma=0.8$. Lower: The band structure with
open boundary along $x$ direction with corresponding system parameters.

To show the topological properties of the system. We can consider the extreme
case where $\lambda=0$. In this situation, four zero-energy corner states are
all localized on the diagonal lines of the lattice and have the analytical
form. One of them can be written as

$\begin{split}&|\Psi\rangle_{\\_}\mathrm{corner}\propto\\\
&\sum_{\\_}nx^{2(n-1)}\big{(}a^{\dagger}_{\\_}{nn,\uparrow}+ixa^{\dagger}_{\\_}{nn,\downarrow}+xd^{\dagger}_{\\_}{nn,\downarrow}+ix^{2}d^{\dagger}_{\\_}{nn,\uparrow}\big{)}|0\rangle\end{split}$
(22)

where $|0\rangle$ is the vacuum state and

$x=\frac{(1+i)}{2\Omega}\big{(}\lambda-\sqrt{2\Omega^{2}+\lambda^{2}}\big{)}$
(23)

It is easy to verify that state $|\Psi\rangle_{\\_}\mathrm{corner}$ is the
eigenstate with zero energy. One can also easily verify that $|x|<1$ is
satisfied in the whole range of $\Omega$, which means that this state is
localized on the corner of lattice and is one of the zero-energy corner
states.

## References

  * (1) W. A. Benalcazar, B. A. Bernevig, T. L. Hughes, Science 357, 61 (2017).
  * (2) W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Phys. Rev. B, 96, 245115 (2017).
  * (3) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (4) L. Trifunovic and P. W. Brouwer, Phys. Rev. X 9, 011012 (2019).
  * (5) R. Banerjee, S. Mandal, and T. C. H. Liew, Phys. Rev. Lett. 124, 063901 (2020).
  * (6) C. Li, S. Zhang, J. Li, and B. Trauzettel, Phys. Rev. Lett. 127, 026803 (2021).
  * (7) F. K. Kunst, G. van Miert, and E. J. Bergholtz, Phys. Rev. B 97, 241405(R) (2018).
  * (8) Z. Yan, F. Song, and Z. Wang, Phys. Rev. Lett. 121, 096803 (2019).
  * (9) Q. Zeng, Y. Yang, and Y. Xu, Phys. Rev. B 101, 241104(R) (2020).
  * (10) M. Ezawa, Phys. Rev. B 98, 045125 (2018).
  * (11) F. Liu, H. Deng, and K. Wakabayashi, Phys. Rev. Lett. 122, 086804 (2019).
  * (12) D. Călugăru, V. Juričić, and B. Roy, Phys. Rev. B 99, 041301(R) (2019).
  * (13) X. Luo, and C. Zhang, Phys. Rev. Lett. 123, 073601 (2019)
  * (14) M. S.-Garcia, V. Peri, R. Susstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Nature (London) 555, 342 (2018).
  * (15) S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, and R. Thomale, Nat. Phys. 14, 925 (2018).
  * (16) C. W. Peterson, W. A. Benalcazar, T. L. Hughes, and G. A. Bahl, Nature (London) 555, 346 (2018).
  * (17) H. Xue, Y. Yang, F. Gao, Y. Chong, and B. Zhang, Nat. Mater. 18, 108 (2019).
  * (18) X. Ni, M. Weiner, Andrea Alu, and A. B. Khanikaev, Nat. Mater. 18, 113 (2019).
  * (19) X. D. Chen, W. M. Deng, F. L. Shi, F. L. Zhao, M. Chen, and J. W. Dong, Phys. Rev. Lett. 122, 233902 (2019).
  * (20) F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, L. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, and T. Neupert, Nat. Phys. 14, 918 (2018).
  * (21) R. Noguchi, M. Kobayashi, Z. Jiang, K. Kuroda, T. Takahashi, Z. Xu, D. Lee, M. Hirayama, M. Ochi, T. Shirasawa, P. Zhang, C. Lin, C. Bareille, S. Sakuragi, H. Tanaka, S. Kunisada, K. Kurokawa, K. Yaji, A. Harasawa, V. Kandyba, A. Giampietri, A. Barinov, T. K. Kim, C. Cacho, M. Hashimoto, D. Lu, S. Shin, R. Arita, K. Lai, T. Sasagawa, and T. Kondo, Nat. Mater. 20, 473 (2021).
  * (22) F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, and T. Neupert, Sci. Adv. 4, eaat0346 (2018).
  * (23) J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, and P. W. Brouwer, Phys. Rev. Lett. 119, 246401 (2017). Phys. Rev. Lett. 119, 246401 (2017).
  * (24) Z. Song, Z. Fang, and C. Fang, Phys. Rev. Lett. 119, 246402 (2017).
  * (25) M. Geier, L. Trifunovic, M. Hoskam, and P. W. Brouwer, Phys. Rev. B 97, 205135 (2018).
  * (26) E. Khalaf, Phys. Rev. B 97, 205136 (2018).
  * (27) G. van Miert and C. Ortix, Phys. Rev. B 98, 081110(R) (2018).
  * (28) W. A. Benalcazar, T. Li, and T. L. Hughes, Phys. Rev. B 99, 245151 (2019).
  * (29) S. M. Young and C. L. Kane, Phys. Rev. Lett. 115, 126803 (2015).
  * (30) C. Fang and L. Fu, Phys. Rev. B 91, 161105(R) (2015).
  * (31) Y. X. Zhao and A. P. Schnyder, Phys. Rev. B 94, 195109 (2016).
  * (32) H. Chen, X. J. Liu, and X. C. Xie, Phys. Rev. A 93, 053610 (2016).
  * (33) K. Shiozaki, M. Sato, and K. Gomi, Phys. Rev. B 91, 155120 (2015).
  * (34) Q. Z. Wang and C. X. Liu, Phys. Rev. B 93, 020505(R) (2016).
  * (35) S. Li, Y. Liu, S. S. Wang, Z. M. Yu, S. Guan, X. L. Sheng, Y. Yao, and S. A. Yang, Phys. Rev. B, 97, 045131 (2018).
  * (36) T. Tan, C. Li, and W. Yao, Phys. Rev. B, 204, 245419 (2021).
  * (37) Z. Wang, B. J. Wieder, J. Li, B. Yan, and B. A. Bernevig, Phys. Rev. Lett. 123, 186401 (2019).
  * (38) X. Zhang, Z. K. Lin, H. X. Wang, Z. Xiong, Y. Tian, M. H. Lu, Y. F. Chen and J. H. Jiang, Nat. Commun. 11, 65 (2020).
  * (39) X. Zhou, C. H. Hsu, C. Y. Huang, M. Iraola, J. L. Mañes, M. G. Vergniory, H. Lin, and N. Kioussis, npj Computational Materials 7, 202 (2021).
  * (40) W. Yao, S. A. Yang, and Q. Niu, LABEL:https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.102.096801Phys. Rev. Lett. 102, 096801 (2009).
  * (41) W. P. Su, J. R. Schrieffer, and A. J. Heeger, Phys. Rev. Lett. 42, 1698 (1979).
  * (42) M. Atala, M. Aidelsburger, J. T. Barreiro, D. Abanin, T. Kitagawa, E. Demler, and I. Bloch, Nat. Phys. 9, 795 (2013).
  * (43) M. Aidelsburger, M. Atala, S. Nascimbène, S. Trotzky, Y.-A. Chen, and I. Bloch, Phys. Rev. Lett. 107, 255301 (2011).
  * (44) O. Mandel, M. Greiner, A. Widera, T. Rom, T. W. Hänsch, and I. Bloch, Phys. Rev. Lett. 91, 010407 (2003).
  * (45) P. Soltan-Panahi, J. Struck, P. Hauke, A. Bick, W. Plenkers, G. Meineke, C. Becker, P. Windpassinger, M. Lewenstein, and K. Sengstock, Nat. Phys. 7, 434 (2011).
  * (46) B. Yang, H.-N. Dai, H. Sun, A. Reingruber, Z.-S. Yuan, and J.-W. Pan, Phys. Rev. A, 96, 011602 (2017).
  * (47) L. Huang, Z. Meng, P. Wang, P. Peng, S. L. Zhang, L. Chen, D. Li, Q. Zhou and J. Zhang Nat. Phys. 12, 540 (2016).
  * (48) 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).
  * (49) M. C. Liang, Y. D. Wei, L. Zhang, X. J. Wang, H. Zhang, W. W. Wang, W. Qi, X. J. Liu, and X. Zhang, arXiv:2109.08885.
  * (50) L. Tarruell, D. Greif, T. Uehlinger, G. Jotzu, and T. Esslinger, Nature, 483, 302 (2012).