Two-dimensional higher-order topology in monolayer graphdiyne

Abstract

Based on first-principles calculations and tight-binding model analysis, we propose monolayer graphdiyne as a candidate material for a two-dimensional higher-order topological insulator protected by inversion symmetry. Despite the absence of chiral symmetry, the higher-order topology of monolayer graphdiyne is manifested in the filling anomaly and charge accumulation at two corners. Although its low energy band structure can be properly described by the tight-binding Hamiltonian constructed by using only the p_z orbital of each atom, the corresponding bulk band topology is trivial. The nontrivial bulk topology can be correctly captured only when the contribution from the core levels derived from p_x,y and s orbitals are included, which is further confirmed by the Wilson loop calculations. We also show that the higher-order band topology of a monolayer graphdyine gives rise to the nontrivial band topology of the corresponding three-dimensional material, ABC-stacked graphdiyne, which hosts monopole nodal lines and hinge states.

Introduction

Bulk-boundary correspondence is a fundamental property of topological phases. In conventional topological insulators, the gapped bulk states in d-dimensions support metallic states in (d − 1)-dimensional surfaces. Such a conventional bulk-boundary correspondence, however, is violated in a class of topological crystalline insulators, so-called higher-order topological insulators (HOTIs). The gapless excitations of a HOTI in d-dimensions are localized, instead, in a subspace with a dimension lower than (d − 1), such as corners or hinges, when the global shape of the material preserves the crystalline symmetry relevant to the nontrivial band topology.^{1,2,3,4,5,6,7,8,9,10,11} Recently, rhombohedral bismuth has been identified as the first example of 3D HOTIs hosting helical hinge states.¹ Also there are other candidate materials of 3D HOTIs including SnTe with strain along the (100) direction,² transition metal dichalcogenides MoTe₂ and WTe₂, hosting helical hinge states,³ and also Bi_2−xSm_xSe₃ with chiral hinge states.⁴

In 2D, one can also consider second-order topological insulators (SOTIs) hosting localized corner states.^{7,8,9,10,11,12} Contrary to 3D HOTIs, however, there are only few theoretical proposals for candidate materials realizing 2D SOTIs.^{9,10,11,13,14} For instance, phosphorene is one candidate for 2D HOTIs proposed based on simple tight-binding model calculations.¹⁵ However, in this system, the corner states originate from the charge polarization,¹⁵ which is a 1D topological invariant. Hence, in a strict sense, phosphorene is not a genuine 2D HOTI.¹³ Another candidate is twisted bilayer graphene, which can also exhibit higher-order band topology at small twist angles. However, to realize the corner states, the twist angle, as well as the chemical potential should be controlled simultaneously.

In general, among the higher-order TIs, only d-dimensional kth-order TIs with 1 < k < d have stable band topology with Wannier obstruction while d-dimensional dth-order TIs with zero-dimensional corner charges are obstructed atomic insulators. Thus, in 3D, a SOTI has Wannier obstruction while a third-order TI (TOTI) is atomic. On the other hand, in 2D, HOTIs are always obstructed atomic insulators. In d-dimensional dth-order TIs, additional chiral or particle-hole symmetry is required to protect corner states which appear as in-gap states separated from bulk states.^8–11 Since chiral or particle-hole symmetry usually does not exist in insulators, realizing 2D SOTIs in materials is generally considered difficult.

However, even if the chiral symmetry of 2D SOTIs is broken, the material still inherits the nontrivial band topology of 2D SOTIs. The resulting state is called a 2D TOTI.^16,17,18 The key idea is that even if the in-gap states are merged into the bulk states, the nontrivial band topology manifests in the filling anomaly.^5,11,12 Namely, the half-filling condition can never be satisfied as long as the symmetry associated with the second-order band topology is preserved. The half-filling condition can be satisfied only when extra electrons or holes are added, which naturally leads to charge accumulation or depletion at corners. Therefore, the hallmark of 2D HOTIs lacking chiral symmetry (or simply 2D TOTIs), is the existence of filling anomaly at half-filling and the accumulation or depletion of a half-integral charge at corners after adding extra electrons or holes.^16,17,18

In this work, we theoretically propose monolayer graphdiyne (MGD) as the first realistic candidate material for 2D HOTIs. Namely, MGD is a 2D SOTI when chiral symmetry exists, but in real systems, it appears as a 2D TOTI due to the lack of chiral symmetry. Based on first-principles calculations and tight-binding model analysis, we show that a MGD is a HOTI characterized by a 2D topological invariant w₂, that is quantized when the system is invariant under inversion P symmetry. Although corner states are not guaranteed to be in the gap due to the lack of chiral symmetry in MGD, the nontrivial higher-order band topology is manifested in filling anomaly and the corner charge accumulation. Also, it is found that although the low energy excitation can be well-described by using the state derived from the p_z orbitals of each atom, the relevant band topology of such a simplified model is trivial. To describe the nontrivial band topology, it is essential to include the contribution from the core levels derived from p_x,y and s orbitals. Finally, it is shown that the nontrivial w₂ of MGD gives rise to the topological band structure of the corresponding three-dimensional material, ABC stacked graphdiyne, hosting nodal lines with Z₂ monopole charges.^19,20 Due to the nonzero w₂ in a range of 2D momentum subspaces with fixed vertical momentum k_z, the ABC stacked graphdiyne host hinge states for the corresponding range of k_z.

Results

Topological invariant of monolayer graphdiyne

MGD is an experimentally realized planar carbon system composed of sp²-sp carbon network of benzene rings connected by ethynyl chains.^21,22 Figure 1a describes the unit cell of MGD composed of 18 carbon atoms. Since spin-orbit coupling (SOC) is negligible, MGD can be regarded as a spinless fermion system. The lattice has inversion P, time-reversal T, a two-fold rotation about the x-axis C_2x, a six-fold rotation about the z-axis C_6z, and a mirror M_z: z → −z symmetries. Also, due to the bipartite lattice structure, MGD has chiral (or sublattice) symmetry when only the nearest neighbor hopping between different sublattices is considered, similar to the case of graphene. Fig. 1b shows the band structure along high symmetry directions obtained by first-principles calculations. Here, the blue color indicates the bands derived from p_z orbitals while the red color denotes the bands derived from s, p_x, p_y orbitals. Since p_z orbitals are odd whereas s, p_x, p_y orbitals are even under M_z symmetry, the energy spectrum from p_z orbitals is not hybridized with that from s, p_x, p_y orbitals. One can see that the low-energy band structure has approximate chiral symmetry and can be effectively described by using only p_z orbitals. However, to capture the higher-order band topology, the core electronic states derived from s, p_x, p_y orbitals should be included as shown below.

Fig. 1: Lattice structure, band structure, and the Wilson loop spectrum.

a A schematic figure describing the unit cell and the relevant atomic orbitals of a monolayer graphdiyne (MGD) composed of 18 carbon atoms. Two sublattices are marked with yellow and blue colors, respectively. Since MGD has a bipartite lattice structure, chiral symmetry exists when only the nearest neighbor hoppings are considered. b The band structure of MGD obtained by first-principles calculations, which shows approximate chiral symmetry near the Fermi level. Since p_z orbitals are odd while s, p_x, p_y orbitals are even under the mirror M_z:z → −z operation, the energy spectrum from p_z orbitals (blue) is not hybridized with that from s, p_x, p_y orbitals (red). c The Wilson loop spectrum computed including s, p_x, p_y, p_z orbitals. The spectrum exhibits a crossing point at k_y = 0, θ = π, which indicates w₂ = 1.

In general, a 2D P-symmetric spinless fermion system carries three Z₂ topological invariants, w_1x, w_1y, and w₂.²⁰ In terms of the sewing matrix G_mn(k) = 〈u_m-k∣P∣u_nk〉, w_1a (a = x, y) is given by its 1D winding number as

$${w}_{1a}=-\frac{i}{2\pi}{\oint_{C}}{\nabla }_{{{\bf{k}}}_{a}}{\log}\,{\mathrm{det}}\,G({\bf{k}})\ {\rm{mod}}\, 2,$$

(1)

which is equivalent to the quantized Berry phase Φ_a, in a way that Φ_a = πw_1a (mod 2π).^23,24 w_1a can be determined by using relation ${(-1)}^{{w}_{1a}}={\prod }_{i=1}^{2}\xi ({\Gamma }_{ia})$ where ξ(Γ_ia) (i = 1, 2) is the product of parity eigenvalues of occupied states at the time-reversal invariant momentum (TRIM) Γ_ia, and Γ_1a, Γ_2a are two TRIMs along the reciprocal lattice vector G_a. Since w_1a is equivalent to the quantized charge polarization along the G_a direction, it can be considered as a 1D topological invariant.

On the other hand, w₂ is a 2D topological invariant given by the 2D winding number of the sewing matrix G(k) (modulo 2),⁶ and measures the higher-order band topology of P-symmetric spinless fermion systems.^5,6,13 w₂ can be determined by

$${(-1)}^{{w}_{2}}=\prod _{i=1}^{4}{(-1)}^{[{N}_{\text{occ}}^{-}({\Gamma }_{i})/2]},$$

(2)

where ${N}_{\,\text{occ}\,}^{-}({\Gamma }_{i})$ is the number of occupied bands with odd parity at the TRIM Γ_i and the square bracket [α] indicates the greatest integer value of α^{25,26,27,28,29,30}.

Using the band structure from first-principles calculations and the corresponding parity eigenvalues at TRIMs, one can easily show that MGD is a higher-order topological insulator with w_1x = w_1y = 0 and w₂= 1. One interesting feature of the MGD band structure is that we get w₂ = 0 when only the bands derived from p_z orbitals are considered. Namely, although the low-energy band structure itself can be well-described by using only p_z orbitals, the higher-order topology can be correctly captured only when the core electronic states derived from s, p_x, p_y orbitals are included. This result can be confirmed by tight-binding analysis, as well as first-principles calculations. Since the bands derived from p_z orbitals do not hybridize with those derived from other orbitals because of M_z symmetry, one can compute the relevant topological invariants separately.

The orbital dependence of w₂ can also be confirmed by Wilson loop calculations.^13,20,31,32 The Wilson loop is defined as a path ordered product of the exponential of Berry connections,^33,34

$${W}_{({k}_{1}+2\pi ,{k}_{2})\leftarrow ({k}_{1},{k}_{2})}=\mathop{{\mathrm{lim}}}\limits_{N\to \infty }{F}_{N-1}{F}_{N-2}\cdots {F}_{1}{F}_{0},$$

(3)

where ${[{F}_{i}]}_{mn}=\langle {u}_{m}(\frac{2\pi (i\,+\,1)}{N},{k}_{2})| {u}_{n}(\frac{2\pi i}{N},{k}_{2})\rangle$, and its spectrum is gauge-invariant.³³ Computed along the k₁ direction parallel to the reciprocal lattice vector G₁ from (k₁, k₂), the set of Wilson loop eigenvalues $\{{e}^{i{\nu }_{i}({k}_{2})}\}$ indicates the position of Wannier centers at given k₂, and the corresponding total charge polarization is given by ${p}_{1}({k}_{2})=\frac{1}{2\pi }{\sum }_{i=1}^{{N}_{{\mathrm{occ}}}}{\nu }_{i}({k}_{2})$. The Wilson loop spectrum of MGD shows that the charge polarization is zero in both k₁ and k₂ directions, which is consistent with the parity eigenvalue analysis (Fig. 1c). Also, the single linear-crossing at (k_y, ν) = (0, π) indicates w₂ = 1.²⁰ Projecting the Wilson loop operator into the p_z (s, p_x, and p_y) orbital basis, we observe w₂ = 0 (w₂ = 1), which confirms that the higher-order band topology originates from the core electronic levels. The orbital dependence of w₂ can be further confirmed by computing the Wannier center of each orbital explicitly and also by computing the nested Wilson loop that measures the higher-order band topology.^5,7 In the presence of PT symmetry, the Wilson loop spectrum can be divided into two sub-bands that are separated by a gap, each centered at ν = 0 or ν = π, respectively. When the Wilson loop operator at given k₂ is diagonalized as

$${W}_{({k}_{1}+2\pi ,{k}_{2})\leftarrow ({k}_{1},{k}_{2})}={\sum_{i}}\left|{\nu }_{i}({k}_{2})\right\rangle {e}^{i{\nu }_{i}({k}_{2})}\left\langle {\nu }_{i}({k}_{2})\right|,$$

(4)

then, using the set of eigenvectors corresponding to the sub-band centered at ν = π, the nested Wilson loop operator $\tilde{W}$ along the G₂ vector is defined as

$${\tilde{W}}_{({k}_{2}+2\pi )\leftarrow ({k}_{2})}=\mathop{{\mathrm{lim}}}\limits_{N\to \infty }{\tilde{F}}_{N-1}{\tilde{F}}_{N-2}\cdots {\tilde{F}}_{0},$$

(5)

where ${[{\tilde{F}}_{i}]}_{mn}=\langle {\nu }_{m}(\frac{2\pi (i\,+\,1)}{N})| {\nu }_{n}(\frac{2\pi i}{N})\rangle$. Since the determinant of the nested Wilson loop operator ${\mathrm{det}}(\tilde{W})$ is identical to ${(-1)}^{{w}_{2}}$,²⁰ we find that ${\mathrm{det}}(\tilde{W})$ is nontrivial. Projecting the Wilson loop operator into the p_z orbital basis, we observe that the corresponding nested Wilson loop spectrum is trivial. On the other hand, within the s, p_x, and p_y orbital basis, the nested Wilson loop spectrum is nontrivial, which confirms that the higher-order band topology originates from the core electronic levels.

Filling anomaly and corner charges

The higher-order band topology of MGD with w₂ = 1 induces a pair of anomalous corner states.³ To see this, we consider a finite-size structure invariant under P symmetry. When the system is chiral symmetric, there are two robust zero-energy in-gap states related by P whose locations are fixed at M_y-invariant corners. One way to see this is to use the topological classification proposed by Geier et al.⁹ When a 2D system has commuting mirror and chiral symmetries, it belongs to BDI${}^{{M}_{++}}$ class that supports nontrivial topology, and the corresponding insulator has zero modes at mirror invariant corners. In MGD, the M_y symmetry commutes with chiral symmetry since M_y symmetry does not mix different sublattices. On the other hand, when chiral symmetry and M_x symmetry anti-commute, the system belongs to BDI${}^{{M}_{+-}}$ class, which is topologically trivial, so there are no protected corner charges at M_x invariant corners.

To be specific, we can find the location of corner states by studying the low-energy Hamiltonian of MGD near the Γ point:

$$H({k}_{x},{k}_{y})=\hslash v({k}_{x}{\Gamma }_{1}+{k}_{y}{\Gamma }_{2})+\Delta {\Gamma }_{3},$$

(6)

where Γ₁ = σ_xτ_x, Γ₂ = σ_yτ_x, Γ₃ = τ_z, and Γ₄ = σ_zτ_x. Corresponding symmetry representations are T = σ_xτ_zK, P = τ_z, M_x = σ_xτ_z, M_y = σ_x, and chiral symmetry S is described by S = Γ₅ = τ_y. Considering a disk-shaped geometry, the low-energy Hamiltonian in real space becomes $H=-i{\tilde{\Gamma }}_{1}(\theta ){\partial }_{r}-i{\tilde{\Gamma }}_{2}(\theta ){r}^{-1}{\partial }_{\theta }+\Delta (r){\Gamma }_{3}$, where ${\tilde{\Gamma }}_{1}(\theta )\,=\,{\mathrm{cos}}\theta {\Gamma }_{1}\,+\,{\mathrm{sin}}\theta {\Gamma }_{2}$ and ${\tilde{\Gamma }}_{2}(\theta )=-{\mathrm{sin}}\theta {\Gamma }_{1}\,+\,{\mathrm{cos}}\theta {\Gamma }_{2}$. We assume that Δ(r) > 0 inside the bulk and Δ(r) < 0 outside the bulk. The edge Hamiltonian is obtained by applying the projection operator $P(\theta )=\frac{1}{2}(1+i{\tilde{\Gamma }}_{1}(\theta ){\Gamma }_{3})$, leading to ${H}_{\text{edge}}=-i{R}^{-1}{\tilde{\sigma }}_{z}{\partial }_{\theta }$ where ${\tilde{\sigma }}_{z}=P{\tilde{\Gamma }}_{2}(\theta ){P}^{-1}$. On the edge, position-dependent mass terms that open the edge gap are allowed. Mass terms that are commuting with P(θ) and anti-commuting with ${\tilde{\Gamma }}_{2}(\theta )$ and S are ${H}_{m}={m}_{4}(R,\theta ){\Gamma }_{4}+i{m}_{25}(R,\theta ){\tilde{\Gamma }}_{2}{\Gamma }_{5}$. The corresponding edge Hamiltonian is ${H}_{\text{edge}}=-i{R}^{-1}{\tilde{\sigma }}_{z}{\partial }_{\theta }+{\tilde{m}}_{1}(R,\theta ){\tilde{\sigma }}_{x}$, where ${\tilde{m}}_{1}={m}_{4}+{m}_{25}$, ${\tilde{\sigma }}_{x}=P{\Gamma }_{4}{P}^{-1}=Pi{\tilde{\Gamma }}_{2}{\Gamma }_{5}{P}^{-1}$. Here, PT symmetry requires ${\tilde{m}}_{1}(\theta +\pi )=-{\tilde{m}}_{1}(\theta )$ while M_y requires ${\tilde{m}}_{1}(x,y)=-{\tilde{m}}_{1}(x,-y)$. Thus, the mass gap is closed at M_y-invariant corners leading to localized corner charges.

Due to the degeneracy of the two in-gap states, the half-filling condition cannot be satisfied as long as P is preserved, which is known as the filling anomaly.^5,11,12 When chiral symmetry is broken, the two degenerate in-gap states can be shifted from zero energy and even merged to the valence (conduction) bands generating a single hole (electron) at half-filling. Namely, the presence of an odd number of holes (electrons) in the valence (conduction) band is the manifestation of the higher-order band topology of generic 2D HOTIs lacking chiral symmetry.^5,16,17,18 To resolve the filling anomaly while keeping P, one needs to add an extra electron (or hole), which leads to the accumulation or depletion of a half-integral charge at each M_y-invariant corner.

Explicitly, let us illustrate how the number of states can be counted in a finite-size HOTI with w₂ = 1 following ref. ⁵ When the total number of states is N and chiral symmetry exists, there are N∕2 − 1 occupied and unoccupied states, respectively, since two in-gap states appear at zero-energy. When chiral symmetry is broken while inversion is preserved, the two in-gap states related by inversion are absorbed into either the valence or the conduction bands simultaneously. Then the number of states below the gap is n_gap = N∕2 ± 1. Therefore n_gap = ($\frac{N}{2}+$ an odd integer) when w₂ = 1, whereas n_gap = ($\frac{N}{2}+$ an even integer) when w₂= 0. When the valence (conduction) band of a system with w₂ = 1 is completely occupied (unoccupied) by adding electrons or holes additionally, a half-integral electric charge should be accumulated (or depleted) at two corners related by inversion.

For the finite-size MGD composed of n × n unit cells shown in Fig. 2a, the total number of states are N = 72n² when 2s, 2p_x, 2p_y, and 2p_z orbitals are considered. In MGD, there is a single sp electron localized at every edge of a hexagonal unit cell. If an edge of a unit cell is not shared with an adjacent unit cell, the relevant sp orbital becomes a half-filled zero-energy non-bonding state when chiral symmetry exists. When we consider the geometry of a finite-size structure shown in Fig. 2a, MGD has 8n − 2 non-bonding states along the boundary. When chiral symmetry is broken, 8n − 2 states can be merged into the valence band, so that ${n}_{\text{gap}}=36{n}^{2}+4n-1=\frac{N}{2}+4n-1$. This is what is observed in first-principles calculations and the corresponding density of states is shown in Fig. 3. Thus, when all the states below the gap are occupied, and thus the filling anomaly is lifted, a half-integral charge (modulo an integer charge) is accumulated at two M_y-invariant corners. This can be easily understood because an odd number of non-bonding states exists only at the M_y-invariant corners (See Fig. 2a).

Fig. 2: Geometry of a finite-size structure and corner charges.

a A finite-size structure of MGD preserving P and M_y symmetries designed to observe the higher-order band topology protected by P symmetry. The red horizontal line indicates the M_y-invariant line. b Electron density distribution for a finite-size MGD composed of 9 × 9 unit cells where hydrogen atoms are attached at every corner except at the two M_y-invariant corners. To resolve the filling anomaly, we add one electron to MGD additionally and fill the valence band completely. The accumulated or depleted charges with respect to the half-filled configuration are plotted. Here, a half-integral charge accumulated at each M_y-invariant corner appears due to the nontrivial 2D topological invariant w₂ = 1.

Fig. 3: Filling anomaly and density of states.

a A schematic figure displaying the energy spectrum of a finite-size MGD composed of n × n unit cells, where the total number of states is N = 72n². Here, a thick gray strip near the Fermi level E_F indicates the non-bonding states arising from carbon atoms along the boundary. The upper panel corresponds to the case at half-filling where an odd number of holes is below the gap, which demonstrates the higher-order band topology with w₂ = 1. The lower panel shows the case when two hydrogen atoms are added to the finite-size MGD. Here, the total number of states $N^{\prime}$ increases by two ($N^{\prime} =N+2$). The hybridization between two hydrogen states and two non-bonding states generates two bonding and two anti-bonding states, so that the number of states below the gap changes from $\frac{N}{2}+4n-1$ to $\frac{N^{\prime} }{2}+4n-2$. Then, w₂ also changes from 1 to 0. Thus, to maintain the value of w₂, the number of hydrogen atoms attached for passivation should be an integer multiple of four. b Density of states (DOS) of a finite-size MGD without hydrogen passivation. Here, the carbon atom at each corner has a non-bonding state. There are (4n − 1) holes below the gap at half-filling. The green color indicate the gapped region. c DOS of a finite-size MGD with hydrogen passivation where (8n − 4) hydrogen atoms are attached along the boundary. Here, only two carbon atoms at M_y-invariant corners have non-bonding states. There is a single hole below the gap at half-filling. Since both systems in b and c have an odd number of holes, both are HOTIs with w₂ = 1.

To observe the corner states more clearly, we add a hydrogen atom to the carbon atom at every corner, except at two M_y-invariant corners, and keep P and M_y symmetries. Then the number of added hydrogen atoms is an integer multiple of four, which does not alter the band topology of MGD. On the other hand, if two hydrogen atoms are added at the M_y-invariant corners, the total number of bands becomes $N^{\prime} =N+2$, while n_gap remains the same as described in Fig. 3a. Then ${n}_{\text{gap}}=\frac{N^{\prime} }{2}+4n-2$ and w₂ becomes trivial. Therefore, the number of added hydrogen atoms should be an integer multiple of four to keep w₂ unchanged. Maximally, (8n − 4) hydrogen atoms can be added to the finite-size MGD with two non-bonding states remaining at M_y-invariant corners. The corresponding density of states is shown in Fig. 3c. Here, a half electric charge is accumulated at each M_y-invariant corner when the states below the gap are fully occupied as shown in Fig. 2b.

Wannier function description

Here, we show that the higher-order band topology is associated with the chemical bonding of the s, p_x, and p_y orbitals across the unit cell boundary. In other words, Wannier functions are located away from the atomic sites and cannot be moved to the atomic sites by a continuous deformation that preserves inversion symmetry. Here, we follow the idea proposed by van Miert and Ortix³⁵ to analyze the Wannier centers in inversion-symmetric systems.

Let us take the inversion-invariant unit cell shown in Fig. 1a whose center is located at the position m. There are four Wyckoff positions in a unit cell labeled as A, B, C, and D—located at 0, $\frac{1}{2}{{\bf{a}}}_{{\bf{1}}}$, $\frac{1}{2}{{\bf{a}}}_{{\bf{2}}}$, and $\frac{1}{2}{{\bf{a}}}_{{\bf{1}}}+\frac{1}{2}{{\bf{a}}}_{{\bf{2}}}$ from the center of the unit cell, respectively—which are invariant under inversion up to lattice vectors a₁ and a₂. The symmetric Wannier functions centered at Wyckoff positions transform under inversion as follows:

$$\begin{array}{lcc}{w}_{A\pm ,{\bf{m}}}({\bf{-x}})\,=\,\pm {w}_{A\pm ,-{\bf{m}}}({\bf{x}}),\\ {w}_{B\pm ,{\bf{m}}}({\bf{-x}})\,=\,\pm {w}_{B\pm ,-{\bf{m}}-{{\bf{a}}}_{{\bf{1}}}}({\bf{x}}),\\ {w}_{C\pm ,{\bf{m}}}({\bf{-x}})\,=\,\pm {w}_{C\pm ,-{\bf{m}}-{{\bf{a}}}_{{\bf{1}}}-{{\bf{a}}}_{{\bf{2}}}}({\bf{x}}),\\ {w}_{D\pm ,{\bf{m}}}({\bf{-x}})\,=\,\pm {w}_{D\pm ,-{\bf{m}}-{{\bf{a}}}_{{\bf{2}}}}({\bf{x}}),\end{array}$$

(7)

where w_α±,m(x) denotes a Wannier function with the coordinate x, centered at the Wyckoff position α of the unit cell with the center at m. The sign ± indicates eigenvalues of the inversion operation with respect to the Wyckoff position. We define the number of symmetric Wannier functions centered at each Wyckoff position with inversion eigenvalues ±1 as N_A±, N_B±, N_C±, and N_D±. Let us note that two Wannier functions with inversion eigenvalues +1 and −1 at A can be expressed as a linear combination of two Wannier functions with inversion eigenvalues +1 and −1 at any other Wyckoff position. Thus, only the difference ν_α = N_α+ − N_α− at each Wyckoff position α is a well-defined ${\mathbb{Z}}$ invariant³⁵.

$$\begin{array}{ccc}{\nu }_{A}&=&{N}_{A+}-{N}_{A-},\\ {\nu }_{B}&=&{N}_{B+}-{N}_{B-},\\ {\nu }_{C}&=&{N}_{C+}-{N}_{C-},\\ {\nu }_{D}&=&{N}_{D+}-{N}_{D-},\end{array}$$

(8)

which are related to the inversion eigenvalues at TRIMs in the Brillouin zone as³⁵

$$\begin{array}{lcc}{\nu }_{A}\,=\,-{\Gamma }_{1}-\frac{1}{2}{\Gamma }_{-1}+\frac{1}{2}{M}_{-1}+\frac{1}{2}{X}_{-1}+\frac{1}{2}{Y}_{-1}=-1,\\ {\nu }_{B}\,=\,\frac{1}{2}{\Gamma }_{-1}-\frac{1}{2}{M}_{-1}+\frac{1}{2}{X}_{-1}-\frac{1}{2}{Y}_{-1}=-1,\\ {\nu }_{C}\,=\,\frac{1}{2}{\Gamma }_{-1}-\frac{1}{2}{M}_{-1}-\frac{1}{2}{X}_{-1}+\frac{1}{2}{Y}_{-1}=-1,\\ {\nu }_{D}\,=\,\frac{1}{2}{\Gamma }_{-1}+\frac{1}{2}{M}_{-1}+\frac{1}{2}{X}_{-1}-\frac{1}{2}{Y}_{-1}=-1.\end{array}$$

(9)

In monolayer graphdiyne (MGD), we find that ν_A = ν_B = ν_C = ν_D = − 1 from the inversion eigenvalues at TRIM. In the viewpoint of chemistry, ∣ν_A∣ = 1 comes from the electrons on the benzene ring around the unit cell center, ∣ν_B∣ = ∣ν_C∣ = ∣ν_D∣ = 1 originates from the sp bonding across the unit cell boundary.

To relate these indices with the two-dimensional topological invariant w₂ defined before, we define electric dipole moments pi = 1,2 and an electric quadrupole moment q₁₂ for the unit cell as follows,

$$\begin{array}{lcc}{p}_{x}\,\equiv\,\sum _{\begin{array}{c}i\in {\rm{occ}},\\ {\rm{unit}}\ {\rm{cell}}\end{array}}{X}_{i}=\frac{1}{2}({N}_{B}+{N}_{C})=0\ {\rm{mod}}\ 1,\\ {p}_{y}\,\equiv\,\sum _{\begin{array}{c}i\in {\rm{occ}},\\ {\rm{unit}}\ {\rm{cell}}\end{array}}{Y}_{i}=\frac{1}{2}({N}_{C}+{N}_{D})=0\ {\rm{mod}}\ 1,\\ {q}_{xy}\,\equiv\,\sum _{\begin{array}{c}i\in {\rm{occ}},\\ {\rm{unit}}\ {\rm{cell}}\end{array}}{X}_{i}{Y}_{i}=\frac{1}{4}{N}_{C}=\frac{1}{4}\ {\rm{mod}}\ \frac{1}{2},\end{array}$$

(10)

where X_i and Y_i are the a₁ and a₂ components of the Wannier center for the ith occupied Wannier function. p_x and p_y are nothing but the polarizations along a₁ and a₂, respectively. w₂ is related to p_x,y, q_xy via w₂ = 4(p_xp_y − q_xy) mod 2,^13,20 so

$$\begin{array}{lcc}{w}_{2}\,=\,({N}_{B}+{N}_{C})({N}_{C}+{N}_{D})-{N}_{C}\\ \,=\,{N}_{B}{N}_{C}+{N}_{C}{N}_{D}+{N}_{D}{N}_{B}\ {\rm{mod}}\ 2\\\,=\,{\nu }_{B}{\nu }_{C}+{\nu }_{C}{\nu }_{D}+{\nu }_{D}{\nu }_{B}\ {\rm{mod}}\ 2,\end{array}$$

(11)

because ν = N₊ − N₋ = N₊ + N₋ − 2N₋ = N − 2N₋ = N mod 2. This shows that ${w}_{2}=1\ {\rm{mod}}\ 2$ when ∣ν_B∣ = ∣ν_C∣ = ∣ν_D∣ = 1 due to the sp bonding across the unit cell boundary.

Discussion

Since MGD has w₂ = 1, when layers of MGD are stacked vertically, the resulting 3D insulator with PT symmetry becomes a 3D weak topological insulator, dubbed a 3D weak Stiefel-Whitney insulator,²⁰ when inter-layer coupling is weak. When inter-layer coupling becomes large enough, an accidental band crossing can happen at a TRIM at which a pair of nodal lines with Z₂ monopole charge are created. The existence of two nodal lines at k_z = ± k_c(k_c > 0) was predicted in ABC-stacked graphdiyne,¹⁹ and their Z₂ monopole charge was also confirmed²⁰ based on first-principles calculations. In ABC-stacked graphdiyne, the 2D subspace with a fixed k_z carries w₂ = 1 (w₂ = 0) when ∣k_z∣ < k_c (k_c < ∣k_z∣ < π) since the band inversion is happened at k = (0, 0, π). Since a 2D HOTI with w₂ = 1 possesses corner charges, similar corner charges are also expected in ABC-stacked graphdiyne in the subspace with a fixed k_z where w₂ = 1,³ which leads to hinge modes of the 3D structure shown in Fig. 4a.

Fig. 4: Hinge modes of ABC-stacked graphdiyne with monopole nodal lines.

a A schematic figure describing the 3D geometry of ABC-stacked graphdiyne which is finite in the xy plane but periodic in the vertical direction. The structure must preserve C_2x symmetry to carry hinge states. b Hinge modes, located in the region k_c < ∣k_z∣ < π, obtained by tight-binding Hamiltonian using only p_z orbitals. In real materials including all core electronic levels, the hinge modes should appear in the region 0 < ∣k_z∣ < k_c where w₂ = 1.

To demonstrate the hinge modes of ABC-stacked graphdiyne, we study a tight-binding model by using only p_z orbitals to reduce numerical costs. The energy spectrum of a finite-size system with PT and C_2x symmetries is shown in Fig. 4b, which clearly shows the presence of hinge modes in the momentum region where w₂ = 1 as expected. When only p_z orbitals are included, although the low-energy band structure from the tight-binding model is consistent with that from first-principles calculations, the topological property is different. Namely, the tight-binding model predicts w₂ = 0 (w₂ = 1) when ∣k_z∣ < k_c (k_c < ∣k_z∣ < π), which is opposite to the result from the first-principles calculations. Such a discrepancy can be remedied when the core electronic levels are included in the tight-binding calculation, which is confirmed by separate calculations.

To sum up, we have shown that the nontrivial band topology of a 2D HOTI lacking chiral symmetry can be revealed by using the filling anomaly, which can be explicitly demonstrated by counting the number of states of a finite-size system. This idea can be straightforwardly generalized to any d-th order topological insulator in d-dimensions hosting zero-dimensional corner states. We believe that our theoretical study provides a general way to extend the scope of HOTI materials to a wider class.

Methods

DFT calculations

We have carried out the first-principles density-functional theory calculations using the Vienna ab initio simulation package.³⁶ The projector augmented-wave method³⁷ and the exchange-correlation functional of the generalized gradient approximation in the Perdew, Burke, and Ernzerhof³⁸ scheme were used. The self-consistent total energy was evaluated with an 12 × 12 × 1k-point mesh, and the cutoff energy for the plane-wave basis set was 500 eV. For the finite-size structure, we used only Γ point. A single unit cell contains 18 carbon atoms, where a benzene ring is connected to six neighboring ones by linear chains. In the optimized structure, we have L = 9.46 Å, b₁ = 1.43 Å, b₂ = 1.40 Å, b₃ = 1.23 Å, and b₄ = 1.34 Å.

Tight-binding calculation

Low-energy band structure of ABC-stacked graphdiyne can be properly described by p_z orbital basis. The intralayer coupling is described by the nearest neighbor hopping parameters expressed as the transfer integral ${V}_{{\mathrm{pp\pi}} }={V}_{{\mathrm{pp\pi}}}^{0}{e}^{(-R-{a}_{0})/{r}_{0}}$, where ${V}_{{\mathrm{pp\pi}}}^{0}\approx -2.7\,{\mathrm{eV}}$ is the hopping parameter for graphene’s nearest neighbor interaction (a₀ = 0.142 nm and r₀ = 0.09 nm).¹⁹ The interlayer coupling is described by two hopping parameters of two shortest vertical bonds which are expressed as the transfer integral ${V}_{{\mathrm{pp\sigma}}}={V}_{{\mathrm{pp\sigma}}}^{0}{e}^{(-R-{d}_{0})/{r}_{0}}$, where V_ppσ≈ 0.48 eV is the hopping parameter at the interlayer distance of graphite (d₀ = 0.335 nm).¹⁹