Two-dimensional higher-order topology in monolayer graphdiyne

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.


I. 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. 1 Also there are other candidate materials of 3D HOTIs including SnTe with strain along the (100) direction, 2 transition metal dichalcogenides MoTe 2 and WTe 2 , hosting helical hinge states, 3 and also Bi 2−x Sm x Se 3 with chiral hinge states. 4 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. 15 However, in this system, the corner states originate from the charge polarization, 15 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. 13 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 zerodimensional 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][9][10][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 2 , 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 non-trivial w 2 of MGD gives rise to the topological band structure of the corresponding three-dimensional material, ABC stacked graphdiyne, hosting nodal lines with Z 2 monopole charges. 19,20 Due to the nonzero w 2 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 .

A. Topological invariant of monolayer graphdiyne
MGD is an experimentally realized planar carbon system composed of sp 2 -sp carbon network of benzene rings connected by ethynyl chains. 21,22 Fig. 1(a) 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. 1(b) 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.
In general, a 2D P -symmetric spinless fermion system carries three Z 2 topological invariants, w 1x , w 1y , and w 2 . 20 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 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 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 2 is a 2D topological invariant given by the 2D winding number of the sewing matrix G(k) (modulo 2), 6 and measures the higher-order band topology of P -symmetric spinless fermion systems. 5,6,13 w 2 can be determined by where N − occ (Γ 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 2 = 1. One interesting feature of the MGD band structure is that we get w 2 = 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 2 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 where , and its spectrum is gauge-invariant. 33 Computed along the k 1 direction parallel to the reciprocal lattice vector G 1 from (k 1 , k 2 ), the set of Wilson loop eigenvalues {e iν i (k 2 ) } indicates the position of Wannier centers at given k 2 , and the corresponding total charge polarization is given by . The Wilson loop spectrum of MGD shows that the charge polarization is zero in both k 1 and k 2 directions, which is consistent with the parity eigenvalue analysis ( Fig. 1(c)). Also, the single linear-crossing at (k y , ν) = (0, π) indicates w 2 = 1. 20 Projecting the Wilson loop operator into the p z (s, p x , and p y ) orbital basis, we observe w 2 = 0 (w 2 = 1), which confirms that the higher-order band topology originates from the core electronic levels. The orbital dependence of w 2 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 P T 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 2 is diagonalized as then, using the set of eigenvectors corresponding to the sub-band centered at ν = π, the nested Wilson loop operatorW along the G 2 vector is defined as where 20 we find that det(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.

B. Filling anomaly and corner charges
The higher-order band topology of MGD with w 2 = 1 induces a pair of anomalous corner states. 3 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.. 9 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: where Considering a disk-shaped geometry, the low-energy Hamiltonian in real space becomes H = 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 (θ) = 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Γ 2 (θ) and S are H m = m 4 (R, θ)Γ 4 + im 25 (R, θ)Γ 2 Γ 5 . The corresponding edge Hamiltonian 25 ,σ x = P Γ 4 P −1 = P iΓ 2 Γ 5 P −1 . Here P T symmetry requiresm 1 (θ + π) = −m 1 (θ) while M y requiresm 1 (x, y) = −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 2 = 1 following Ref. 5. 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 zeroenergy. 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 = ( N 2 + an odd integer) when w 2 = 1, whereas n gap = ( N 2 + an even integer) when w 2 = 0. When the valence (conduction) band of a system with w 2 = 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. 2(a), the total number of states are N = 72n 2 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. 2(a), 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 gap = 36n 2 + 4n − 1 = 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. 2(a).] 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 = N +2, while n gap remains the same as described in Fig. 3(a). Then n gap = N 2 +4n−2 and w 2 becomes trivial. Therefore, the number of added hydrogen atoms should be an integer multiple of four to keep w 2 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. 3(c). 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. 2(b).

C. 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 35 to analyze the Wannier centers in inversion-symmetric systems.
Let us take the inversion-invariant unit cell shown in Fig. 1(a) where w α±,m (x) denotes a Wannier function with the coordinate x, centered at the Wyckoff po- which are related to the inversion eigenvalues at TRIMs in the Brillouin zone as 35 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 2 defined before, we define electric dipole moments p i=1,2 and an electric quadrupole moment q 12 for the unit cell as follows, where X i and Y i are the a 1 and a 2 components of the Wannier center for the ith occupied Wannier function. p x and p y are nothing but the polarizations along a 1 and a 2 , respectively. w 2 is related to p x,y , q xy via w 2 = 4(p x p y − q xy ) mod 2, 13,20 so This shows that w 2 = 1 mod 2 when |ν B | = |ν C | = |ν D | = 1 due to the sp bonding across the unit cell boundary.

III. DISCUSSION
Since MGD has w 2 = 1, when layers of MGD are stacked vertically, the resulting 3D insulator with P T symmetry becomes a 3D weak topological insulator, 20 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 2 monopole charge are created. The existence of two nodal lines at k z = ±k c (k c > 0) was predicted in ABC-stacked graphdiyne, 19 and their Z 2 monopole charge was also confirmed 20 based on first-principles calculations. In ABC-stacked graphdiyne, the 2D subspace with a fixed k z carries w 2 = 1 (w 2 = 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 2 = 1 possesses corner charges, similar corner charges are also expected in ABC-stacked graphdiyne in the subspace with a fixed k z where w 2 = 1, 3 which leads to hinge modes of the 3D structure shown in Fig. 4(a).
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 P T and C 2x symmetries is shown in Fig. 4(b), which clearly shows the presence of hinge modes in the momentum region where w 2 = 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 2 = 0 (w 2 = 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.

A. DFT calculations
We have carried out the first-principles density-functional theory calculations using the Vienna ab initio simulation package. 36 The projector augmented-wave method 37

B. 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 7eV is the hopping parameter for graphene's nearest neighbor interaction (a 0 = 0.142nm and r 0 = 0.09nm). 19 The interlayer coupling is described by two hopping parameters of two shortest vertical bonds which are expressed as the transfer integral V ppσ = V 0 ppσ e (−R−d 0 )/r 0 , where V ppσ ≈ 0.48eV is the hopping parameter at the interlayer distance of graphite (d 0 = 0.335nm). 19

V. NOTES
During the preparation of our manuscript, we have found a related work. 39 In [39], the authors also observed HOTI phase in MGD. However, some of their theoretical arguments are inconsistent with our results. For instance, in [39], it is stated that crystalline symmetries are not essential for the protection of corner states, which is in sharp contrast to our paper where the role of P symmetry is emphasized. Also, in [39], it is mentioned that double band inversion of p z orbitals at the Γ point is the origin of the HOTI phase. This is different from our result in that higher-order band topology originates from the core levels derived from s, p x , p y orbitals.

VI. DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

VII. ACKNOWLEDGEMENTS
We acknowledge the helpful discussions with Dongwook Kim and Youngkuk Kim.

VIII. COMPETING INTERESTS
The Authors declare no competing interests.    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 2 = 1.