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 pz 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 px,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. 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 tightbinding 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. 13 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 kthorder 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 halfintegral 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 welldescribed 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 2 of MGD gives rise to the topological band structure of the corresponding threedimensional 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 .

RESULTS
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 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.
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 W ðk1þ2π;k2Þ ðk1;k2Þ ¼ lim Þi, 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 fe iνiðk2Þ g indicates the position of Wannier centers at given k 2 , and the corresponding total charge polarization is given by 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. 1c). 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 PT symmetry, the Wilson loop spectrum can be divided into two subbands 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 Since the determinant of the nested Wilson loop operator detðWÞ is identical to ðÀ1Þ w2 , 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.
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 finitesize 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 Γ 1 = σ x τ x , Γ 2 = σ y τ x , Γ 3 = τ z , and Γ 4 = σ z τ x . Corresponding symmetry representations are T = σ x τ z K, P = τ z , M x = σ x τ z , M y = σ x , and chiral symmetry S is described by S = Γ 5 = τ y . Considering a disk-shaped geometry, the low-energy Hamiltonian in real space becomes H ¼ ÀiΓ 1 ðθÞ∂ r À iΓ 2 ðθÞr À1 ∂ θ þ ΔðrÞΓ 3 , whereΓ 1 ðθÞ ¼ cosθΓ 1 þ sinθΓ 2 andΓ 2 ðθÞ ¼ ÀsinθΓ 1 þ cosθΓ 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 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 is H edge ¼ ÀiR À1σ Here, PT 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 higherorder 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 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 = ( 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. 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 2 = 1.
For the finite-size MGD composed of n × n unit cells shown in Fig. 2a, 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. 2a, MGD has 8n − 2 nonbonding 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 yinvariant 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).
To observe the corner states more clearly, we add a hydrogen atom to the carbon atom at every corner, except at two M yinvariant 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 0 ¼ N þ 2, while n gap remains the same as described in Fig. 3a. Then n gap ¼ N 0 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. 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 35 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, 1 2 a 1 , 1 2 a 2 , and 1 2 a 1 þ 1 2 a 2 from the center of the unit cell, respectivelywhich are invariant under inversion up to lattice vectors a 1 and a 2 . The symmetric Wannier functions centered at Wyckoff positions transform under inversion as follows: w A ± ;m ðÀxÞ ¼ ± w A ± ;Àm ðxÞ; w B ± ;m ðÀxÞ ¼ ± w B ± ;ÀmÀa1 ðxÞ; w C ± ;m ðÀxÞ ¼ ± w C ± ;ÀmÀa1Àa2 ðxÞ; w D ± ;m ðÀxÞ ¼ ± w D ± ;ÀmÀa2 ðxÞ; 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 α− 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 yinvariant 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 yinvariant corner appears due to the nontrivial 2D topological invariant w 2 = 1.
at each Wyckoff position α is a well-defined Z invariant 35 .
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 pi = 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.

DISCUSSION
Since MGD has w 2 = 1, when layers of MGD are stacked vertically, the resulting 3D insulator with PT symmetry becomes a 3D weak 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 2 = 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 0 increases by two (N 0 ¼ 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 N 2 þ 4n À 1 to N 0 2 þ 4n À 2. Then, w 2 also changes from 1 to 0. Thus, to maintain the value of w 2 , 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 2 = 1.
topological insulator, dubbed a 3D weak Stiefel-Whitney insulator, 20 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. 4a.
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 2 = 1 as expected. When only p z orbitals are included, although the lowenergy 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.

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 and the exchange-correlation functional of the generalized gradient approximation in the Perdew, Burke, and Ernzerhof 38 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 = 1.43 Å, b 2 = 1.40 Å, b 3 = 1.23 Å, and b 4 = 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 ppπ ¼ V 0 ppπ e ðÀRÀa0Þ=r0 , where V 0 ppπ % À2:7 eV is the hopping parameter for graphene's nearest neighbor interaction (a 0 = 0.142 nm and r 0 = 0.09 nm). 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Àd0Þ=r0 , where V ppσ ≈ 0.48 eV is the hopping parameter at the interlayer distance of graphite (d 0 = 0.335 nm). 19 Fig. 4 Hinge modes of ABC-stacked graphdiyne with monopole nodal lines. a A schematic figure describing the 3D geometry of ABCstacked 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.