Higher-Order Topological Crystalline Insulating Phase and Quantized Hinge Charge in Topological Electride Apatite

In higher-order topological insulators, bulk and surface electronic states are gapped, while there appear gapless hinge states protected by spatial symmetry. Here we show by ab initio calculations that the La apatite electride is a higher-order topological crystalline insulator. It is a one-dimensional electride, in which the one-dimensional interstitial hollows along the $c$ axis support anionic electrons, and the electronic states in these one-dimensional channels are well approximated by the one-dimensional Su-Schrieffer-Heeger model. When the crystal is cleaved into a hexagonal prism, the 120$^\circ$ hinges support gapless hinge states, with their filling quantized to be 2/3. This quantization of the filling comes from a topological origin. We find that the quantized value of the filling depends on the fundamental blocks that constitute the crystal. The apatite consists of the triangular blocks, which is crucial for giving nontrivial fractional charge at the hinge.

Studies on the topology of band structures in k-space reveal various nontrivial phases of insulators and semimetals [1][2][3]. In recent years, higher-order topological phases have attracted high attention, where higher-order topological states appear due to nontrivial bulk topology. There are two types of higher-order topological phases: those supporting currents [4][5][6][7][8] and those supporting charges [9][10][11][12][13]. One of the realistic materials having chiral hinge states is bismuth. Such chiral hinge states have been confirmed by the scanning tunneling microscope (STM) [7]. The other example is Bi 4 Br 4 , where the bulk is completely insulating unlike bismuth [8]. On the other hand, there has been no proposal for a realistic threedimensional material having higher-order topological states supporting quantized charges.
Here, we focus on electride as a system that can realize a higher-order topological phase. In electrides, electrons enter the voids surrounded by cations and are responsible for stabilizing the structure as anions [14][15][16][17]. Interstitial electrons are not strongly stabilized by the electric field from positive ions, as compared to electrons belonging to atomic orbitals stabilized by the strong electric field from the nuclei. Therefore, the work function of an electride is generally small [15,18], and electrides are actively studied in the chemistry field such as catalyst [19]. Compared with covalent crystals originating from p orbitals, there is less structural instability in electrides when making the cleavage plane [20].
In recent years, it has been proposed that electrides have a high affinity with topological phases in physics [20]. Since the interstitial electrons lead to a small work function, the bands originating from interstitial electrons appear near the Fermi level. Therefore, electrides are likely to have band inversion. Actually, various topological phases with/without the spinorbit interaction (SOI) are proposed in the electrides [20][21][22][23]. For example, a two-dimensional electride Sc 2 C is a nontrivial insulator having a quantized Zak phase π [20]. Because anionic electrons exist in the two-dimensional regions between Sc 2 C layers, floating topological charges protected by the bulk topology appear on the surface.
In this letter, we propose an apatite electride as a higherorder topological crystalline insulator (TCI). The apatite electride is a one-dimensional electride, having anionic electrons in one-dimensional hollows surrounded by cations [24]. We show that quantized topological charges protected by the C 6 symmetry appear at the hinge. generalized gradient approximation (GGA). We use VASP for the lattice optimization. The energy cutoff is 50 Ry for the calculation by VASP [25]. We use the ab initio code OpenMX for the calculation of the electronic structure [26]. and the energy cutoff for the numerical integrations is 150 Ry for the calculation by OpenMX. The 6 × 8 × 6 regular k-mesh is employed for the bulk. We construct Wannier functions from the Kohn-Sham bands, using the maximally localized Wannier function [27,28]. The density of states on the (0110) surface is calculated by the recursive Green's function method [29]. In the slab (wire) calculations, we tune the on-site potential of the Wannier function at the surface (hinge) to satisfy the charge neutrality condition, respectively. We confirm the effect of the charge neutrality and the long-ranged Coulomb interaction by first-principles calculations (see Supplementary Material [30]). Figure 1(a) shows the crystal structure of the apatite electride A 6 B 4 (SiO 4 ) 6 . Crystals with apatite structure are so stable that they appear in nature. The most well-known one is hydroxyapatite Ca 10 (PO 4 ) 6 (OH) 2 , that is the main constituent of bones and teeth. It has a hexagonal crystal structure having a cleavage plane of hexagonal shape around the (0001) direction [31]. There are no atoms on the boundary of the unit cell except for the top and bottom planes. Since the apatite has a one-dimensional structure, it is used in the study of ionic conductors. Experimentally, the apatite electride A 6 B 4 (SiO 4 ) 6 can be synthesized from A 6 B 4 (SiO 4 ) 6 O 2 by releasing two oxygen atoms located in a region surrounded by cations. The apatite electride has a one-dimensional hollow along the zaxis as shown in Fig. 1(b). The positive ions in the hollows are located at z = 1/4 and 3/4, where the lattice constant along the z-direction is set to be unity. Since the Wannier orbitals in the interstitial region are s-orbital-like and there is no nucleus inside the interstitial region, the SOI of the interstitial electron is very weak. In the following, the apatite refers to the apatite electride, unless otherwise specified. Figure 1(c) is the spatial symmetry of the apatite. The space group of the apatite is No. 176 (P6 3 /m). One of the space inversion centers lies at (x, y, z) = (0, 0, 0) of the unit cell. It also has C 6 screw rotation symmetry with 1/2 translation along the c-axis with its screw axis at the hollow (x, y) = (0, 0). The unit cell consists of two triangular prisms, which are transformed to each other by space inversion. Each triangular prism region has the C 3 rotational axis around the c-axis at the center. The system also has mirror symmetry perpendicular to the c-axis.
Figures 2(a) and (b) show electronic band structures of La and Y apatites, respectively. The corresponding Brillouin zone is shown in Fig. 2(c). While the A site is La 3+ /Y 3+ , the B site is completely dissolved with La 3+ /Y 3+ and another positive ion or void and has +2.5e charge on average. In our numerical calculation, such mixed states are handled by the virtual crystal approximation. The La apatite is an insulator, while the Y apatite has a semimetallic band structure. Reflecting a quasi one-dimensional crystal structure, they have a large band dispersion along the k z direction. Both the valence band and the conduction band near the Fermi level originate from the interstitial electron as shown in Figs. 2 (d)-(f). The corresponding Wannier bands for the interstitial electron are shown in Fig. 2 [34,35] are topological electrides in which either the bottom of the conduction band or the top of the valence band originates from the interstitial orbital. No topological material has been known so far in which the band inversion occurs between the interstitial orbitals, and apatite is a first example of such a case. We discuss the irreducible represen- We show our results on surface states of the La apatite. We consider the slab with the (0110) surfaces as shown in Fig 3(a), which preserves space inversion and C 2 symmetry. Figure 3(b) is the surface states calculated using the recursive Green's function. The effect of the lattice optimization is shown in Supplementary Material [36]. The surface states are gapped similarly to the bulk states. This gapped surface state means that the Zak phase of the bulk is zero [3,38,39]. The detail of the Zak phase is discussed in Supplementary Material [40]. Two spinless states are occupied in the hollow in the bulk, and one spinless state is occupied at the surface ( Fig. 3(a)), excluding the double degeneracy due to the spins. Because the filling of the surface bands is an integer, the surface bands are gapped, as seen from Fig. 3(c).
Next we show the hinge band structure of the La apatite. Here, we consider a hexagonal prism, which preserves the space inversion and the C 6 symmetry ( Fig. 4(a)). Figures 4(b),(c) shows the one-dimensional band structure of the La apatite prism. While the band structure is gapped in the bulk and the surface, as is shown in Fig. 3(b), topological gapless hinge states appear at the Fermi level. Such hinge states are protected by nontrivial topology of the bulk La apatite, and is characterized as a higher-order topological insulator protected by a rotational symmetry, proposed in Ref. [12]. In Ref. [12], to describe this higher-order topology, integer topological invariants for the high symmetry Π point is defined in terms of the C n eigenvalues at high-symmetry points Π compared to those at the Γ point where C n is a local symmetry at the high-symmetry point Π, Π (n ) p (= e 2πi(p−1)/n ) is the C n eigenvalue labeled by p at the high-symmetry point denoted by Π, and #Π (n ) p is the number of the occupied bands having the eigenvalue Π (n ) p . In Ref. [12], it was proposed that the topological class of the two-dimensional TCI having the C 6 symmetry is characterized by χ (6) = ([M (2) 1 ], [K (3) 1 ]). By using the connections between the positions of the Wannier orbitals and the irrep at high-symmetry points, the corner charge is given in terms of the topological numbers by [12] Q H(6) where we adopt the convention e = −|e| < 0, following Ref. [12]. The rotational symmetry along the c-axis in the apatite is C 6 screw rotation 6 3 . Since the apatite is insulating, the eigenvalues for the rotational symmetries are same between k z = 0 and k z = π of the high-symmetry points, and we can limit our discussion to the irrpdf at k z = 0 for characterization of the hinge charge. Therefore, the 6 3 helically symmetric system is essentially the same as the 6-fold symmetric system in terms of the topological charges. Nonetheless, this formula (2) does not apply to our system. Equation (2) is derived for a crystal consisting of hexagonal blocks ( Fig. 4(d)), with the sixfold rotation axis at the center of the hexagon, and Eq. (2) represents a corner charge for a 2D crystal with a hexagonal shape, composed of the hexagonal blocks ( Fig. 4(g)). Meanwhile, in the present case, if we focus on the crystal shape along the ab plane, the fundamental unit is a regular triangle, rather than a hexagon. The whole crystal forming a hexagonal prism is composed of the triangular blocks, and the unit cell consists of two triangular blocks ( Fig. 4(e)). This choice of the fundamental units of the crystal makes a difference in the corner charge. In the present case, the Wannier orbital is at the hollow sites, located at the Wyckoff position 1a. In this case the representations are trivial, and the same between the high-symmetry points, and [M (2) 1 ] = 0 and [K (3) 1 ] = 0. This leads to an absence of the corner charge if we use Eq. (2) which is natural because the Wyckoff position 1a is at the center of the hexagon. Nonetheless, it is not valid, because the fundamental blocks are triangular and this Wyckoff position 1a is at the corners of the triangle (Fig. 4(f)), giving rise to a nonzero quantized corner charge as we show later. Thus, even for trivial representations, one can have a nonvanishing quantized corner charge. For a triangular block, we show maximal Wyckoff positions in Fig. 4(e). For each Wyckoff position, one can calculate the corner charge and irreducible representations at highsymmetry points similarly as in Ref. [12]. Through some calculations whose details are given in Supplementary Material [41], the corner charge for C 6 -symmetric 2D systems is given by where ν is a number of electrons per unit cell, i.e. the number of occupied bands. It means that the filling of the corner states is Q corner /|e|, and it is well defined only when the bulk and the surface are gapped. In the present case, the surface is gapped only when ν is even. This equation is quite different from (2) in that it depends on ν. We now apply this theory to the apatite. We consider a hexagonal prism, having a structure within the ab plane as shown in Fig. 4(a). Let us consider the whole system as a two-dimensional system, and consider only the interstitial electrons. Then we have [M (2) 1 ] = 0, [K (3) 1 ] = 0, and ν = 2, giving rise to Q T(6) corner = 2|e| 3 mod e. From the charge neutrality condition, the filling of the topological hinge state is 2 3 .
To summarize, we show that the apatite electrides are higher-order topological insulators, having topological gap-less hinge states, while the bulk and the surface are gapped. The hinge states reside at interstitial sites forming onedimensional hollows in the electride. In the bulk and the surface, these interstitial sites are integer-filled, but at the 120 • hinge, they are 2/3 filled due to the topological properties of the bulk occupied bands. Realization of the topological bands stemming from the band inversion, together with the floating hinge states at the hinges originates from the characteristics as elelctrides.
This work was supported by JSPS KAKENHI Grant Number 18H03678, and by the MEXT Elements Strategy Initiative to Form Core Research Center (TIES).
[28] I. Souza, N. Marzari, and D. Here the bands are shown as semimetallic bands, and the same argument can apply to insulators as well. If a material in which the valence band and the conduction band are composed of atomic orbitals, as shown in Supplementary Fig. S1, it is not an electride. Then, if a conduction band originating from the interstitial electron is inverted with the valence band, a part of the occupied band originates from the interstitial electrons ( Supplementary  Fig. S1(b)). Such a system is a topological electride, where interstitial electrons are responsible for band inversion. We note that if the region of the band inversion in the k space is small, the occupied bands are mostly atomic orbitals, and the anionic electrons do not largely contribute to the stability of the material. We find that Na 3 Bi is an example of this type, shown in Supplementary Fig. S1(b) (Supplementary Figs. S2(a)-(e)). Na 3 Bi is a well-known Dirac semimetal [S1]. Although the conduction band of Na 3 Bi might seem to originate from the s orbital of Na, it actually originates from the interstitial electrons. From the above scenario, the interstitial states are occupied near the Γ point, and the band crossings in Supplementary Fig. S1(b) form Dirac cones in this Dirac semimetal. As a result, the Fermi-arc surface states, connecting the Dirac points, are floating states originating from the interstitial region. We also found that Ca 3 P 2 is a similar example (Supplementary Figs. S2(f)-(h)). Ca 3 P 2 is a nodal-line semimetal [S2, S3]. As in the case of the La apatite, interstitial electrons are spreading in a one-dimensional hollow surrounded by Ca 2+ cations. In fact, many of the well-known topological materials turn out to be topological electrides.
In the systems shown in Supplementary Fig. S1(c), the structure is stabilized by interstitial electrons. Examples of such a type are Y 2 C (nodal-line semimetal), A 2 B (A=Ca, Sr, Ba, B=As, Sb, Bi) (nodal-line semimetal), HfX (quantum spin Hall), and LaX (X=Cl, Br, I) (quantum anomalous Hall) [S4]. The La apatite is an example of band inversion between two bands, both with interstitial electrons (Supplementary Fig. S1(d)). Sc 2 C (insulator with π Zak phase) also belongs to this type of the topological electride in that the interstitial electron is completely occupied [S4].

Supplementary Material 2. ELECTRONIC BAND STRUCTURE AND IRREDUCIBLE REPRESENTA-TION OF LA-APATITE
Supplementary Figures S3(a),(b) show the electronic band structure of the La apatite with and without two O 2− ions at the hollow, respectively. In the space group P6 3 /m, all the states at k z = π are doubly degenerate. This is physically reasonable because bonding and antibonding states cannot be formed in the c direction at k z = π. In the calculation of the polarization, the Zak phase is zero even if the contribution of the valence bands below −4 eV are included. In fact, even in slab calculations with structural optimization (Supplementary Henceforth, we focus on interstitial electrons near the Fermi level of the La apatite ( Fig. 2(g) in the main text). The irrpdf of the two occupied bands at the Γ point are Γ 1+ and Γ 2+ from the lowest band. Since the occupied bands originate from the interstitial s-like orbital at (x, y, z) = (0, 0, 0), eigenvalues for point-group operations not including a translation such as I and C 3 are always trivial. The irrpdf at the M points of the two occupied bands are M 1+ and M 2+ at all three M points and those at the K-points are K 1 and K 2 at the two K points, due to the C 6 symmetry.

Supplementary Material 3. DETAILS OF THE EFFEC-TIVE MODEL
We construct an effective model of the interstitial electron near the Fermi level of the La apatite ( Fig. 2(g) in the main text), based on the ab initio calculation. The effective model is composed of four equivalent Wannier functions localized along the hollow x = y = 0, and they are located at the Wyckoff positions 4e (0,0,±z 0 ) and (0,0,±z 0 + 1/2) (z 0 = 0.163), where the electric potential is lowered by the six La 3+ ions located at z = 1/4 and 3/4. This model is well approximated by the Su-Schrieffer-Heeger (SSH) model along the c direction ( Supplementary Figs. S4(a),(b)), and has negligible hopping in the ab directions. We note that due to the screw symmetry with a half translation along the c axis, the unit cell for the SSH model is doubled, as shown in Supplementary Fig. S4(a). In this doubled unit cell, the two occupied bands arise from strong covalent bonds at z = 0 and 1/2, because of |t 1 | > |t 2 | in Supplementary Fig. S4(a). This model with Wannier functions at 4e reproduces the bulk bands of apatite well (Fig. 2(g) in the main text), but it is not suitable for calculations on surfaces and hinges, because the Wyckoff positions 4e are on the cleavage planes. Therefore, in our calculations on surfaces and hinges as we show later, we also construct an effective model by putting twelve equivalent Wannier functions for interstitial electrons at the Wyckoff position 12i near the six La 3+ ions. The bulk band structure is shown in the inset of Fig. 2(g) in the main text.

Supplementary Material 4. ZAK PHASE IN APATITE
The surface band is closely related to the polarization P(k ) and the polarization is associated with the Zak phase θ(k ) as P = e 2π θ(k ). To define the Zak phase, we pick up one crystallographic plane, and its normal vector is specified by a reciprocal lattice vector G. Then the wavevector k is decomposed into the components parallel and perpendicular to G. n = G/|G|, k = k + k ⊥ n, k · n = 0. Then the Zak phase is where u n (k) is the periodic part of the bulk Bloch function in the n-th band with the gauge choice u n (k) = u n (k + G)e iG·r , and the sum is taken over the occupied states. The Zak phase is quantized to 0 or π when the system has both space inversion and time-reversal symmetry, and it takes a constant value for insulators. In the La apatite, regarded as a spinless system thanks to the weak SOI, the Zak phase along the [0110] direction is 0 (mod 2π). The Zak phase is related with the parity inversion for the occupied bands. Since the interstitial electron is s-orbital and is centered at the edge of the unit cell, the parities of the two occupied states at the Γ point are both even and those at the M point are both odd. Therefore, the product of parity for the Γ and M points is even in total and the Zak phase is 0. This is consistent with the absence of topological in-gap topological surface states in Fig. 3b. The Zak phase also corresponds to the center of the Wannier function for the occupied band. The centers of the Wannier functions for the occupied bands are (0, 0, 0) and (0, 0, 1/2), which are at boundaries of the unit cell, contributing to the Zak phase along the [0110] by π. Therefore, the Zak phase is 2 · π ≡ 0 (mod 2π).

Supplementary Material 5. DETAIL OF THE CALCU-LATION OF THE CORNER CHARGE
Here, we relate the corner charge with topological indices in TCIs in class AI with C 6 crystalline symmetry, when the crystal is composed of triangular blocks as we discuss in the context of apatites in the main text. The previous study has assumed that the crystal is composed of a hexagonal block, which is invariant under C 6 symmetry [S5]. In our case, however, the crystal is composed of triangular units, and natural termination of the crystal is different from that of crystals with hexagonal blocks. Each unit cell consists of two triangular blocks.
We use a notation of C n -symmetric configurations with Wannier orbitals used in Ref. S5. Let h (n) mW to be a C nsymmetric configuration with m filled bands, having Wannier centers at the maximal Wyckoff positions W. At each Wyck- off position, we put Wannier orbitals following a certain irreducible representation of the local point-group symmetry at W. In the present case for systems with C n symmetry, local symmetry at W is rotational symmetry C n , where n is a factor of n. Thus, Wannier orbitals to be put at Wyckoff positions W have a C n eigenvalue r = e iθ , θ ≡ 2i jπ/n (mod 2π), where j is an integer. Within the unit cell there are n/n locations corresponding to this Wyckoff position with C n symmetry, and therefore the number of occupied bands m is given by m = n/n . Such a configuration is denoted by h (n) mW (mW| θ ). In the present case, since we limit ourselves to systems with time-reversal symmetry, the Wannier orbitals with θ and −θ are always degenerate, where θ 0, π (mod 2π). In such a case, the configurations with θ and −θ are jointly written as h (n) 2m W (mW| ±θ ), where m bands with θ and m bands with −θ are degenerate, forming 2m bands in total. These notations represent configurations for Wannier orbitals, and we use the same notation for the corresponding Hamiltonian. Now we consider all the Wannier configurations with C 6rotation symmetry, both for the cases with hexagonal blocks and with triangular blocks. There are three types of maximal Wyckoff positions: the first one at the Wyckoff position a (Supplementary Fig. S5(a)), the second one at b and b ( Supplementary Fig. S5(b)), and the third one at c, c and c ( Supplementary Fig. S5(c)). The local symmetry at Wyckoff positions a, b and c are C 6 , C 3 and C 2 , respectively, and therefore the possible configurations are h (6) 1a (1a| 0 ), h (6) 1a (1a| π ), h (6) 2a (1a| ±π/3 ), h (6) 2a (1a| ±2π/3 ), h (6) 2b (2b| 0 ), h (6) 4b (2b| ±2π/3 ), h (6) 3c (3c| 0 ), and h (6) 3c (3c| π ). Among the eight configurations, it suffices to consider only the three configurations h (6) 1a (1a| 0 ), h (6) 2b (2b| 0 ), h (6) 3c (3c| 0 ), and they are generators of the algebra considered, as we explain later.

Corner Charge for Wannier Configurations
In the following, we calculate the edge filling and the corner charges for the Wannier configurations, following Ref. S5. To calculate the corner charge, we consider a crystal of finite size, respecting the C n symmetry. In the present case of n = 6 we consider a crystal with a regular hexagonal shape. In Supplementary Figs. S5A, B and C, we show the three configurations h (6) 1a (1a| 0 ), h (6) 2b (2b| 0 ), h (6) 3c (3c| 0 ), respectively. In these cases, at each Wyckoff position, there exists a Wannier orbital with C n eigenvalue equal to one, i.e. s-orbital. In this calculation, we need to put electrons into the system so as to respect charge neutrality of the system. Therefore, we put ν electrons per unit cell, where ν is the number of occupied bands. In other words,  Fig. S5(b)), all the Wannier orbitals are inside the crystal, and we can safely put all the electrons into all the Wyckoff positions b. On the other hand, in h (6) 1a (1a| 0 ) (Supplementary Fig. S5(a)) and in h (6) 3c (3c| 0 ) (Supplementary Fig. S5(c)) , some Wannier orbitals are at the boundary of the crystal, and such orbitals are no longer identical with those in the bulk. Therefore, we need to specify how we put electrons into these Wannier orbitals at the crystal boundary. In fact, as discussed in Ref. S5, the bulk Wannier orbitals can be considered as bonding states composed of states at several sites inside of each block slightly away from the Wyckoff position considered. In such a way we can "divide" the Wannier orbitals into the sites inside each block, and one can safely define a lattice model in a finite-sized crystal. In the present case, the s-like orbitals are totally symmetric bonding states formed by the states at these sites close to one of the Wyckoff positions, and they are fully occupied in the bulk, leading to a finite gap, while the edge and corner sites (coming from the "division" of the Wannier orbitals) might support partially occupied states at or near the Fermi energy.
From this argument, following Ref. S5, we first classify the units into a bulk area, edge areas, and corner areas. We put the electrons into the crystal so as to respect charge neutrality. We fill the Wannier orbitals inside of the crystal first. After filling up all the Wannier orbitals inside the crystal, there remain some electrons left unallocated. At this point, the bulk area becomes charge neutral, but the edge and corner areas might still be lacking some electrons and not charge neutral. In particular, the edge areas might already be partially occupied by bulk electrons, but its filling can be still less than that of the bulk area. We use the unallocated electrons to fill up all the edge areas up to the filling same as the bulk, i.e. so that the charge neutrality holds at the edge units. We then de-fine an edge filling ν edge (mod 1) to be the number of these additional electrons put onto the edge areas per unit cell. At this point, there may still remain some unallocated electrons, and let I be the number of the remaining electrons. Then at each corner, the corner is positively charged with the charge Q corner = I|e|/n, if such electrons are put away. In other words, if we put these electrons into the n corners, the corner states at each corner are filled by I/n = Q corner /|e| electrons. Here, we adopt a convention of an electron charge being e = −|e|(< 0), similarly to Ref. S5.
In Supplementary Figs. S5(a), (b) and (c), in each triangular unit, we write the number of bulk electrons, occupying the bulk Wannier orbitals, and a way of classifying bulk, edge and corner areas. In Supplementary Fig. S5(b), the edge filling and corner charge are both zero from the previous argument. On the other hand, in Supplementary Fig. S5(a), the bulk filling is ν = 1, and the edge area accommodates 1/2 bulk electron per unit cell. Therefore, we shall put extra 1 2 electron per unit cell at the edge, and we get ν edge = 1 2 . In Supplementary  Fig. S5(c), we also get ν edge = 1 2 . In both cases, the edge is metallic. In Supplementary Fig. S5(a) the number of remaining unallocated electrons is two, and the corner charge is Q corner = 2|e|/6 = |e|/3, and in Supplementary Fig. S5(c) its number is three, and we get is Q corner = 3|e|/6 = |e|/2.
We note that there is some flexibility in how to separate the corner area and the edge area. Instead of the choices of the edge and corner areas in Supplementary Fig. S5(a), (b) and (c), we can take alternative choices shown in Supplementary  Fig. S5(d), (e) and (f), respectively. This changes the amount of the corner charge Q corner by an amount of edge filling times charge |e|, which is |e| 2 in h 1a and h 3c but zero in h 2b . Thus, the corner charge has such kind of ambiguity in general, but this ambiguity disappears when we take some linear combitations of these irreducible representations to make the edge insulaing.
We also note that in a macroscopic system, when the edge is partially filled, i.e. metallic, a slight change in the edge will give rise to a huge change in a corner charge. For example, by changing the details of the on-site energy and hopping on the edges and at the corners, the corner charge largely changes. Thus the corner charge is not determined uniquely, independently from the details of boundary conditions at the edges and corners. On the other hand, in linear combinations of the Wannier configurations with integer filling of the edges, the edge spectrum has a gap, and the corner charge does not change under a small change of the Hamiltonian. Therefore the corner charge is well defined. Thus the corner charge calculated above does not have a physical meaning when the edge is metallic, but they restore their meaning when we make linear combinations of the Wannier orbitals so that the edge becomes insulating.
Next, we build a formula for the fractional corner charge, in a similar manner as in Ref. S5. First, we consider a Hamiltonian h expressed as a linear combination of Wannier configurations: h 1a ≡ h (6) 1a  (6) h 2b (6) h 3c (6) h 1a (6) h 2b (6) h 3c 3c (3c| 0 ), and h 3c ≡ h (6) 3c (3c| π ). Thus, we write where j = 1a, 1a , 2a, 2a , 2b, 4b, 3c, 3c and Z ≥ is a set of nonnegative integers. As discussed in Ref. S5, the topological invariants of h are also represented as a linear combination of the topologiclal invariants for the constituent Wannier configurations. Thus, a set of topological invariants χ can be written as where χ ( j) is a set of topological invariants for the Wannier configuration h j . As we explain later, when we consider a triangular block, we need to add the filling ν, i.e. the number of occupied bands, to the list of topological invariants in χ, as is different from the previous study. In our case , ν), and the values for topological invariants χ for the eight Wannier configurations are summarized in Table S1. In this case, the corner charge of the Hamiltonian h is represented as the summation of the corner charge Q j for the Wannier configuration h j : Derivation of the generators Among the eight Wannier configurations listed in Table S1, one can pick up the three configurations h 1a , h 2b , and h 3c as generators, and we need to consider only these three generators for calculation of the corner charge. Any insulators can be deformed to a linear combinations of the eight Wannier configurations with its coefficients α = t symmetry WC invariants P Q corner ν edge [M (2) 1 ] [K (3) 1 ] ν C 6 (hexagonal) h (6) 1a (1a| 0 ) 0 0 1 (0, 0) 0 0 (α 1a , α 1a , α 2a , α 2a , α 2b , α 4b , α 3c , α 3c ). We get By definition, this map C is a linear map from the set A = Z 8 ≥ of the vector α to the set C = Z 2 × Z ≥ of the topological invariants χ. We extend these sets A and C toÃ = Z 8 andC = Z 3 , respectively, in order to make these sets Abelian groups. We also extend the domain and the range of the map C toÃ andC, respectively, and call the new map asC. Obviously, KerC 0, and this means that there are multiple α's giving the same χ. By identifying all the values of α giving the same χ, we can reduce the number of α to be considered. Then, we getÃ This means that we can choose three configurations as generators ofÃ/KerC. The linear combinations of the generators α| generator have one-to-one correspondence with χ ∈C . If we define the restriction ofC to the generators asC| generator , and its matrix representation as C χ | generator , the one-to-one correspondence is explicitly written as follows: On the other hand, the map B defined by (S4) maps α inÃ to Q corner (modulo 1) in B = R/Z. Here the vector Q = (Q j ) in Eq. (S4) is given by (0, 0, 0, 0, 2e 3 , e 3 , e 2 , e 2 ) for hexagonal blocks and ( 2e 3 , 2e 3 , e 3 , e 3 , 0, 0, e 2 , e 2 ) for triangular blocks. By straightforward calculation we can check that through this map, KerC is mapped to Q corner = 0. This means that if two Wannier configurations α 1 , α 2 have same χ, i.e., (α 1 − α 2 ) ∈ KerC, they have the same Q corner . Therefore, the configurations α| generator is sufficient for calculation of the corner charge Q corner , and we can calculate Q corner directly from the topological invariants χ as follows: In the following we choose h (6) 1a (1a| 0 ), h (6) 2b (2b| 0 ), h (6) 3c (3c| 0 ), because they are represented as Wannier configurations with sorbitals at each Wyckoff position. We note that the generators taken in Ref. S5 are h (6) 1a ≡ h (6) 1a (1a| π ), h (6) 4b ≡ h (6) 4b (2b| ±2π/3 ), and h (6) 3c ≡ h (6) 3c (3c| π ), and are different from ours. Derivation of the formula of the corner charge from generators Thus from these arguments, we can safely restrict ourselves to the generators h 1a = h (6) 1a (1a| 0 ), h 2b = h (6) 2b (2b| 0 ), and h 3c = h (6) 3c (3c| 0 ). Equation (S3) is explicitly calculated from the Supplementary Table S1: Therefore, we can calculate the coefficients α j from the topological invariants χ as follows: (S10) In the previous study [S5] of hexagonal blocks, we have (Q H 1a , Q H 2b , Q H 3c ) = (0, |e| 3 , |e| 2 ). But in our case of triangular blocks, (Q T 1a , Q T 2b , Q T 3c ) = ( |e| 3 , 0, |e| 2 ) (Supplementary Figs. S5(a),(b),(c)) , as shown in Table S1. Henceforth, the superscripts H and T represent hexagonal and triangular blocks, respectively. The difference in the corner charge comes from the different choice of the fundamental units, which results in the different termination of the system. Finally, by combining Eq.(S4) and Eq. (S10), we get a formula for the corner charge for triangular blocks:  edge is insulating, having an integer filling ν edge . From the above formula for the triangular blocks, we get Because [M (2) 1 ] and [K (3) 1 ] are always even integers for any Wannier configurations, as seen from ImC =C = 2Z×2Z×Z , we conclude that only when ν is even, the corner charge for the triangular blocks is well defined.
It is noted that in the previous study with the hexagonal block [S5], one does not need to include ν into the formula Eq. (S12) for the corner charge, while in the present theory with the triangular block, ν enters the formula Eq. (S11) of the corner charge. It is understood if we consider the Wyckoff position a. The Wyckoff position a has the full C 6 symmetry, meaning that the irreducible representations at all the k points are the same [M (2) 1 ] = 0 = [K (3) 1 ]. For the hexagonal blocks, it has neither an edge charge nor a corner charge, and we do not have to take into account this case with Wannier centers at the Wyckoff position a, whereas for the triangular blocks, it has an edge charge and the corner charge, and we need to take this configuration into account. The additional term proportional to ν is needed to count the number of bands which have Wannier centers at the Wyckoff position a.

Supplementary Material 6. MODEL CALCULATION OF THE ELECTRONIC BAND STRUCTURE
We check the TCI phase in the apatite by a tight-binding model calculation with hydrogen atoms put at the locations corresponding to the Wannier centers in the apatite as shown in Supplementary Figs. S6(a)-(d). The bulk band structure is shown in Supplementary Fig. S6(e). Two Wannier functions for the occupied bonding orbitals are located at the C 6 centers. Supplementary Figs. S6(f),(g) are the corresponding band structures for the surface and the hinge, respectively, which well describes their characteristic in the apatite. Indeed the hinge is 2/3-filled, in accordance with the ab initio calculation.
Here, band structure for a model corresponding to the apatite with Coulomb interaction is calculated by the firstprinciples calculation. The reason for using the ab initio calculation is to introduce charge neutral condition and Coulomb interaction at the surface and hinge. Supplementary  Figure S7(c) shows the band structure of our interacting model corresponding to the apatite. The occupation number as a spinless system is 1/3 for one site. Calculation is performed with an atomic number of Z = 2/3 using a virtual crystal approximation. The gap is formed primarily by the hopping in the ab plane, and the splittings inside the valence and the conduction bands are due to the hopping for the caxis. Other hoppings except for the hopping shown in the Supplementary Figures S6(b)-(d) are 1 meV or less and can be ignored. The band structure for a slab with (0110) surfaces is also insulating as shown in Supplementary Fig. S7(d).
Such an insulating phase is consistent with the Zak phase 0 (mod 2π). The valence band and the conduction band near the Fermi level originate from the surface orbitals (Supplementary Fig. S7(e)). The band structure of the hinge is shown in Supplementary Fig. S7(f). Unlike the insulating bulk and surface, the topologically protected 2/3-filled states appear on the hinge (Supplementary Figs. S7(g),(h)).
In the tight-binding model without the interaction, we can-not uniquely determine the position of the Fermi level. The Fermi level must exist in the energy gap of the bulk and surface band structures but the position of the Fermi level is determined by the long-range Coulomb interaction. In the main text of this letter, we choose the Fermi level as the bottom of the bulk and surface bands.