Higher-order topological insulator with glide reflection symmetry

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.

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.

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}$.

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.

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)

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.

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.

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.

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.

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$.

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.

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.

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.

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