Takagi Topological Insulator on the Honeycomb Lattice

Introduction. The symmetry-protected topological phases, such as topological (crystalline) insulators (TIs) and superconductors (TSCs), have been one of the most active fields of physics during the last fifteen years Volovik (2003); Hasan and Kane (2010); Qi and Zhang (2011); Fu (2011); Chiu et al. (2016); Kruthoff et al. (2017); Benalcazar et al. (2017); Liu et al. (2019); Xie et al. (2021). Based on the topological $K$ theory, the topological band theory has been established to classify and characterize various topological states Atiyah (1966); Kitaev (2010); Schnyder et al. (2008). Symmetry plays an fundamental role in the classification of topological phases. Considering three discrete symmetries, namely time reversal $\mathcal{T}$, charge conjugation $\mathcal{C}$ and chiral symmetry $\mathcal{S}$, physical systems can be classified into ten symmetry classes, termed Altland-Zirnbauer (AZ) classes Altland and Zirnbauer (1997); Kitaev (2010); Hořava (2005); Atiyah (1966); Schnyder et al. (2008); Zhao and Wang (2014), among which the eight ones with at least $\mathcal{T}$ or $\mathcal{C}$ are called real AZ classes. The topological classifications in the framework of the eight real AZ classes correspond to the real $K$ theory. Using the real $K$ theory, gapped systems including topological insulators and topological superconductors were first classified Schnyder et al. (2008); Kitaev (2010); Ryu et al. (2010), and then gapless systems were classified as well Matsuura et al. (2013); Zhao and Wang (2013); Chiu and Schnyder (2014); Shiozaki and Sato (2014). All the classification tables exhibit an elegant eightfold periodicity along the dimensions for the eight real AZ classes.

After internal symmetries like $\mathcal{T}$ and $\mathcal{C}$, more and more spatial symmetries were involved to enrich symmetry-protected topological matter. It was noticed that combined symmetries $\mathcal{P}\mathcal{T}$ and $\mathcal{C}\mathcal{P}$ correspond to the orthogonal $K$ theory with $\mathcal{P}$ the spatial inversion, since they leave every $k$ point fixed in the reciprocal space. Hence, the topological classification table was worked out Zhao et al. (2016). A remarkable feature is that groups $\mathbb{Z}$, $\mathbb{Z}_{2}$ and $0$ in the table appear in the reversed order in dimensionality, compared with previous tables for the real AZ classes. $\mathcal{P}\mathcal{T}$ and $\mathcal{C}\mathcal{P}$ are fundamental in nature, and therefore the classification table has been applied to explore topological phases in various physical systems, such as quantum materials Zhao and Lu (2017); Kruthoff et al. (2017); Ahn et al. (2019), topological superconductors Timm et al. (2017); Yu et al. (2021); Tomonaga et al. (2021); Lapp et al. (2020), and photonic/phononic crystals and electric-circuit arrays  Zhang et al. (2013); Yang et al. (2015); Imhof et al. (2018); Ozawa et al. (2019); Ma et al. (2019); Serra-Garcia et al. (2018); Yu et al. (2020); Peterson et al. (2018); Lapp et al. (2020), and can generate unique topological structures with many novel consequences, such as non-Abelian topological charges, cross-order boundary transitions, and nodal-loop linking structures  Zhao and Lu (2017); Yu et al. (2015); Ahn et al. (2019); Sheng et al. (2019); Wu et al. (2019); Wang et al. (2019); Li et al. (2020).

Remarkably, from the classification table, the symmetry class with $(\mathcal{P}\mathcal{T})^{2}=1$ corresponds to the $\mathbb{Z}_{2}$ classification for $d=1$ and $d=2$. As revealed in Ref.Zhao and Lu (2017), $(\mathcal{P}\mathcal{T})^{2}=1$ leads to real band structures in contrast to conventional complex band structures. Then, the $\mathbb{Z}_{2}$ topological invariant $w_{1}$ for $d=1$ can be formulated as the quantized Berry phase in units of $\pi$ modulo $2\pi$. The case of $d=2$ is much fascinating. The topological invariant is the Euler number, a real version of the Chern number, for two valence bands. The Euler number is valued in $\mathbb{Z}$, but only its parity is stable if more trivial valence bands are added into consideration. The parity, namely the Euler number modulo $2$, is just the Stiefel-Whitney number $w_{2}$ in two dimensions, which determines whether the real vector bundle can be lifted into a spinor bundle.

The topological invariant $w_{2}$ gives rise to novel topological phases with extraordinary properties. In $3$D, it characterizes a real Dirac semimetal, which can be transformed into a nodal ring with symmetry-preserving perturbations. Then, the nodal ring is characterized by two topological charges $(w_{1},w_{2})$. In $2$D, it describes a topological insulator. The common topological wisdom is that the bulk topological invariant determines a unique form of the boundary modes, namely the well-known one-to-one bulk-boundary correspondence. However, a remarkably discovery in Ref.Wang et al. (2020) is that $w_{2}$ corresponds to multiple forms of boundary modes, extending the one-to-one correspondence to one-to-many. The $2$D topological insulator can host various second-order phases with odd $\mathcal{P}\mathcal{T}$-related pairs of corner zero-modes, which are mediated by first-order phases with helical edge states. Similarly, the $3$D semimetal can host second-order hinge Fermi arcs and first-order surface Dirac states as well. Recently, graphynes have been proposed as the material candidates which can realize both the $2$D topological insuslator and the $3$D topological semimetal Sheng et al. (2019); Chen et al. (2021, 2022).

As aforementioned, the second-order phases of the $2$D topological insulator feature odd $\mathcal{P}\mathcal{T}$-related pairs of corner zero modes. It is interesting to look for its $3$D analog, which has been presented in Ref.Dai et al. (2021). Referring to the topological classification table for $\mathcal{P}\mathcal{T}$ and $\mathcal{C}\mathcal{P}$ symmetries, we notice that although the classification for $(\mathcal{P}\mathcal{T})^{2}=1$ is trivial, with an additional chiral symmetry $\mathcal{S}$ with $\{\mathcal{P}\mathcal{T},\mathcal{S}\}=0$ the classification is preserved as $\mathbb{Z}_{2}$ in $2$D and, more importantly, becomes nontrivial as $\mathbb{Z}_{2}$ in $3$D. It is found that the corresponding topological invariants can be formulated in terms of Takagi’s factorization. The topological invariant in $2$D is equivalent to $w_{2}$, while that in $3$D is a new topological invariant. Either in $2$D or in $3$D the bulk topological invariant can be manifested as odd $\mathcal{P}\mathcal{T}$-related pairs of corner zero-modes. Now, with the chiral symmetry, the two zero-modes in each pair are eigenstates with opposite eigenvalues of the chiral symmetry.

In this article, we review the elements of $2$D $\mathcal{P}\mathcal{T}$-protected topological insulators with or without chiral symmetry. The essential physics is demonstrated by the Honeycomb-lattice model, with only the nearest-neighbor hopping amplitudes. We show that under certain dimerization patterns the model is a topological insulator with nontrivial Stiefel-Whitney number or the Takagi topological invariant, and therefore presents all the nontrivial topological phenomena. Particularly, under various $\mathcal{P}\mathcal{T}$-invariant geometries, there are always odd $\mathcal{P}\mathcal{T}$-related pairs of corner zero-modes for the second-order topological phase. Before diving into the details, it is noteworthy that the dimerized honeycomb model can be regarded as an abstraction from the graphynes Sheng et al. (2019); Chen et al. (2021, 2022).

The honeycomb-lattice model.

Let us start with presenting the honeycomb-lattice model, the lattice structure is shown in Fig. 1(a). The Hamiltonian in momentum space is given by

where $\chi^{(i)}_{\bm{k}}=t_{i}^{\prime}e^{-i\bm{k}\cdot\bm{a}_{i}}$ with $i=1,2,3$. Here, $\bm{a}_{i}$ are the bond vectors connecting the centers of nearest-neighbor unit cells, as indicated in Fig.1(a) with $\sum_{i}\bm{a}_{i}=0$. The Hamiltonian has inversion symmetry with $\hat{\mathcal{P}}=\sigma_{1}\otimes I_{3}\hat{I}$, spinless time-reversal symmetry with $\hat{\mathcal{T}}=\hat{\mathcal{K}}\hat{I}$, and therefore spacetime-inversion symmetry with $\hat{\mathcal{P}}\hat{\mathcal{T}}=\sigma_{1}\otimes I_{3}\hat{\mathcal{K}}$, where $\hat{\mathcal{K}}$ is the complex conjugation and $\hat{I}$ is the inversion of momenta. Note that $\sigma$’s are the Pauli matrices acting on the sublattice space [see Fig.1(a)], and $I_{3}$ is the $3\times 3$ identity matrix. Each inversion center is taken as the center of a hexagon in real space. The sublattice symmetry operator is $\hat{\mathcal{S}}=I_{3}\otimes\sigma_{3}$. Since the inversion exchanges sublattices, both $\mathcal{P}$ and $\mathcal{P}\mathcal{T}$ anti-commute with $\mathcal{S}$, namely, $\{\hat{\mathcal{P}},\hat{\mathcal{S}}\}=\{\hat{\mathcal{P}}\hat{\mathcal{T}},\hat{\mathcal{S}}\}=0$.

To obtain the nontrivial topological phases, we calculate the determinant of the Hamiltonian (1) at $\Gamma$ point Gam in the Brillouin zone as

Since the bulk topological criticality generally corresponds gap-closing point, we can obtain the topological phase-transition points by letting $\mathrm{det}[\mathcal{H}(\Gamma)]=0$, which gives

Interestingly, if (3) holds, the system is generally reduced to a topologically equivalent graphene model with two Dirac points in the first Brillouin zone Haldane (1988). When $t_{1}^{2}t^{\prime}_{1}+t_{2}^{2}t^{\prime}_{2}+t_{3}^{2}t^{\prime}_{3}<2t_{1}t_{2}t_{3}+t^{\prime}_{1}t^{\prime}_{2}t^{\prime}_{3}$, the system steps into a topological phase, while conversely the system becomes a trivial phase, which can be checked by computing Stiefel-Whitney number or Takagi’s factorization.

Topological invariants The topology can be determined by various formulas of the topological invariant. We now briefly review them. First, as given in Ref.Zhao and Lu (2017), the topological invariant can be determined by the Wilson loop

(with $P$ indicating the path order) along large circles parametrized by $k_{x}$. $C_{k_{y}}$ is the contour at a fixed $k_{y}$ and $\mathcal{A}(k_{x},k_{y})$ is the non-Abelian Berry connection for the valence bands. The topological information is encoded in the phase factors $\theta(k_{y})\in\left(-\pi,\pi\right]$ of the $N$ eigenvalues $\lambda_{m}(k_{y})$ of $\mathcal{W}(k_{y})$ for valence bands:

Different from the conventional TIs and Chern insulators, the Wilson loop spectral flow for real phases are mirror symmetric with respect to the $\theta=0$ axis [see Fig. 1(b)]. This is because $\mathcal{W}(k_{y})$ is equivalent to a mapping from $k_{y}\in S^{1}$ to $O(N)$ up to a unitary transformation Zhao and Lu (2017). The topological information can be pictorially derived from counting how many times $\zeta$ the trajectories cross $\theta=\pi$ as

For honeycomb lattice with $t_{1}^{2}t^{\prime}_{1}+t_{2}^{2}t^{\prime}_{2}+t_{3}^{2}t^{\prime}_{3}<2t_{1}t_{2}t_{3}+t^{\prime}_{1}t^{\prime}_{2}t^{\prime}_{3}$, a single crossing exsit as shown in Fig.1(b), namely, $w_{2}=1$, which indicates the model is in a topological nontrivial phase.

As aforementioned, our system is protected by spacetime inversion symmetry $\mathcal{P}\mathcal{T}$ and sublattice (chiral) symmetry $\mathcal{S}$. These symmetries constraint the classifying space of $\mathcal{H}(\bm{k})$ to be symmeric unitary matrices. Thus the $\mathbb{Z}_{2}$ invariant from the Takagi’s factorization can be defined Dai et al. (2021), which leads to an alternative formulation for $w_{2}$. We now prove the equivalence of the two formulas. For technical simplicity, we assume the momentum space as a sphere $S^{2}$, which is sufficient to present the essential ideas.

In general, $\mathcal{S}$ requires the Hamiltonian $\mathcal{H}(\bm{k})$ to be block anti-diagonal and $\mathcal{P}\mathcal{T}$ requires the upper-right block to be symmetric. Thus the flattened Hamiltonian $\tilde{\mathcal{H}}(\bm{k})$ is given by

where $\mathcal{Q}(\bm{k})=\mathcal{U}(\bm{k})\mathcal{U}^{T}(\bm{k})$ is a unitary symmetric matrix for each $\bm{k}$ and $M$ denotes the number of valence (conduction) bands. $\mathcal{U}(\bm{k})\in U(M)$ is the Takagi factor. The classifying space for this symmetric class is ${US}(M)=U(M)/O(M)$ Dai et al. (2021). Here, $\pi_{2}[{US}(M)]=\mathbb{Z}_{2}$ corresponds to the topological invariant of our system. Consider a $2$D sphere $S^{2}$, which is divided into north and south hemispheres $D^{2}_{N,S}$, overlapping along the equator $S^{1}$. The Takagi factors $\mathcal{U}_{N/S}$ over $D_{N/S}^{2}$, respectively, can be transformed to each other by a gauge transformation $\mathcal{O}_{S^{1}}$ over the equator $S^{1}$, as shown in Fig. 2. $\mathcal{O}_{S^{1}}$ is given by

$\pi_{1}[{O}(M)]=\mathbb{Z}_{2}$ for $M>2$ leads to obstructions for a global Takagi’s factorization over $S^{2}$.

The conduction and valence wavefunctions of $\tilde{\mathcal{H}}(\bm{k})$ can be given by

where $n\in\{1,2,\cdots,M\}$. $\varphi_{n}=(0\leavevmode\nobreak\ 0\leavevmode\nobreak\ \cdots\leavevmode\nobreak\ 0\leavevmode\nobreak\ 1\leavevmode\nobreak\ 0\leavevmode\nobreak\ 0\leavevmode\nobreak\ \cdots\leavevmode\nobreak\ 0)^{T}$ is a unit vector with “$1$” locating at the $n$-th position.

Performing a unitary transformation $\mathcal{U}_{R}=e^{-i{\pi}/{4}}e^{i{\pi\sigma_{1}}/{4}}$ on this system, the Hamiltonian and valence wavefunctions both become real. Meanwhile, $\mathcal{P}\mathcal{T}$ and $\mathcal{S}$ are transformed to $\hat{\mathcal{K}}$ and $\sigma_{2}$, respectively. Over the intersection $S^{1}$, transition function $t_{S^{1}}$ of real valence wavefunctions can be given by

Thus, we know the transition function $t_{S^{1}}$ of real valence wavefunctions is equal to the gauge transformation $\mathcal{O}_{S^{1}}$. As noted in Ref.Zhao and Lu (2017), $w_{2}$ is just the parity of the winding number of the transition function for valence bands. Thus, we see the equivalence of two $2$D topological invariants.

Physical consequence. According to analytical and numerical methods, we reveal that three pairs of hopping parameters $t_{i}$ and $t^{\prime}_{i}$ (with $i=1,2,3$) jointly determine the configuration of topological boundary modes. To facilitate understanding the relation between distinct boundary modes and parameters, we define a boundary effective mass term $m_{i}$ for each edge:

where the subscript $i$ denotes the hopping along the primitive vector $\mathbf{a}_{i}$ direction ($\mathbf{a}_{3}=-\mathbf{a}_{1}-\mathbf{a}_{2}$). The above Eq.(10) can be derived from the boundary effective Hamiltonian Edg . Hence, if $m_{i}=0$, the corresponding edges are gapless, which is also the boundary critical condition to separate two second-order topological phases.

To demonstate the boundary modes, we consider a rhombic-shaped $2$D sample with armchair termination, i.e., by opening boundary along $\mathbf{a}_{1}$ and $\mathbf{a}_{2}$ direction, as shown in Fig. 3. If $m_{1}=0$ and $m_{2}\neq 0$, the helical edge modes along periodic $\mathbf{a}_{1}$ can be obtained, as shown in Fig. 3(b). However, once $m_{1,2}\neq 0$, the helical edge modes will be gapped and the localized corner modes will emerge. More specifically, for the case with $\mathrm{sgn}(m_{1})=\mathrm{sgn}(m_{2})$ ($\mathrm{sgn}(m_{1})=-\mathrm{sgn}(m_{2})$), the corner modes will locate at $120^{\circ}$ ($60^{\circ}$) corners, as shown in Fig.3(a) and (c) respectively. The $\mathcal{P}\mathcal{T}$ symmetry requires that the corner zero-modes always come in pairs and the chiral symmetry sets the midgap modes exactly at zero energy Fin .

To keep the completeness of honeycomb unit cell in a rhomboid sample, one only has three kinds of armchair edges, namely the edge parallel to $\mathbf{a}_{i}$ direction with $i=1,2,3$. If the edge connected by the same corner has the same mass term $m_{i}$ sign, the corner zero-modes will be localized at the obtuse angle of the rhomboid, otherwise at acute corners. We shall theoretically explain these numerical results in the next section. It is emphasized that in the whole process of the edge-phase transitions, the bulk gap is always open and the symmetries are preserved, therefore, the bulk invariant $\nu$ is unchanged. Thus the conventional bulk-boundary correspondence is not appliable for TTI, namely, the bulk invariant can not uniquely determine the boundary modes, but dictates an edge criticality, as the concept mentioned in previous work Wang et al. (2020).

As promised in introduction, we now proceed to tune the boundary modes with fixed parameters. In the rhombic case, all samples terminate with armchair edges and exhibit parameter-depended boundary modes. As long as $\mathcal{P}\mathcal{T}$ and $\mathcal{S}$ are not violated, the finite samples can be cut with not only rhomb as shown in Fig. 3, but also hexagon(see Fig.5 ). Beside armchair edges, zigzag edges can serve as termination too. Creatively, with fixed parameters but different boundary selections, one can also find various distinguishable boundary modes. For instance, helical edge states emerge on the zigzag edges in a rectangle sample as shown in Fig. 4(b), with the same hopping parameters as Fig. 4(a). This result further proves that the bulk topological invariant can not uniquely determine the topological boundary modes. We also study lots of other patterns with the same parameters, and abundant topological boundary modes consisting of corner zero modes and gapless edge modes can be obtained (see Appendix. D). They are all boundary-selection-depended. Hence, we propose that one can obtain needed topological boundary modes by choosing particular boundary geometry, without tuning parameters, which is usually difficult to perform in real systems.

Novelly, in both situations discussed above, for second-order topological phases, the number of the zero-energy corners must be odd pairs. For example, for a hexagonal sample, we can only find one or three pairs of zero-mode corners, as shown in Fig. 5(a) and (b). Similarly, the Octagonal sample also has one or three pairs of zero-mode corners as shown in Fig. 5(c) and (d).

We find that the peculiarity of odd-pair-zero-modes is universal and it can be generalized to a higher-dimensional situation, such as 3D Dai et al. (2021). The 3D model is constructed by stacking the 2D honeycomb TTI discussed above in a staggered manner to preserve the sublattice symmetry $\mathcal{S}$. The details of the construction of the model can be found in Appendix. C . The inversion center is chosen as the center of hexagonal in one layer. Thus the anti-commuting relation of $\mathcal{P}\mathcal{T}$ and $\mathcal{S}$ are preserved. We cut a finite hexagonal prisms sample that keeps the symmetries. One can find only odd pairs (one or three) of corner states related by $\mathcal{P}\mathcal{T}$ appear. Note that the corner zero-modes can be driven to other corners by tunning the hopping parameters like in a 2D situation.

Analytic method. We first proceed to solve the boundary criticality along the periodic $\mathbf{a}_{1}$ direction and openning boundary with $\mathbf{a}_{2}$. Replace $e^{-i\bm{k}\cdot\mathbf{a_{2}}}$ by $S$ in Hamiltonian (1), with $S$ ladder operator and $S|i\rangle=|i+1\rangle,S^{\dagger}|i\rangle=|i-1\rangle$. Then $e^{-i\bm{k}\cdot\mathbf{a}_{3}}$ can be represented by $e^{i\bm{k}\cdot\mathbf{a}_{1}}S^{\dagger}$ since $\mathbf{a}_{3}=-(\mathbf{a}_{1}+\mathbf{a}_{2})$. After a series of tedious derivation (see Appendix. A), we obtian the effective Hamiltonian of the bottom boundary:

Following the same argument with aforementioned bulk criticality, we can obtain the boundary criticality by letting $\mathrm{det}[\mathcal{H}_{B}(\bm{k}\cdot\mathbf{a}_{1})]=0$, which leads to

The above Eq. (12) holds only at $\bm{k}\cdot\mathbf{a}_{1}=0$. Thus, when the system is in a topological nontrivial case, Eq. (12) related edge criticality separates two different second-order topological phases with corner zero modes. Likely, we can obtain similar results for periodic $\mathbf{a}_{2}$ and $\mathbf{a}_{3}$ directions. For convenience, we can define boundary effective mass by the left of Eq. (12) for edges parallel to $\mathbf{a}_{1}$. Or, generally Eq. (10) for edges parallel to $\mathbf{a}_{i}$. Different from the armchair edges, the zigzag terminations has an additional boundary criticality, namely, the effective masses can be defined by

The derived details can be found in the Appendix. B.

With the effective mass orderly distributing on each edge Edg , the existence of corner zero-modes is reduced to a Jackiw-Rebbi problem Jackiw and Rebbi (1976). The corners with opposite effective masses on both sides can have zero-modes.

Discussion. In this article, we present a simple 2D realizable honeycomb-lattice model to demonstrate the essential physics of the Takagi topological insulator. It is found that with unchanged topological invariant, one can tune topological boundary modes by not only parameters, but also boundary selections. It goes beyond the common wisdom about bulk-boundary correspondence, and gives rise to much richer boundary physics.

Our model with novel physics is closely related to real systems. It is easier to realize our model by photonic/phononic crystals, electric-circuit arrays and mechanics systems, since only have nearest-neighbor hopping amplitudes are included into the model. Several special cases of our model have been recently realized in photonic/phononic crystals Yang et al. (2020); Noh et al. (2018), where hopefully the general form of our model can be further experimentally examined.