{'fraction_non_alphanumeric': 0.043074075632019516, 'fraction_numerical': 0.014785681235014511, 'mean_word_length': 4.168152173913043, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
General corner charge formula in two-dimensional C_n-symmetric higher-order topological insulators

In this paper, we derive a general formula for the quantized fractional corner charge in two-dimensional C_n-symmetric higher-order topological insulators. We assume that the electronic states can be described by the Wannier functions and that the edges are charge neutral, but we do not assume vanishing bulk electric polarization. We expand the scope of the corner charge formula obtained in previous works by considering more general surface conditions, such as surfaces with higher Miller index and surfaces with surface reconstruction. Our theory is applicable even when the electronic states are largely modulated near system boundaries. It also applies to insulators with non-vanishing bulk polarization, and we find that in such cases the value of the corner charge depends on the surface termination even for the same bulk crystal with C_3 or C_4 symmetry, via a difference in the Wyckoff position of the center of the C_n-symmetric crystal.

In this paper, we derive a general formula for the quantized fractional corner charge in twodimensional Cn-symmetric higher-order topological insulators. We assume that the electronic states can be described by the Wannier functions and that the edges are charge neutral, but we do not assume vanishing bulk electric polarization. We expand the scope of the corner charge formula obtained in previous works by considering more general surface conditions, such as surfaces with higher Miller index and surfaces with surface reconstruction. Our theory is applicable even when the electronic states are largely modulated near system boundaries. It also applies to insulators with non-vanishing bulk polarization, and we find that in such cases the value of the corner charge depends on the surface termination even for the same bulk crystal with C3 or C4 symmetry, via a difference in the Wyckoff position of the center of the Cn-symmetric crystal.
As discussed in the previous works, the quantized fractional corner charge in C n -symmetric crystalline insulators are associated with rotational eigenvalues of the bulk wavefunctions 32,36,40,42 ; nonetheless, they are limited to special cases. For example, in C 6 -symmetric systems, Ref. [32] considers the case with hexagonal unit cells that cover the whole crystal, while in Ref. [40] systems covered by triangular building blocks are considered, and so the formulas of the corner charge for these two cases are different. In Ref. [41], although a general formula of the quantized corner charge is discussed, the authors assume that the surface is flat and the electronic charge distribution is the same between the finite C 6 -symmetric system and an infinite system. Most importantly, the formula in Ref. [41] treats the contributions of electrons and ions to corner charges on an equal footing.
Here, a question arises whether the formula for the corner charge in previous works 32,36,40,42 holds for any insulating C n symmetric systems; this question is important for application to real materials, which may have deformations of nuclei positions and of electronic states near the boundaries. In the previous works, the cornercharge formula was derived for systems in which nuclei positions and electronic states are perfectly periodic up to the boundaries, and this result was complemented by an argument that surface decoration will leave the corner charge unaffected. Nevertheless, the surface-decoration argument has not been addressed in detail, and it is not clear whether the resulting corner-charge formula holds in general C n -symmetric systems, including those where the nuclei positions and electronic states may be modulated near the boundary.
In the present paper we derive a general corner-charge formula on general grounds with minimal assumptions. Our study shows that even when the positions of the nuclei and the electronic states are largely modulated near the boundaries, the corner charge is quantized to fractional values as long as the system is C n -symmetric and the edge charge density is zero. One of the merits of our method is that it naturally includes surfaces with higher Miller indices and surfaces with surface reconstructions. Our theory also includes insulators with non-vanishing bulk polarization but without edge charges. In particular, in such cases with C 3 or C 4 symmetries, we find that the value of the corner charge depends on the Wyckoff position of the center of the C n -symmetric crystal. It is in contrast with previous works, where the corner-charge formula assumes zero polarization. This paper is organized as follows. In Sec. II, we derive the general corner charge formula and clarified the assumptions required for the quantization of the corner charge. We show that as long as the edge is gapped and charge neutral in C n -symmetric finite system, the corner charge is fractionally quantized modulo |e|/n. We also discuss the relation of the general formula and the bulk rotational eigenvalues. In Sec. III, we consider systems with nonzero electric polarization, which is beyond the scope of the previous corner charge formula. In Sec. IV, we consider more general surface conditions, such as general edge orientations and surface reconstructions. Conclusion and discussion are given in Sec. V.

II. CORNER CHARGE FORMULA
In this section, we construct a general theory to calculate the corner charges in an insulator with C n -rotation symmetry, based on the notion of filling anomaly. We derive its fractional quantization under the following minimal assumptions. We assume that the system is topologically trivial in the sense that it is adiabatically connected to an atomic limit 32,41 . It assures that the system has no topological edge states in the gap.
To explain another important assumption, we show the definition of the filling anomaly adopted in this paper: • In some bulk insulators, charge neutrality is incompatible with existence of a gap at the Fermi energy for the whole system including the boundries; in such cases, by adding or subtracting a few electrons from charge neutrality, the Fermi energy is shifted and the system becomes gapped including the boundaries. This deviation of the electron number from charge neutrality is called filling anomaly in this paper.
Based on this definition of the filling anomaly, we assume that the finite-sized C n -symmetric crystal considered is gapped. These two assumptions, namely topological triviality and existence of the gap, guarantee that the occupied states are described in term of localized orbitals. We note that the gap in this finite crystal should survive in the large system size limit. For example, in a hexagonal-shaped crystal with C 6 symmetry, suppose the edge and bulk are gapped, and the corners support one in-gap bound state per each corner, with the Fermi energy at these in-gap corner states. Then a tiny hybridization gap O(e −L/λ ) appears within the sixfold multiplet of the corner states, where L is the distance between neighboring corners and λ is a penetration depth of the corner states along the edges. In this case, the occupied states at the corners cannot be described as localized orbitals, but are extended over the system size, and it violates our assumptions. Thus, these assumptions justify the description of occupied electronic states as localized orbitals.
In the previous works 32, 36,[40][41][42] , the formula of the quantized corner charge is derived for cases with the electronic charge distribution being perfectly periodic even near the system boundaries. It is then argued that the quantization will remain even with surface decorations onto the system while preserving the C n symmetry and the bulk gap. Nonetheless, details of surface-decoration argument have not been addressed, and it is not clear in which kinds of systems the corner charge remains quantized to be the value determined by bulk electronic states. In real materials, electronic states are determined as solutions of the Schrödinger equation, and nuclei positions are determined as minima of the total energy; this necessarily leads to deviations of electronic wavefunctions and nuclei positions near the boundaries. In many cases, details of the nuclei positions and electronic states near the boundaries are complicated and not easily discussed on general grounds. Thus, a general proof for C n -symmetric insulators is desired. In the following, we give a general proof for this problem.
In this paper, we assume that the charges of the respective ions are integer multiples of the electronic charge |e|. It naturally holds in real materials. On the other hand, some theoretical models may have non-integer ionic charges, but such models are outside of our theory.

A. Derivation of the general corner charge formula
We consider a crystal of finite size which preserves the C n symmetry. We assume that the system is insulating, in the sense that the Fermi level is within the energy gap for both the bulk and the boundaries. In some cases, degenerate corner states exist at the Fermi energy, which violates the above assumption, when charge neutrality is assumed. However, we can add or remove several electrons so that no states are at the Fermi energy. This deviation from charge neutrality is called filling anomaly 32 .
In order to calculate the filling anomaly, we divide the system conceptually into two C n -symmetric regions (I) and (II), as shown in Fig. 1. The regions (I) and (II) can be regarded as surface and bulk regions, respectively, and we explain this division in the following. The width of the surface region (I), d edge , is taken to be large enough so that both the electronic states and the ionic positions in the region (II) are the same as those at the bulk.
To calculate the filling anomaly, we proceed with the following steps. First, we classify the ions into two groups according to which region they belong to, so that no ions are at the junctures of two regions.
Second, we classify the electronic states. We assume that the occupied states of this system can be described by a C n -symmetric basis set {φ i } composed of exponentially localized basis orbitals 43 . Well within the crystal, the electronic states can be well described by bulk Wannier functions W nl , and we take these Wannier functions W nl to be among {φ i }. The remaining φ i 's are additional localized basis orbitals ϕ (add) i near the boundary. We classify them into regions (I) and (II) so that the set of orbitals in region (I) and that in region (II) are C n symmetric. By construction the numbers of ions and electronic localized basis orbitals in the region (I) are integer multiples of n. Here, we impose a condition that all the additional localized basis orbitals ϕ (add) i are classified into the region (I), and we choose the region (I) large enough so that this condition is satisfied. Thus, even if the electronic states near the boundary are largely deformed, our theory is still applicable, by incorporating these deformed wavefunctions into region (I). We note that this construction of localized basis functions is similar to that used in the argument in Ref. [43] to show that the bulk polarization in terms of the Zak phase gives surface charge density. We note that some degrees of freedom remain in this division into two regions, but they do not affect the following argument.
Third, we calculate the total charge of the finite crystal in terms of mod n|e|. First, the total charge contribution from the region (I) should be 0 (mod n|e|) because the number of the nuclei and the electronic basis orbitals within region (I) is an integer multiple of n. Next, we calculate total charges of electrons and ions in the region (II). First, in the region (II), all the orbitals are Wannier orbitals W nl , and therefore positions of the ions and electrons are classified in terms of Wyckoff positions. Because the whole crystal is C n symmetric, the center of the crystal should be located at the Wyckoff position with C n local symmetry. Let 1a be the Wyckoff position with C n symmetry located at the center of the crystal. Then, the ions not located at 1a contribute nM |e| (M ∈ Z) to the net charge in the region (II), while those at 1a contribute nM |e| + n 1a |e| is the charge of the ion at the 1a position. Therefore, the total ionic charge in the region (II) is equal to n (ion) 1a |e| mod n|e| and so is that in the whole crystal. Similarly, the total electron charge is equal to −n 1a is the number of electron Wannier orbitals at the 1a position. From the above discussion, the total charge of the system in terms of modulo n|e| is calculated as follows, where n (ion) 1a |e| is an ionic charge at 1a, and n (e) 1a is the number of electronic Wannier functions at 1a. This formula represents a filling anomaly in the crystal.
Fourth, we discuss how Q tot is distributed in the crystal. As a first step, we define the moving average of the total charge density ρ(x) as follows: Here, ρ(x) is the total charge density including ions and electrons, Ω is a C n -symmetric unit cell whose center is at the origin and |Ω| is the area of Ω. From the definition, ρ(x) have the following properties: FIG. 2. (Color online) A conceptual picture of the region division for calculating the corner charge for a C4-symmetric system as an example. The total charge in each region is calculated by integrating the moving average of the total charge density ρ(x) in each region. In the bulk region (C), ρ(x) = 0, and the total charge in the bulk region Q bulk is zero. The edge region (B) includes an integer number of edge unit cells. The charge of each edge unit cell is calculated by the product of the edge charge density σ = P · n and the edge period a edge . The total charge is represented as summation of the edge charges and the corner charges: Qtot = n(Q edge + Qcorner). Therefore, when Q edge vanishes, Qcorner ≡ Qtot/n (mod e).
We note that as long as ρ(x) satisfies Eqs. (3), (4) and (5), we can also use other definitions for moving averages. In a sufficiently large system, we expect that electronic states in the middle part of the edges far away from the corners are almost the same as those on an edge of a semiinfinite system with an open boundary condition in one direction. To calculate the corner charge based on this expectation, we divide the system into three types of regions: the (A) corner, (B) edge, and (C) bulk regions in a C n -symmetric way as schematically shown in Fig. 2. The edge region (B) should be chosen to satisfy the following conditions. First, it should be sufficiently far away from the corners, such that the electronic states in the edge region (B) are the same as those on an edge of the corresponding semi-infinite system. Second, it should have a width equal to an integer multiple of the period along the edge, N edge a edge , and a sufficiently long depth into the bulk, d edge , such that in the remaining bulk region (C), the electronic states are the same as those for a bulk infinite system. Here, we assume that one edge region includes N edge edge unit cells. Since ρ(x) takes nonzero values only in the (A) corner and (B) edge regions, its integral can be calculated as the sum of the total charges in these regions: where Q edge and Q corner are the total charges at one of the (B) edge regions and one of the (C) corner regions, respectively. Then Q edge can be calculated as: Here, N edge is the number of edge unit cells within one edge region, σ edge is the edge charge density, and a edge is the period of the edge unit cell. Finally, we calculate the corner charge. We assume that the edge is charge neutral, σ edge = 0. Then the extra charge Q tot is distributed into the n corners. Since the crystal is assumed to be C n -symmetric, the distributed charges are equal among the corners, and this gives a corner charge formula: Here, Q (n) c 1a means the corner charge when the center of the crystal is located at the 1a Wyckoff position.
If the edge charge density σ edge = 0, the extra charge Q tot is distributed not only to the corners, but also to the edges. Since the separation of charges between the corners and the edges depends on the choice of the corner regions, the corner charge is not well defined.

B. Comparison with previous works
Our discussion so far clarifies that the formula (8) for the corner charge applies to a general system, as long as the system is C n -symmetric and is insulating over the whole system, irrespective of the system details. In particular, our theory applies to (i) systems with boundaries with high Miller indices, (ii) those with surface reconstructions, and (iii) systems with nonzero P but with zero edge charge density σ edge = 0. In the following, we briefly explain these three cases. Details of (i) and (ii) are given in Sec. IV and those of (iii) in Sec. III.
First is the case (i) with boundaries with higher Miller indices. In the previous works, the arguments are limited to boundaries with low Miller indices such as (10) and (01) edges, but our theory applies to general Miller indices. For example, a C 4 -symmetric square system with edges along (34) and (43) directions is covered by our theory. Second, we can treat cases with (ii) boundaries with surface reconstructions. It is nontrivial because the surface reconstruction may make the perioidicity at edges different from that in the bulk. In some cases, as we discusss later, the edges without the surface reconstruction are gapless and σ edge is zero, but those with the surface reconstruction may become gapped, which makes our theory applicable.
Third, our theory covers also the cases (iii) with P = 0. Our assumption is σ edge = 0, which does not necessarily mean P = 0. In contrast, in the previous works 32, 36,40,41 , the bulk polarization P is assumed to be 0, which ensures that the edge charge density σ edge is 0. We find that this assumption of P = 0 in the previous works is too strong.
The bulk polarization P is related to the edge charge density 43-46 via σ edge ≡ P · n (mod |e|/a edge ), (9) where n is the unit normal vector of the edge. It is not an equality, because the bulk polarization P is defined only in terms of modulo some unit [43][44][45][46] , which we call polarization quantum. Thus, our assumption σ edge = 0 means P · n ≡ 0 (mod |e|/a edge ), which does not necessarily mean P = 0.

C. Corner-charge formula in terms of rotation eigenvalues
In this section, we rewrite our formula for the corner charge in terms of rotation eigenvalues of the Bloch eigenstates at high-symmetry wavevectors following Refs. [32,36,40], in three different systems with and without time-reversal symmetry, and with and without spin-orbit coupling, separately. We need to rewrite it, because eigenstates in a crystal are Bloch wavefunctions, extended over the crystal, and one cannot directly calculate n (e) 1a from the Bloch wavefunctions. Then, we discuss a relation between our formula and the formulas in the previous works 32, 36,40 . The notable difference from the previous works is that the Wyckoff positions of the ions are not always 1a, and the corner charge generally depends on the filling ν as in Ref. [40]. The detail of the calculation is shown in the Appendix.
By following Refs. [32,36,40], we introduce topological invariants which distinguish the Wannier functions with different Wyckoff positions. With a sufficient number of topological invariants, there can be a one-to-one correspondence between the values of the topological invariants and Wyckoff positions of the occupied Wannier orbitals. We assume that the occupied bands can be written in terms of Wannier functions. Then, the number of Wannier orbitals at each Wyckoff position is determined from the topological invariant for the occupied bands. Figure 3 shows the Wyckoff positions under C n rotational symmetry (n = 3, 4, 6). In each system with C n symmetry, there are four kinds of Wyckoff positions. We define the number of occupied Wannier functions at Wyckoff position mX as n For the calculation of the corner charge, the value of n (e) 1a is needed. As we show in the appendix, mn (e) mX modulo n (X = a) is directly calculated from the rotational topological invariants defined from rotation eigenvalues of the Bloch wavefunctions, and n mX . In the appendix, we show the detailed calculation of n (e) mX .

class A
First, we consider a system called class A. Unlike the class AI systems 32,40 or class AII systems 36,42 considered in previous studies, the class A systems do not have timereversal symmetry. From Eq. (8) and Appendix A, the corner charge is calculated as follows: 3 ] (mod e), (12) Here, the rotation topological invariant [Π (n) p ] represents the difference of the number of states with C n eigenvalue e 2π(p−1)i n for spinless systems, and e 2π(p−1)i n e πi n for spinful systems, between the rotation-invariant k-points Π and Γ. We assumed (C n ) n = 1 and (C n ) n = −1 for spinless and spinful systems, respectively.
In addition to Eqs. (11)- (13), where the center of the C n -symmetric crystal is at 1a, in appendix A, we also discuss the formula when the center of the crystal is at other Wyckoff positions, i.e., 1b and 1c for C 3 and 1b for C 4 . The formulae are summarized as follows: 3 ] (mod e). (16)

class AI
Here, we consider a spinless system with time-reversal symmetry called class AI. In this case, there are some constraints on the topological invariants [Π (n) p ], and the corner charge formula is somewhat simplified.
In particular, these formulae are reproduced to the results in Ref. [32] only if all the ions are located at the center of the C n -symmetric unit cell (i.e., Wyckoff position 1a) leading to n (ion) 1a = ν and if the electric polarization is zero. The formulae are reduced to Eq. (11) of Ref. [32], i.e., 1 ] (mod e).
As shown in appendix B, we also discuss the formula when the center of the C n -symmetric crystal is not at the Wyckoff position 1a. The formulae are summarized as follows: also the case in compounds containing elements with odd atomic numbers; a change in crystal terminations results in a change of the corner charge by e. For the above reasons, we deal with corner charges defined by modulo e, even in the class AII. As with class AI systems, the topological invariants [Π (n) p ] have some constraints due to time-reversal symmetry, and the corner charge formula is somewhat simplified.
In appendix C, we also discuss the formula when the center of the C n -symmetric crystal is not at the Wyckoff position 1a. The formulae are summarized as follows:

III. CORNER CHARGE WITH NONZERO BULK POLARIZATION
In this section, we consider cases with nonzero bulk polarization P . In the previous works 32, 36,40,41 , the bulk polarization P is set to be 0, to ensure that the edge charge density per edge unit cell is 0 modulo |e|. However, as we explain in the following, even if the bulk polarization P is nonzero, as long as the edge charge density σ edge is 0, our general formula is applicable.
To illustrate this, we consider a C 4 -symmetric tightbinding model with nonzero bulk polarization as an example. We show that the fractional corner charge appears even though the bulk polarization is nonzero, as long as the edge charge density is 0. Moreover, we also show that two finite C 4 -symmetric flakes with different terminations can have different corner charges due to difference in the C 4 centers. This property is specific to cases with nonzero polarization.
Our model is a simple spinless tight-binding model. It has background positive ions with 2|e| charge at 1a, being the center of the square unit cell (see Fig. 4(a),(b)). The Hamiltonian is given asĤ Fig. 4(a). We set t s = 2 and t w = 0.8 in the numerical calculation. The bulk bands are gapped, and we set the Fermi energy E F = 0 within the gap. The parameters are chosen so that the occupied Wannier orbitals are located at 2c Wyckoff positions, in the middle of the red bonds. In this example, the bulk polarization is P = (e/2a)(1, 1) (mod e/a), where a is the lattice constant. Therefore, if we cut the system along the x and y axis as shown in Fig. 4(c), nonzero fractional edge charge appears, σ edge a edge = P · (1, 0)a = e 2 , where the size of the edge unit cell a edge is taken to be equal to a. This model has half-filled edge states within the bulk gap, and they are metallic (see Fig. 4(d)). Thus, this is out of the scope of the corner charge formula. On the other hand, if the system is cut along the (11) direction as shown in Fig. 4(e), the edge charge is 0 modulo e: Then the edge is insulating as shown in Fig. 4(f), and we can apply our general formula.
In C 4 -symmetric systems, there are two C 4 -symmetric Wyckoff positions, 1a and 1b. As shown in Fig. 3(b), their positions are r 1a = 0 and r 1b = 1 2 a 1 + 1 2 a 2 , respectively. Therefore, which formula to be used depends on whether the C 4 center of the finite crystal is 1a or 1b. If the center of the finite crystal is at 1b, we need to use a formula with 1a replaced by 1b in Eq. (8): The difference between Eqs. (8) and (32) is expressed as Here, q X = (n mX )|e| (mX = 1a, 1b, 2c, 4d). The r.h.s of Eq. (33) is equal to P ·(a, 0) and P ·(0, a). Therefore, when the quantized polarization along (10) or (01) direction is nonzero, the corner charge is different between a C 4 -symmetric crystal with its center at 1a and one at 1b.
We note that in Ref. [41], systems without polarization are studied, and is shown, meaning that the corner charge is independent of whether the crystal center is at 1a or 1b. It is a special case of our result. This difference between Q c 1a and Q c 1b can be attributed to a difference in the shape of the corners depending on the C 4 center. For example, when the shape of the crystal is as shown in Fig. 5(a), the C 4 center is at 1a, and the corner charge is calculated as Q 1a )|e|/4 = |e|/2. On the other hand, when the crystal shape is as shown in Fig. 5(c), the C 4 center is at 1b, and the corner charge is calculated as Q 1b )|e|/4 = 0. Correspondingly, in the system shown in Fig. 5(a), where the C 4 center is at 1a, four in-gap corner states appear (see Fig. 5(b)). Here, we briefly explain the relation of the in-gap corner states and the corner charge. In the case of a charge-neutral filling indicated by the line (i) in Fig. 5(b), the corner states are at the Fermi energy. In this case, we cannot define its corner charge. By adding or extracting two electrons as shown by lines (ii) and (iii), the whole system becomes gapped and the total charge deviates from charge neutrality by ±2|e|. As discussed in Sec. II, the nonzero total charge is equally distributed to the four corners, resulting in a corner charge |e|/2 (mod |e|). In contrast, in the system in Fig. 5(c), where the C 4 center is at 1b, there are no corner states (see Fig. 5(d)). This agrees with our result in Eq. (33).
In C 3 -symmetric systems, there are three C 3symmetric Wyckoff positions, 1a, 1b and 1c. As shown in Fig. 3(a), their positions are r 1a = 0, r 1b = 1 3 a 1 + 1 3 a 2 and r 1c = 2 3 a 1 + 2 3 a 2 , respectively. Therefore, similar to the case of C 4 -symmetric systems, the formula for the corner charge depends on whether the C 3 center of the finite crystal is 1a, 1b or 1c, and let Q represent the corresponding corner charges respectively. Then, we obtain As discussed in Ref. [41], the r.h.s of Eq. (35) is equal to both the electric polarization |Ω|P · b1 2π and |Ω|P · b2 2π , where b i (i = 1, 2) is the reciprocal lattice vector defined as a i · b j = 2πδ ij . Therefore, when the quantized polarization along b 1 or b 2 is nonzero, the corner charge depends on the C 3 -center of the finite crystal. In particular, when the polarization is zero, qa 3 = q b 3 = qc 3 (mod e) holds as shown in Ref. [41], and the corner charge is independent of the C 3 -center.
In a C 6 -symmetric system, there is only one C 6symmetric Wyckoff position, 1a. In this case, the corner charge formula is unique.
We discuss the role of a choice of the unit cell in our theory. In the definition of the bulk polarization according to the modern theory of polarization 43 , a choice of the unit cell plays a crucial role, and so does it in the present theory of the corner charge. It is noted that a choice of the unit cell is a gauge degree of freedom, which does not affect the system itself. The edge charge density σ edge and the corner charge Q c are observables, and therefore they cannot depend on a choice of the unit cell.
We stress that in our theory, the unit cell is always the primitive one in the bulk. Meanwhile, one can optionally take a unit cell larger than the primitive one. It has been discussed previously in the context of reducing the cases with nontrivial polarization into those with P ≡ 0, in order to apply the previous theories based on P = 0 36,41 .
Nonetheless, we note that taking a larger unit cell than the primitive one is misleading in some cases. For example, in the C 4 -symmetric model in Fig. 4(a), the polarization P is equal to P = e 2a (1, 1) (mod e/a). Meanwhile, if we take a larger unit cell of 2a × 2a instead of the primitive one, a × a, P is treated in terms of modulo e/(2a), and P ≡ 0 (mod e/(2a)). Thus, with an enlarged unit cell, one cannot distinguish between the two possible cases for C 4 -symmetric cases, i.e. P = e 2a (1, 1) (mod e/a) and P = 0 (mod e/a). In terms of the edge charge density, the polarization quantum corresponds to the charge |e| per one period along the edge. Thus in the present model, when n = (1, 0), the edge charge per one edge periodicity is σ edge a edge ≡ P · na edge = e 2 (mod e) when the edge periodicity taken to be a edge = a. Meanwhile by enlarging the unit cell to 2a×2a, the edge charge density becomes σ edge a edge = e, where a edge = 2a, and it is congruent to zero modulo e, corresponding to the trivial bulk polarization. Here we note that the edge charge density σ edge is a physical observable, whereas the bulk polarization P is not an observable in a strict sense, due to its ambiguity modulo a polarization quantum. Thus it is not a contradiction that the bulk polarization P changes from nontrivial to trivial only by a change of the unit cell, being a gauge transformation.
Our theory, where the unit cell is always taken to be the primitive cell, is simple and convenient, because we do not need to switch the choice of the unit cell in the middle of the calculation. It is natural to work always on the primitive unit cell, because the choice of the unit cell, being a gauge degree of freedom, should not affect physical observables. One needs to be careful in dealing with the bulk polarization P , because it is not an observable in a strict sense as we mentioned earlier. There might be a belief that the corner-charge quantization requires P = 0, as has been adopted in previous papers, but our theory has revealed that the correct condition necessary for the corner-charge quantization is σ edge = 0 and not P = 0.

IV. GENERAL EDGE ORIENTATIONS AND SURFACE RECONSTRUCTIONS
In this section, we consider edge charge densities for edges with general orientations. When the edge charge density is nonzero, the corner charge is nonzero. Nonetheless, even in these cases, surface reconstructions may make the edge density to vanish, and the corner charge becomes quantized. In this section, we discuss this in general systems.

A. surface reconstruction
Here we consider systems with surface reconstruction. Conceptual pictures of the surface reconstruction is shown in Fig. 6. We assume that the system has nonzero (e,f) The ribbon system cut along (11) direction and its energy spectrum. Unlike that along (01) direction, the edge energy spectrum is gapped. The finite system with the C4 center at 1a and its energy spectrum. In the numerical calculation, the system size is N = 8, and it contains 113 unit cells in total. The charge-neutral filling shown by line (i) in the figure makes the system metallic. In an insulating filling as shown by lines (ii) and (iii), the total charge deviates from charge neutrality by ±2|e|. As discussed in the main text, the nonzero total charge is equally distributed to the four corners, resulting in a corner charge |e|/2 (mod |e|). (c,d) The finite system with the C4 center at 1b and its energy spectrum. In the numerical calculation, the system size is N = 8, and it contains 112 unit cells in total. quantized electric polarization in the bulk. Then the edge have nonzero fractional edge charge σ edge . In Fig. 6(a) we show its band structure when the edge states are metallic. However, it can under go Peierls transition or surface reconstruction to multiply edge periodicity and open a gap. In Fig. 6(b), due to the surface reconstruction, the surface periodicity is doubled, and the edge charge becomes an integer in the enlarged surface unit cells. Then, the edge energy spectrum can be gapped and insulating. In this case, one can apply our general formula.

B. General edge orientations
Here, we consider the surface with a high Miller index. First, we consider the C 4 -symmetric system with basis vectors of a 1 = (a, 0) and a 2 = (0, a), as shown in Fig. 3(b). In the C 4 -symmetric system, the bulk electric polarization is quantized to be P ≡ (0, 0) or ( e 2a , e 2a ) (mod e a (m 1 , m 2 )), (m 1 , m 2 ∈ Z) 32,41,43 . In the following, we consider the case with P ≡ ( e 2a , e 2a ). Next, for the surface with a Miller index (α, β) (α, β ∈ Z, and α and β are coprime), the surface normal vector n is defined as follows: Therefore, the edge charge in one edge unit cell is calculated as follows: Here, we set the size of the edge unit cell to be the minimal one, a edge = |βa 1 − αa 2 | = α 2 + β 2 a, as shown in Fig. 7. If both α and β are odd integers, σa edge ≡ 0 (mod e), and otherwise σa edge ≡ e/2 (mod e). Obviously, the edge can be charge neutral when σa edge ≡ 0 (mod e).
Even when σa edge ≡ e/2 (mod e), after surface reconstruction with doubling the surface period a edge = 2a edge , we get σa edge ≡ 0 (mod e), and the edge can be charge neutral. In particular, when the system has edge states in the gap, and the fractional edge charge is accommodated in the edge states, the edge states are half-filled and it always has Peierls instability, meaning that the edge becomes insulating via spontaneous lattice deformation with doubling the edge unit cell at sufficiently low temperature. Next, we consider the C 3 -symmetric system with basis vectors of a 1 = ( a 2 , √ 3 2 a) and a 2 = (− a 2 , √ 3 2 a), as shown in Fig. 3(a). In the C 3 -symmetric system, the bulk electric polarization is P = 0, e 3Ω (a 1 + a 2 ) or −e 3Ω (a 1 + a 2 ) (mod e Ω R) 32,41,43 , where Ω = √ 3 2 a 2 is the area of the bulk unit cell and R = m 1 a 1 + m 2 a 2 (m 1 , m 2 ∈ Z) is a lattice vector. In the following, we consider the case with P = ±e 3Ω (a 1 + a 2 ) = (0, ±2e 3a ). The surface normal vector n for the surface with Miller index (α, β) (α, β ∈ Z, and α and β are coprime) is defined as follows: Therefore, the edge charge in one edge unit cell is calculated as follows: Here, we used a edge = |βa 1 − αa 2 | = α 2 − αβ + β 2 a. If α + β ∈ 3Z, σa edge ≡ 0 (mod e), and otherwise σa edge ≡ e/3 or (2e)/3 (mod e). Obviously, edge can be charge neutral when σa edge ≡ 0 (mod e). Even when σa edge ≡ e/3 or (2e)/3 (mod e), after surface reconstruction with tripling the surface period a edge = 3a edge , we get σa edge ≡ 0 (mod e), and the edge can be charge neutral. Finally, we consider the C 6 -symmetric system with basis vectors of a 1 = (a, 0) and a 2 = (− a 2 , √ 3 2 a), as shown in Fig. 3(c). In the C 6 -symmetric system, the bulk electric polarization is always P = 0 (mod e Ω R) 32,41,43 , where Ω = √ 3 2 a 2 is the area of the bulk unit cell and R = (m 1 a − m2 2 a, √ 3 2 m 2 a) (m 1 , m 2 ∈ Z) is a lattice vector. Therefore, σ edge = P · n is 0 for surfaces in any direction (mod e/a edge ).
To summarize, in either case, even when P = 0 in the C n -symmetric systems, σa edge ≡ 0 (mod e) can be achieved with proper surface reconstruction thanks to the fractional quantization of P .

A. Tight-binding models
Here, we will explain how to apply our formula to tightbinding models. In principle, tight-binding models contain only electron degrees of freedom. In order to calculate the corner charge, it is necessary to give appropriate information on the positions of ions.
First, it should be pointed out that the sites of the tight-binding models are different from the positions of ions. The apatite electride is one of the examples showing that these two are independent 40 . The electronic states near the Fermi level of the apatite electride are well described by the Wannier functions localized at the one-dimensional hollows in the crystal. Therefore the tight-binding model made from these Wannier orbitals has its sites inside the hollows. Since there are no ions in the hollows, the ion positions and the sites are clearly different in this model.
Next, we apply our formula to the well-known 2D Benalcazar-Bernevig-Hughes (BBH) model 24,26 as an example. This model is a spinless tight-binding model on a square lattice with π-flux, showing a higher-order topological insulating phase. As in Ref. [26], we assume that all the positive ions are located at the center of the C 4symmetric unit cell, i.e., n (ion) 1a = ν holds. The k-space Hamiltonian is given as follows: Here, γ and λ are real parameters, and τ j , σ j (j = x, y, z) are Pauli matrices. For simplicity, we take the lattice constants as a = 1. The Hamiltonian has C 4 -symmetrŷ Here, due to the π-flux,Ĉ 4 satisfies (Ĉ 4 ) 4 = −1, and its eigenvalues are C 4 = e (2n+1)πi 4 (n = 0, 1, 2, 3). In the derivation of the corner charge formula in this paper, we assumed that (Ĉ n ) n = 1 for spinless systems. Therefore, our formula cannot be applied as it is, when the rotation eigenvalues are modified due to the π-flux. However, we can formally determine the corner charge by applying the formula for the spinful class A system as if the model described a spinful system because (Ĉ 4 ) 4 = −1.  For simplicity, we consider two limits: λ = 0 and γ = 0. When λ = 0, the Hamiltonian becomes H(k) = γ(τ x σ 0 − τ y σ y ), which does not depend on k. Then topological invariants are trivial: [X Topological invariants can be calculated by taking the simultaneous eigenstates of the Hamiltonian and the rotation operationĈ 4 orĈ 2 = (Ĉ 4 ) 2 at each of the high symmetry points Γ = (0, 0), X = (π, 0), and M = (π, π).
The results are summarized in Table. I, and the corner charge is calculated from the Eq. (12): This result agrees with the results of previous studies 24, 26 . This value of the corner charge is stable as long as |γ| < |λ|, where the bulk gap is open.

B. First-principles calculations
In this subsection, we explain ab initio calculations for two materials, to show agreement with our corner-charge formula.
The first example is graphdiyne, which was studied in Ref. [47] as a C 6 -symmetric higher-order topological insulator with a fractional corner charge. Figure 8(a) shows the crystal structure of graphdiyne, which belongs to the wallpaper group p6mm with C 6 rotation symmetry. Since the spin-orbit coupling does not lead to additional band inversion in graphdiyne and therefore does not influence the corner charge, we calculate the spinless electronic band structure of graphdiyne, which is shown in Fig. 8(b), after a full relaxation on atomic positions until the forces on the ions are less than 0.01 eV/Å via Vienna ab initio simulation package (VASP) [48][49][50] . An 11×11×1 k-mesh is used in the BZ for the self-consistent calculations with a plane-wave energy cutoff of 500 eV.  In a C 6 -preserved nanodisk configuration, the eigenvalues are shown in Fig. 8(c), which have six degenerate states #1162, · · · #1167 around the Fermi energy (marked by red dots). These six states are localized at six corners, as is seen in Fig. 8(d) showing the charge distributions of the #1162 state. Because the number of electrons is 1164, we need to add or subtract three electrons to make these six states fully occupied or empty to make the system gapped. Therefore, the value of the corner charge is ± 3e 6 = e 2 (mod e). This is consistent with the corner charge calculation from our formula in Tab. II.
The other example is C 4 -symmetric two-dimensional BiSb. Figure 9(a) show both the top and side views for the crystal structure of BiSb, which belongs to the wallpaper group p4gm with C 4 rotation symmetry. Figure 9(b) is the spinless electronic band structure also calculated by VASP after a full relaxation on atomic positions. An 8 × 8 × 1 k-mesh is used in the BZ for the selfconsistent calculations with a plane-wave energy cutoff of 320 eV. From the corner charge formula calculation in Tab. III, there is no corner charge in a C 4 -preserved nanodisk of BiSb, as shown in Fig. 9(c-d). Namely, the gap is open everywhere including the corners. It means that the filling anomaly is zero. However, the corner state can exist in a C 4 -breaking configuration. Figure 9(e) shows the eigenstates of a C 4 -broken C 2 -preserved nanodisk, which have 8 states around the Fermi energy. Figures 9(fg) are the charge distributions of the #980 and #981 state. The number of electrons is 980. Therefore, corner states appears, but because C 4 -symmetry is not preserved, the corner charge is not quantized. These corner states are not symmetry-protected and their appearance depends on the surface termination. Their appearance is attributed to dangling bonds at the corners.

VI. CONCLUSION AND DISCUSSION
In the present paper, we derived the general formula for the quantized fractional corner charge under the assumptions that the electronic states can be described by Wannier orbitals and the edges are charge neutral. These assumptions are more general than the previous studies, which assumed vanishing bulk electric polarization, and these formulas are more general to be used in different systems with/without time-reversal symmetry/spin-orbit coupling. We expanded the scope of the corner charge formula by considering more general surface conditions, such as surfaces with higher Miller index and surfaces with surface reconstruction.
Our proof also shows that even when the positions of the nuclei and the electronic states are largely modulated near the boundaries, the corner charge is still quantized to fractional values as long as the system is C n -symmetric and the edge charge density is zero. Thus, unlike previous theories, our theory also includes some insulators with non-vanishing bulk polarization, and in such cases with C 3 or C 4 symmetries, we find that the value of the corner charge depends on the Wyckoff position of the center of the C n -symmetric crystal.
We note that the corner charge formula has a simple form like Eq. (8) only when the system has C n rotational symmetry. The corner charge in systems without rotational symmetry is discussed in Refs. [51,52].
Finally, we briefly comment on the generalization of our results to three-dimensional systems. There are two directions for generalization: hinge charge calculation and corner charge calculation. First, let us consider a case where a three-dimensional system can be regarded as a stack of C n -symmetric two-dimensional systems, and cut the system into n-gonal prism. Then, the hinge charge per unit length is equal to the corner charge of the two-dimensional layer divided by the hinge period a hinge . In this case, our formula can be used to determine the hinge charge. However, it is not clear that our formula can be generalized to the calculation of the hinge charge of general three-dimensional systems, and it is future work.
The generalization to the calculation of corner charges is a more non-trivial problem. In order to discuss corner charges in three-dimensional systems, it is necessary to discuss the charge neutrality of the hinge in addition to the charge neutrality of the bulk and the surface. However, it is not clear whether the charge neutrality of the hinge can be determined only from the bulk quantity. Furthermore, three-dimensional systems can be cut out in various ways, and the same crystal may have different corner charges. For example, in the case cubic symmetry, both octahedron and cube maintain the cubic symmetry. These two may have different corner charges due to the different number of corners.
Note added. A recent publication 53 computes the filling anomaly for two-and three-dimensional C 4 symmetric lattices, and our results agree with it for twodimensional C 4 symmetric lattices.
ACKNOWLEDGMENTS R.T. thanks K. Naitou for useful comments. This work was supported by JSPS KAKENHI Grant Numbers JP18J23289, JP18H03678, and JP20H04633.
for spinful systems, between the rotation-invariant k-points Π and Γ. First, we assume that the eigenstates can be described by well-localized Wannier orbitals. Thanks to the discussion based on lattice homotopy, we can limit ourselves to Wannier orbitals localized at one of the high-symmetry Wyckoff positions 54 . Thus, we calculate the topological invariants for such Wannier configurations W j 32,36 , and Tables IV, V and VI summarize the values of the topological number [Π (n) p ] for Wannier configurations for C 3 -, C 4 -and C 6 -symmetric systems, respectively. Different from the Refs. [32,36], we calculate the general corner charge formula including the case where the ion position is not limited to 1a Wyckoff positions. Since the eigenstates are expressed in terms of Wannier orbitals, the topological invariant χ for the occupied bands is represented as the summation of the topological invariants of the Wannier configurations χ (j) : Here, M is the number of Wannier configurations with nonzero topological invariants, C χ = (χ (1) , χ (2) , · · · , χ (M ) ) and α = (α 1 , α 2 , · · · , α M ) T , which represents the numbers of occupied bands corresponding to the Wannier configurations W j . By definition, n (e) mX is represented by the sum of α j belonging to the Wyckoff position mX. Since C χ is generally not a square matrix and is not always invertible, α j itself cannot be uniquely determined from χ. However, in a C N symmetric system, n (e) mX , expressed in terms of α j , can be determined modulo N (N = 3, 4, 6) through direct calculations, which is sufficient for the calculation of the corner charge.
Label Wannier configurations invariants In class AI systems, as compared to class A systems, several Wannier configurations are connected by the time-reversal symmetry Θ (Θ 2 = 1). Then, the Wannier configurations are changed as shown in Tables VII, VIII and IX for C 3 , C 4 and C 6 symmetric systems, respectively.

C4-symmetry
In C 4 -symmetric systems, we need to consider five sets of Wannier configurations in Table VIII. Table VIII From the definition, n = ν and vanishing polarization.

C6-symmetry
In C 6 -symmetric systems, we need to consider four sets of Wannier configurations in Table IX. Table IX is   Here, we briefly explain why [K   In C 3 -symmetric systems, we need to consider four sets of Wannier configurations in Table X. Table X is summarized as follows: From the definition, n 1 ] (mod e).
In C 4 -symmetric systems, we need to consider two sets of Wannier configurations in Table XI. Table XI