Fractional corner charges in spin-orbit coupled crystals

We study two-dimensional spinful insulating phases of matter that are protected by time-reversal and crystalline symmetries. To characterize these phases we employ the concept of corner charge fractionalization: Corners can carry charges that are fractions of even multiples of the electric charge. The charges are quantized and topologically stable as long as all symmetries are preserved. We classify the different corner charge configurations for all point groups, and match them with the corresponding bulk topology. For this we employ symmetry indicators and (nested) Wilson loop invariants. We provide formulas that allow for a convenient calculation of the corner charge from Bloch wavefunctions and illustrate our results using the example of arsenic and antimony monolayers. Depending on the degree of structural buckling, these materials can exhibit two distinct obstructed atomic limits. We present density functional theory calculations for open flakes to support our findings.


I. INTRODUCTION
The classification of insulating phases of matter by topological invariants has been refined in important ways recently. A detailed understanding of symmetries is paramount for the classification of topological phases. An example foundational to the field are the Z 2 topological insulators in two dimensions (2D) and three dimensions (3D) [1][2][3][4] , whose topology is protected by time-reversal symmetry (TRS) and manifests in edge and surface states, respectively, which are immune to Anderson localization. Two characteristics of these phases can equivalently be used to define them as topological: (i) the bulk cannot be adiabatically deformed into the trivial phase (an atomic limit or the vacuum) while retaining the protecting symmetry and (ii) boundary modes protected by the symmetry appear.
Recently, spatial and in particular space group symmetries have been used to define topological properties [5][6][7][8][9][10][11] , characterizing topological crystalline insulators. As this extends the number of known topological phases significantly, it also calls for a sharper definition of what is topological, as the criteria (i) and (ii) do not coincide anymore when spatial symmetries are required for topological protection. For one, the topological bulk-boundary correspondence is extended to include higherorder topological insulators (HOTIs) which exhibit hinge or corner modes in 3D and corner modes in 2D [12][13][14][15][16][17][18][19] . Secondly, the atomic limit as a trivial reference point is not unique 20 . With spatial symmetries, several atomic limits exist that cannot be adiabatically deformed into one another. The physical reference point for an actual material is the atomic limit that corresponds to the physical location of ions. A situation where Wannier charge centers of the occupied bands of an insulator correspond to a different atomic limit is referred to as an obstructed atomic limit 20 (OAL).
This leaves three types of (spatial) symmetry-protected phases that can be distinguished according to the criterion of bulk phase transitions: (1) phases which cannot be adiabatically transformed to any atomic limit, which we refer to as strong topological, (2) phases which correspond to an OAL. Curiously, there is a third category of (3) phases with fragile topology [21][22][23][24][25][26][27][28] , which are not adiabatically deformable to an atomic limit, but can be deformed to one upon the addition of bands that correspond to an atomic limit. Many of the phases of type (2) support point-like boundary states in an open geometry, since the physical boundary of the system may "cut through" the Wannier charge centers in the OAL. At fixed bulk filling, these "dangling" Wannier charge centers become fractionally filled, which allows to define fractional charges, even in the absence of accompanying boundary states [29][30][31][32][33][34] . The notion of a filling anomaly 32 makes this precise: at a filling that corresponds to an insulating band gap with periodic boundary conditions, the system has to be metallic with open boundary conditions when the relevant symmetries are respected, because the boundary modes are then fractionally filled. This notion does not require a spectral symmetry, which is often used to pin boundary modes in the middle of a band gap. A paradigmatic example are one-dimensional (1D) lattices with reflection (inversion) symmetry, as represented by the Su-Schrieffer-Heeger model, for instance. There exist two atomic limits, with Wannier centers at inequivalent high-symmetry positions in the unit cell. One of them, with Wannier center at the unit cell boundary (potentially realized in polyacetylene) leads to half charges at the end of an open chain 29,35,36 . How-ever, it is only with spectral symmetries imposed that systems with end modes are strong topological phases.
Obstructed atomic limits in 2D are the natural extensions of the above-mentioned topological phases in 1D. Twodimensional lattices, however, can exhibit richer classifications: even in the absence of bulk polarization, point-like corner charges can get generated 12,15,37 . From the exclusive perspective of charge fractionalization, these cases -which are protected only by spatial symmetries-broaden our understanding of second-order topological insulators [12][13][14][15][16][17][18]37 , which have the additional requirement of particle-hole or chiral symmetries. The relation between 2D Wannier centers and corner charge was developed in Refs. 30, 32, and 37. In particular, Ref. 32 uses the algebraic structure of the classifications of C n -symmetric insulators in class AI (spinless, time-reversal symmetric insulators) to build topological indices for corner charge. It has also been recently found that fragile topological phases can also host corner charges 25,32 and that 2D secondorder topological insulators also exhibit fractional charges at the core of defects with curvature singularities in both spinless and spinful insulators 32,38,39 .
In this work, we are concerned with 2D TRS spin-orbit coupled crystalline solids that admit a band structure description in terms of free fermions and fall in category (2) above, i.e., OALs. (This excludes in particular Z 2 topological insulators protected by TRS and phases with a mirror Chern number, where the mirror plane is the plane of the 2D solid itself.) Our aim is to classify OALs with fractional corner charges/filling anomalies in all 2D layer groups and to provide topological invariants which allow to infer the presence or absence of such charges from the knowledge of the bulk band structure alone. Such invariants are either computed from the irreducible representations of the Bloch states -in which case we speak of symmetry indicators -or expressed as integrals of a connection obtained from the Bloch states over (subsets of) the Brillouin zone, which will be referred to as Berry phase or Wilson loop type invariants. In particular, we make use of nested Wilson loop invariants 12 which were particularly helpful in diagnosing OALs for which symmetry indicators fall short. With the invariants presented here, we can identify OALs in a computationally more efficient way than by explicitly computing maximally localized Wannier functions. Importantly, the corner charge/filling anomaly of such a 2D system cannot be removed by any symmetry-respecting boundary manipulation, including the "gluing" of additional degrees of freedom to the boundary.
In general, a set of occupied bands (without strong or fragile topology) can be decomposed into subsets of bands stemming from localized orbitals at different Wyckoff positions. The minimal subblocks that cannot be further decomposed are elementary band representations (EBRs), which are a connected set of subbands induced from placing a certain orbital at a given Wyckoff position 20,22,[40][41][42][43][44][45] . Reference 20 introduced EBRs as a means to discern bands that stem entirely from atomic limits from those with strong or fragile topology. Following the approach in Ref. 32 correspondence between bulk invariants and corner charges sourced by OAL Wannier charge centers. In addition to the general classification, we discuss buckled bilayers of Bi, As, and Sb as a family of 2D materials that can realize OALs and corner charges. Based on density functional theory (DFT) calculations, we provide a phase diagram as a function of buckling strength (which may be controlled with a suitable substrate), identifying strong topological insulator (TI) phases and 2D OALs. One of the OALs supports corner charges, while the other one does not.
Our paper is structured as follows. In Sec. II we define the precise meaning of corner charge for our work as well as the role of the sample termination. In Sec. III, we show how OALs can be identified by bulk invariants. Finally, in Sec. IV the material candidates are discussed.

II. CHARGE QUANTIZATION IN TIME-REVERSAL SYMMETRIC SPIN-ORBIT COUPLED INSULATORS
The 2D phases we are interested in support no 1D gapless boundaries. The bulk gap may be populated by boundarylocalized midgap states, but those cannot be stabilized by TRS or crystalline symmetries. In this section, we will show that it is nevertheless often possible to diagnose a phase as an OAL in 2D via its corner charge fractionalization. We establish that this property remains invariant even when all midgap states are pushed out of the bulk gap and arbitrary symmetrypreserving boundary manipulations are allowed. 15,24,30,32,46 A. Quantization of the corner charge We consider 2D spinful insulating systems with TRS T (class AII in the Altland-Zirnbauer classification, T 2 = −1) and the spatial symmetries corresponding to a symmetry group G. We exclude first-order topological insulators so that the models we are studying are generically gapped even in a geometry with open boundary conditions. Additionally, we exclude insulators with bulk (TRS) polarization because those have edge-induced filling anomalies that scale with edge length and therefore result in metallic edges that preclude the existence of stable localized corner charges 32 .
We assume a tight-binding description of the system of interest. Denote by a 1 and a 2 the translation vectors corresponding to the decomposition of the G-symmetric lattice Λ into n-site unit cells S = {r 1 , . . . , r n }, where r i denotes the position of site i in the unit cell as measured from the unit cell origin r 1 ≡ 0. (Note that here and in the following, we only treat unit cells that are mapped to themselves under all available point-group symmetries, and do not cut through atomic sites. By these properties, a finite-size termination which does not cut through unit cells becomes possible.) We have We are considering tight-binding Hamiltonians of the form where µ, ν run over orbital degrees of freedom defined at each lattice site and c † vµ creates an electron in orbital µ at lattice site v. Hermiticity of H as well as the symmetry requirements posed by T and the symmetry group G imply relations among the Hamiltonian elements h vµ,wν which we implicitly assume to be fulfilled here and in the following.
Given a unit cell decomposition of Λ in terms of S, we define a trivial atomic limit by a Hamiltonian that is adiabatically deformable into one for which the implication v ∈ r∈S (xa 1 + ya 2 + r) w ⇒ h vµ,wν = 0 (3) holds for all choices of x and y, that is, there are no couplings between different unit cells.
To calculate corner charges we consider a finite system of |F| unit cells, via restricting H to a subset Λ F ∈ Λ (thereby obtaining H F ), which is given by We choose Λ F so as to retain all point group symmetries contained in G, a subgroup we denote by G F (it does not contain translations or nonsymmorphic symmetries). Then we consider a subset C ⊂ Λ F comprised of a minimal (but larger than 1) number of disjoint boundary regions that form an orbit under G F and contain an integer number of unit cells each. We choose C to cover all boundaries of Λ F . A particular boundary region c ⊂ C has charge where |n denotes an eigenstate of H F that is taken out of the occupied subspace occ bounded by E Fermi and we have |vµ = c † vµ |0 where |0 denotes the electronic vacuum. Since we only consider regions c that are related to each other by elements of G F , they have necessarily the same charge. Now, note that the charge of the full system is an even integer (given by |occ|), as is the charge of the complement of C, as long as we choose the regions in C large enough so as to ensure that the eigenstates localized in the complement are pure bulk-like in character and unaffected by the presence of a boundary. This is always possible when the linear extent by which C penetrates the bulk is much larger than the correlation length set by the bulk gap. We may then view the states contributing to the charge of the complement of C as states of a complete system of reduced size that has periodic boundary conditions and even integer charge. We thus deduce that Q c is quantized in even integer multiples of 1/q, where q = |C| denotes the number of elements in C. See Figs. 1 a and b for an example with threefold rotational symmetry.
We call Q c the corner charge since, in a pristine OAL, its fractional part derives from exponentially localized Wannier orbitals that are "cut through" by corners in the boundary of the system 30-32 : The Wannier orbitals in OALs are localized at maximal Wyckoff positions in the unit cell, and have shapes that respect the little group of their Wyckoff position. When a Wyckoff position lies on the boundary of the unit cell, the latter cuts through the respective Wannier orbital. The corner charge Q c can then be calculated conveniently and is equal to the volume that all occupied Wannier functions integrate to in c (where a single Wannier function is normalized to unit volume).
Note that Wannier orbitals which are cut through by edges instead of corners contribute to the TRS polarization 3,32,47 and thereby correspond to a charge that scales linearly with the extent of the boundary. The corner charge, on the other hand, stays constant as the thermodynamic limit is taken. It is thus well defined only in absence of TRS polarization.
Importantly, not all OALs have a fractional corner charge in all finite geometries. For example, as we will see for the 1a OAL discussed in section IV, sometimes there are no symmetry-preserving terminations that cut through Wannier functions (if only entire unit cells are retained), even though the latter are not centered at the atomic positions of the crystal.
Any trivial atomic limit has Q c ∈ 2Z for any such choice of boundary region: when different unit cells are not coupled to each other, the charge in each unit cell has to be equal to the total charge of the occupied subspace of H F={(0,0)} , which is necessarily an even integer. We may then define corner charge fractionalization as occurring in systems for which Q c mod 2 is equal to non-zero even integer multiples of 1/q (odd integer multiples are forbidden by TRS). Note that in this work we assume all systems with nontrivial corner charge to be given by OALs, which have a representation in terms of exponentially localized Wannier functions 20,48 . However, the corner charge formulas we supply in section III C apply equally well to fragile phases 21,22,[24][25][26][27][28] . These can always be adiabatically continued into OALs when other OALs are added. For a calculation of the corner charge in spinful materials, such a "trivialization" of a fragile phase becomes necessary only in the symmetry class that has C 4 rotational symmetry as its sole crystalline component, since this symmetry does not by itself allow for explicit corner charge formulas in terms of the elementary topological invariants we consider.
The classification of corner charge fractionalization in class AII and symmetry group G is given by the set of inequivalent Q c mod 2 that cannot be changed without breaking G F or closing the bulk gap. We present the classification for all layer groups G in section F of the Appendix.

B. Robustness of the corner charge
We now discuss to what extent symmetry-preserving edge manipulations can change the corner charge Q c defined in Eq. (5). We treat an edge manipulation as the introduction of an additional 1D system along the circumference of the finite 2D sample, and ask how the corner charges of the combined system, defined on the appropriately augmented Hilbert space, can differ from those of the original 2D model. Since charges are additive it is enough to determine the possible charges of the 1D system. In the following, we take Q to be the total charge of the 1D addition. It is even due to the requirement that we may only add complete and non-anomalous gapped 1D systems with TRS. We then use the remaining crystalline symmetries to derive further constraints on the charges Q c that the 1D system contributes to a boundary region c.
We note that the point group symmetries in 2D that G F can contain are mirror and n-fold rotational symmetries, where n ∈ {2, 3, 4, 6}. We first discuss the latter case of C n rotational symmetries. For spinful systems with TRS, we have (C n ) n = −1. Let H 1D denote a general 1D TRS gapped Hamiltonian defined on a Hilbert space of L lattice sites (with L/n an integer), possibly augmented by orbital degrees of freedom (see also Fig. 1 c). A C n rotational symmetry implies that we can choose the order of regions c i ∈ C (which in combination cover all of the L sites of the 1D system) such that in real space the symmetry effects c i → c i+1 mod n , that is, a translation by L/n sites. Now, due to (C n ) n = −1, rotations are equivalent to translations around a 1D circle that encloses a π-flux. Let t be the operator for translations by a single site, i.e., it shifts site r ∈ {1, . . . , L} of the 1D lattice to site r + 1 mod L. It is not a symmetry of H 1D , however, we can obtain a t-symmetric Hamiltonian (on a ring enclosing a π-flux) by adding up L/n copies of H 1D that are subsequently shifted by one lattice site, to arrive at which acts on an L/n-fold enlarged Hilbert space. The occupied subspace of H TRN 1D has a total charge of QL/n and enjoys a translational symmetry that corresponds to L repeated unit cells, with twisted boundary conditions so as to accommodate the π-flux. It is gapped and has TRS just as H 1D , and its charge thus necessarily corresponds to an even integer number of filled Bloch bands, which each hold L states. We conclude that its charge per unit cell Q/n is an even integer. Returning our attention to H 1D , since all boundaries c carry the same charge, this is exactly the corner charge Q c = Q/n. Thus in the case of C n -symmetries there is no 1D addition that can trivialize the fractional corner charges of a 2D OAL.
Next, we turn to mirror symmetries, which for spinful systems satisfy M 2 = −1. In the case of two reflections, say M x and M y , we also have a two-fold rotation symmetry C 2 = M x M y , which by the argument above allows us to conclude that all corner charges Q c contributed by any gapped and TRS 1D addition are necessarily even (note that the minimal nontrivial boundary decomposition has q = 2). When there is only a single mirror symmetry, we cannot argue along these lines, since it does not act on the 1D real space as a translation. In fact, it "translates" different sites along the 1D chain by different amounts. Hence, here the symmetry constraint on Q c is the same as that for the 2D bulk, namely that Q ci ∈ Z, i = 1, 2, (compare this to the Q ci ∈ 2Z we obtain for C 2 symmetry) and a fractional charge of 1 mod 2 can be trivialized.
Finally, we note that in the case where we have C 3 symmetry as well as 3D inversion symmetry I (which is the same as C 2 M z symmetry), we can define an effective 1/6 translation by t 1/6 = IC 2 3 which allows to argue that patches c of size 1/6 of the linear extent of the full 1D system have even integer charge. This is important for the robustness of the Q c = 1/6 corner charges of this symmetry class.
Since any finite-size geometry breaks the remaining nonsymmorphic symmetries a system might have, we do not need to consider their effect on charge fractionalization. We conclude that quantized corner charges can be changed by 1D edge manipulations only in the case of a single mirror symmetry.

III. IDENTIFICATION OF OBSTRUCTED ATOMIC LIMITS
In this section we give a prescription for obtaining the corner charge of the occupied subspace of a bulk model represented by a Bloch Hamiltonian H(k) [the Fourier transform of the translationally invariant tight-binding model given in Eq. (2)], assuming that its occupied subspace realizes an OAL. We take a Wannier center point of view: in particular, we define an OAL by the way it is built up from exponentially localized and symmetric Wannier functions 20,48 .
As shown in section II B, mirror symmetries can protect fractional corner charges only when they are combined to yield a twofold rotational symmetry. The protecting symmetries we consider are therefore C n rotations, with or without an additional 3D inversion symmetry I. The inclusion of I symmetry allows us to extend our discussion to the experimentally relevant case of 2D honeycomb monolayers with nonzero buckling, and our classification (given in the Appendix) to the 80 layer groups instead of the 17 wallpaper groups. We note that inversion effectively replaces C 2 in its role of enforcing a Q c = 0, 1 mod 2 quantization of the corner charge, but due to I 2 = +1 [whereas (C 2 ) 2 = −1] allows for symmetry indicator invariants. Furthermore, in the case of C 4 symmetry, we find that we require an additional inversion symmetry in order to be able to read off the corner charge from the available topological invariants. Inversion symme- try is however, unlike C 4 , not necessary for the topological robustness of the corner charge. We first establish which topological invariants can be defined for each point group in III A. We then go on to calculate these indices for the elementary band representations of each symmetry class in III B. Finally, in III C we present formulas that allow for a determination of the corner charge in all symmetry classes except for the one that as its only crystalline component has C 4 rotational symmetry.
A list of all possible corner charges for the various layer groups is given in Table XI of the Appendix. The remaining symmetry operations a layer group may contain in addition to the ones listed in section III C are either irrelevant for finitesize corner terminations since they involve translations (as in the case of non-symmorphic symmetries), or merely impose constraints for the shape of the finite-size termination, without affecting the corner charge quantization itself (such as mirror symmetries).

A. Bulk topological indices
To identify different EBRs, we employ a combination of symmetry indicator 20,32,40,41,[49][50][51][52][53][54] and Wilson loop 22,40,41,[55][56][57][58][59] topological invariants. These can be evaluated from the crystal's Bloch Hamiltonian H(k) and so do not require a realspace calculation to be performed. The main ingredient for both kinds of invariants is the bundle of occupied Bloch states Sets of allowed eigenvalues for spinful rotational symmetries. a C2 symmetry. b C4 symmetry. c C3 symmetry. The possible eigenvalues of C6 symmetry are not shown, since they do not allow for the definition of symmetry indicators (there is at most one C6symmetric point in any two-dimensional Brillouin zone).
[k is an element of the first Brillouin zone (BZ) of the crystal.] Given a unitary crystal symmetry S that is realized on the Bloch Hamiltonian as SH(k)S † = H(Sk) and acts on the momenta as k → Sk, we can calculate its corresponding symmetry indicator topological invariants from the eigenvalues of the matrices where m, n run over the occupied subspace only andk = Sk are high-symmetry points (HSPs) of the Brillouin zone that are left invariant by the symmetry S (see also Fig. 2). An nfold symmetry acting on spinful fermions satisfies S n = ±1 (positive sign for 3D inversion, negative sign for 2D mirror and rotational symmetries), this, together with TRS, imposes constraints on the possible eigenvalues of S mn and allows for the definition of topological invariants that capture the different symmetry representations of the occupied bands across the BZ. A trivial OAL, being deformable to a momentum-independent Hamiltonian, will have the same representation across HSPs that are invariant under the same symmetry and hence will have trivial symmetry indicator invariants. Nonzero symmetry indicator invariants, on the other hand, indicate that the bands adopt different representations of the symmetry across the BZ and correspond to nontrivial OALs. The Wilson loop (along a closed, non-contractible path γ in the BZ that starts and ends at the momentum k * ) is an operator on the filled band subspace of H(k) defined as where P (k) = m∈occ |u m (k) u m (k)| is the projector onto the subspace of filled bands at momentum k. Note that we choose a gauge where H(k) = H(k + G) for a reciprocal lattice vector G and the product is path-ordered along γ. The Wilson loop operator satisfies W γ W † γ = P (k * ) and so, since any projector satisfies [P (k)] 2 = P (k), its eigenvalues are either zero or of the form e iθ γ α , α = 1 . . . N . In the following, we refer to the set of {θ γ α } α=1...N as the Wilson loop spectrum, suppressing the zero eigenvalues.
The anti-unitary TRS T acts on the Bloch Hamiltonian as T H(k)T −1 = H(−k). For the projectors this implies T P (k)T −1 = P (−k). When γ is mapped onto itself by TRS, and its starting point satisfies k * = −k * up to a reciprocal lattice vector, we then have Due to T being anti-unitary and T 2 = −1, this implies a Kramers degeneracy of the Wilson loop spectrum, i.e., every θ γ α is (at least) two-fold degenerate when γ is mapped onto itself by time reversal. Now, if there is a crystal symmetry S that reverses the direction of γ and leaves the starting point invariant so that k * = Sk * up to a reciprocal lattice vector, we have Since S is unitary, the Wilson loop is unitarily equivalent to its complex conjugate and so its eigenvalues come in complex conjugated pairs. This implies a symmetry of the Wilson loop spectrum around θ = 0, for every θ γ α there is a corresponding −θ γ α mod 2π.
We may furthermore employ nested Wilson loops 12,15 . Let labels a point in the two-dimensional BZ in some basis (chosen such that k i,j = 0 and k i,j = 2π are related by reciprocal lattice vectors). Consider the Wilson loop Hamiltonian Equations (10) and (11) then imply where we defined We see that T implies a TRS of the Wilson loop Hamiltonian, whereas S implies a particle-hole symmetry. These properties are needed for the definition of quantized topological invariants of the nested Wilson loop: We define W b i as the Wilson loop calculated from a gapped set eigenstates b of H Wi (k j ) along a closed, non-contractible path k j : 0 → 2π in the reduced BZ.
We differentiate between three kinds of nested Wilson loops that differ by the choice of the set of eigenstates b: 1) The nested loop W 0 i , which is calculated from the two bands in the spectrum of H Wi (k j ) that at k j = 0, π have a degeneracy pinned to the Wilson eigenvalue 0 [note that any such degeneracy at k j = 0 implies one at k j = π and vice versa due to the absence of Wannier center flow in (obstructed) atomic limits]. 2) The nested Wilson loop W λ i , calculated for the upper or lower half of the bands in the spectrum of H Wi (k j ) that are not pinned at k j = 0, π to a Wilson eigenvalue 0, π (that is, half of the freely dangling Wilson bands, which by the particle hole symmetry come in pairs). 3) The nested loop W π i , which is calculated from the two bands in the spectrum of H Wi (k j ) that at k j = 0, π have a degeneracy pinned to the Wilson eigenvalue π.
The nested Wilson loops of type (1) and (3) cannot be trivialized by transformations that preserve S and T and are adiabatic with respect to the bulk gap. The reason is that the invariants calculated from these loops are equal to the partial polarizations of Wilson bands pinned to eigenvalues 0 or π by S at the transverse momenta k j = 0, π. Wilson gap closings that preserve the energy gap can only occur in pairs (due to S) at intermediate transverse momenta k j , −k j . It is rigorously shown in Appendix A of Ref. 60 that these gap closings together always contribute integer multiples of 2π to the nested partial polarization, and therefore cannot trivialize Wilson loops of type (1) and (3). In this work, we do not consider invariants derived from Wilson loops of type (2).
We will now list the topological invariants that can be defined for a given point group. We find that often the inclusion of I symmetry allows for the replacement of Wilson-loop invariants by symmetry indicator invariants. For the discussion of symmetry indicators, we make use of the definitions and derivations presented in sections A-C of the Appendix.
Note that in the following, and as motivated at the beginning of section II A, we explicitly exclude invariants that characterize topological insulators because they are necessarily gapless along the edges in a 2D geometry with open boundary conditions, and so do not allow for stable quantized corner charges. In addition to removing some invariants from our analysis altogether, this imposes constraints on Wilson loops.
We emphasize that our list of invariants may not be exhaustive. As noted in Ref. 22, it is in general difficult to identify all possibly nontrivial Wilson loop invariants. In this work we only treat "straight" (nested) Wilson loops, which (given a starting point) go around one of the two inequivalent noncontractible loops of the Brillouin zone torus.
1. C2 symmetry a. Symmetry indicator invariants The BZ has four highsymmetry points (HSPs, defined in section B of the Appendix), see also Fig. 2 a. All the points are invariant under C 2 . Thus, they all have C 2 eigenvalues +i, −i (see Fig. 3 a). However, since all the HSPs are also time-reversal invariant momenta (TRIMs), the eigenvalues have to come in complexconjugate pairs, leading to a single available 2D irreducible representation. Therefore, the C 2 eigenvalues on their own do not afford a topological distinction and there are no symmetry indicator invariants.
b. Wilson-loop invariants For every closed highsymmetry line γ (which connects two HSPs) of the 2D BZ that is left invariant by C 2 , we can define a Wilson loop that is TRS and C 2 symmetric. Due to Eqs. (10) and (11) the parities of the numbers of θ γ α = 0 and θ γ α = π eigenvalues in its spectrum cannot be changed under adiabatic deformations of H(k): adiabatic perturbations of the Hamiltonian at most move particle-hole related Kramers pairs of eigenvalues in and out of 0 and π, this does not change the total parity.
A topological invariant of W γ with spectrum {θ γ α } α=1...N is therefore given by where the product is taken over only one eigenvalue of each Wilson loop Kramers pair. We call ν γ = 0 trivial and ν γ = 1 nontrivial. This invariant is equivalent to the TRS polarization 47 and counts the parity number of Wilson loop pairs of eigenvalues equal to π. We also define which counts the number of Wilson loop pairs of eigenvalues equal to 0. The invariants ν γ and µ γ are not independent when the total number of bands is fixed. They obey Therefore, we drop µ γ as it provides redundant topological information. In the following, we will consider Wilson loops that go through high-symmetry points in the 2D BZ. We denote by ν AB the loop that goes from point A to point B and then back to A via the shortest non-contractible loop around the BZ torus.
There are in total four TRIMs and three topologically inequivalent straight and C 2 -symmetric Wilson loops. This can be seen by noting that, holding one of the four C 2 -symmetric momenta fixed as a starting point, there are two incontractible loops around the Brillouin zone torus (which necessarily go through one other C 2 -symmetric momentum). Keeping in mind that path-reversed Wilson loops are not independent [as per Eq. (11)], this naively yields the set of Wilson loop invariants {ν ΓX , ν ΓY , ν XM , ν Y M }. We note however that the path denoted by Y Γ−ΓX −XM is topologically equivalent to the path denoted by Y M and so we have The remaining invariants are further constrained due to the requirement that the Z 2 TI invariant ∆ TI vanishes: we have that We are left with two Wilson loop invariants. Similarly, we may define the quantized invariants ν 0,π x,y and µ 0,π x,y from the nested Wilson loops W 0,π x,y , since these are calculated for particle-hole symmetric sets of bands 18,25 [in contrast to W λ x,y , which does not satisfy Eq. (13)]: the anticommutativity with the Wilson loop Hamiltonian that distinguishes particle-hole symmetry from a reflection symmetry is irrelevant from the point of view of the nested Wilson loop, as long as the latter is defined via a projector onto a particle-hole symmetric set of bands. We may therefore define ν 0,π x,y and µ 0,π x,y just as in Eqs. (15) and (16), where θ γ α this time refers to the spectrum of the nested Wilson loop. As before, we drop the µ invariants since they are not independent when the number of occupied Wilson bands is held fixed. Taking into account the constraints 60 reduces the number of independent invariants to three. The third equation can be seen in the following way: ν π x is nonzero if and only if the occupied subspace hosts an odd number of Wannier Kramers pairs whose centers are shifted by 1/2 in both x and y direction (taking the lattice constant to be 1) with respect to the center of the unit cell, i.e., if there is an odd number of Kramers pairs at Wyckoff position 1b of the crystal (see also Table II). This Wyckoff position stays unchanged when exchanging x and y, we therefore obtain that ν π x is nonzero if and only if ν π y is nonzero. Note that the corresponding statement does not hold for ν 0 x and ν 0 y , since these indicate Wannier Kramers pairs at the 1c and 1d Wyckoff positions, respectively.
We therefore choose the classification to be given by c. With inversion symmetry C 2 + I symmetry is equivalent to I symmetry for all our purposes. Inversion symmetry allows us to replace the Z 2 valued Wilson-loop invariants by 2Z-valued symmetry indicators. The BZ has the I-invariant points Γ, X, Y and M , which support the six inversion eigenvalue invariants where #X I i (#Γ I i ) is the number of occupied states with inversion eigenvalue X I i (Γ I i ), and X I i=1,2 , Γ I i=1,2 = {1, −1}, and similarly for Y and M . Due to the fixed number of occupied bands, we have the constraints The three remaining invariants completely fix 59 the Wilson loops in Eq. (21). Due to TRS they are necessarily even integers. We retain the classification 2. C3 symmetry a. Symmetry indicator invariants The BZ only has the C 3 -invariant points K and K , see also Fig. 2 c. Now, we discuss the invariants that compare the representations at the K (K ) and Γ points of the BZ, where K i=1,2,3 = {e iπ/3 , −1, e −iπ/3 }, and similarly for K (see Fig. 3 c). Unlike M , the HSP K is not a TRIM. Instead, TRS relates K and K . TRS imposes the constraints, The six invariants are subject to the constraints (26) along with due to the constant number of occupied states across the BZ. The symmetry-indicated part of the classification is given by the two invariants 2 ]}. where In addition we have We retain [M I 2 ] as the invariant that determines the classification, in addition to the C 3 invariant [K 3. C4 symmetry a. Symmetry indicator invariants The BZ has four HSPs (Fig. 2 b). Two of them are invariant under C 2 and give rise to trivial indicators due to time reversal symmetry. We can then only build indices that compare the C 4 symmetry representations at M with those at Γ as follows: where the eigenvalues are taken from the set M With the constraints in (34) and (35), we eliminate the redundant invariants [M We conclude that the classification is given by c. With inversion symmetry The invariants given in Eq. (24) (together with the C 4 constraint [X I 2 ] = [Y I 2 ]) allow us to replace ν ΓX , ν Y M . We conclude that the classification with inversion symmetry is given by 4. C6 symmetry a. Symmetry indicator invariants In a C 6 -symmetric BZ, there are two inequivalent HSPs, M , which is invariant under C 2 , and K, which is invariant under C 3 (Fig. 2 d). All other points are related by rotations, and thus provide redundant representations for the purpose of classification. Furthermore, M is both a HSP and a TRIM. Thus, from the analysis of the previous classifications, no invariants can be derived from its representations. Now, we discuss the invariants that compare the representations at the K and Γ points of the BZ, where K i=1,2,3 = {e iπ/3 , −1, e −iπ/3 }. Unlike M , the HSP K is not a TRIM. Instead, TRS relates K and K . TRS imposes the constraints, But the representations at K and K are the same due to C 6 symmetry, The last two sets of constraints leave us with only two nonredundant invariants, [K 1 ] and [K 2 ]. However, due to the constant number of occupied states, we have i #K 2 ] = 0, which makes one of these invariants redundant too. We choose the symmetry-indicated part of the classification to be given by [K b. Wilson-loop invariants C 6 symmetry implies having C 2 symmetry as well and so we can define µ ΓM as an invariant due to Eq. (16). We choose µ ΓM here instead of ν ΓM since it directly indicates Wannier centers at the 3c Wyckoff position (see Fig. 4 d and Table V) of the hexagonal unit cell. We do not consider nested Wilson loops in this symmetry class because the corner charge can be completely determined without them. In conclusion, we have c. With inversion symmetry The BZ has the I-invariant points M , M and M , which support the invariants where In addition we have We retain [M I 2 ] as the invariant that determines the classification. Due to C 6 symmetry 32 and TRS, [M I 2 ] ∈ 4Z. We conclude that χ (6)

B. Decomposition into EBRs
Tables II-V list the EBRs 20,22,40-45 supported by systems with C n rotational symmetry, together with their invariants and corner charges. The minimal block sizes correspond to the multiplicities of the respective Wyckoff positions (multiplied by two to account for spin). If multiple choices for the site-symmetry group 20 representation at a Wyckoff position W are available, we denote the representation with eigenvalues e iα as W | α .
We show in Sec. E of the Appendix how the symmetry indicator invariants for different EBRs can be derived. The (nested) Wilson loop invariants can be obtained by the mapping of Wilson loop spectra to Wannier centers 48 .   Fig. 4 a), and their invariants. All atomic limits can be decomposed into EBRs formed by single Kramers pairs.
EBRs with C4 symmetry induced from the maximal Wyckoff positions listed in the first column (see Fig. 4 b), and their invariants. All atomic limits can be decomposed into EBRs formed by at most two Kramers pairs. Importantly, the (non-elementary) band representation 1b| ± π 4 ⊕ 1b| ± 3π 4 has trivial C4 invariants but nonzero corner charge. In systems with C4 as the sole crystalline symmetry, this obstructs a determination of the corner charge in terms of topological invariants.

C3
[K  Fig. 4 c), and their invariants. All atomic limits can be decomposed into EBRs formed by single Kramers pairs.   Fig. 4 d), and their invariants. All atomic limits can be decomposed into EBRs formed by at most three Kramers pairs.

C. Formulas for corner charges
In this section, we provide explicit formulas for the corner charge in terms of the topological invariants as evaluated on the entire occupied subspace of a given model. For systems with I, C 3 , and C 3 + I symmetry, we can uniquely identify the spinless limit of a given spinful model. In this case we can employ the results of Ref. 32. In the remaining cases we deduce the formulas from the EBR tables given in section III B. Importantly, all corner charges appearing in these formulas as well as in the EBR tables apply only to crystal terminations where Λ F in Eq. (4) has corners at the intersection of 1D edges that are obtained from translating unit cells with crystal lattice vectors 32 , but not necessarily primitive ones.
As noted before, in the case where we only have C 4 symmetry at our disposal, no corner charge formula can be constructed from our invariants. We leave the investigation of this symmetry class to future work, and instead consider the case of C 4 + I symmetry here.

I symmetry
Inversion symmetry becomes equal to C 2 symmetry in the spinless case. This means that, using inversion eigenvalues, we can uniquely read off the C 2 eigenvalues of the spinless version of any model at hand, and may then use the formula presented in Ref. 32 for spinless C 2 symmetry to infer the corner charge of our model. Note that the doubling of the corner charge, which comes with going from spinless to spinful and imposing TRS, is automatically taken into account by the fact that the inversion eigenvalues are equal for Kramers partners. We therefore obtain A nonzero value implies two equal fractional corner charges at I-related sectors with Q c = 1.

C2 symmetry
Comparing with Table II, we have where, if H Wx (k y = 0, π) does not have pinned bands at eigenvalue π, we declare ν π x = 0. We note that ν π x is Z 2 valued, in accordance with the fact that two Wannier Kramers pairs at 1b are trivial in that they can be removed from 1b and moved around the unit cell in a C 2 symmetric fashion. A nonzero value of Q c = 1 implies two equal fractional corner charges at C 2 -related sectors.

C3 symmetry
As shown in Appendix D, there is a one-to-one mapping between the C 3 eigenvalues of the spinless and spinful cases. It implies that A nonzero value implies three equal fractional corner charges at C 3 -related corners, with possibilities Q c = 2 3 or Q c = 4 3 .

C3 + I symmetry
The one-to-one mapping of C 3 eigenvalues from Appendix D, as well the observation that inversion symmetry becomes the same as C 2 symmetry in the spinless case, yields A nonzero value implies six equal fractional corner charges at C 3 , I-related corners, with possibilities Q c = 1 3 , Q c = 2 3 , Q c = 1, Q c = 4 3 , or Q c = 5 3 .

C4 + I symmetry
While I symmetry only allows for the decomposition of the sample into two halves, and therefore for a corner charge quantized in units of 1 mod 2, C 4 symmetry affords a further halving, so that the corner charge is quantized in units of 1/2 mod 2. Any I protected corner charge can in this way be split up into two C 4 + I protected corner charges of half the size. Using Eq. (50), we therefore obtain which we simplified by the . A nonzero value of Q c implies four equal fractional corner charges at C 4 -related corners (this configuration is automatically I symmetric), with possibilities Q c = 1 2 , Q c = 1, Q c = 3 2 .

C6 symmetry
Comparing with Table V, we have where µ ΓM denotes the parity of the number of W ΓM zero eigenvalue pairs. A nonzero value implies six equal fractional corner charges at C 6 -related corners, with possibilities Q c = 1 3 , Q c = 2 3 , Q c = 1, Q c = 4 3 , or Q c = 5 3 .

Summary
S Qc TABLE VI. Summary of corner charge formulas.

IV. MATERIAL CANDIDATES
We propose the group-V buckled honeycomb monolayers of elemental antimony (Sb) and arsenic (As) as material realizations of protected fractional corner charges. Theoretical studies suggest that antimonene and arsenene can serve as an excellent platform for electronics due to high band gap tunability and mechanical stability [61][62][63][64][65] . Moreover, these 2D materials, as well as atomically thin bismuth monolayers (called bismuthene), deposited on a SiC substrate, are promising candidates for a realization of the quantum spin Hall states at room temperature [66][67][68] . Only recently, several experimental reports have demonstrated a successful fabrication of a monolayer structure of antimony [69][70][71] and arsenic 72 .
Free-standing monolayers with nonzero buckling d z have a three-fold rotational symmetry C 3 as well as inversion I symmetry (consult Fig. 5 d). (In practice, we consider weak substrate coupling so that the inversion symmetry is approximately retained.) Applying strain leads to a decreasing d z parameter up to a fully flat structure with six-fold symmetry. In Fig. 5 a, b, c, we present the band gap evolution of Bi, Sb and As as a function of tensile strain, which is modeled by a modification to the in-plane lattice parameter (larger strain corresponds to a longer in-plane distance between atoms). First, we note the qualitative similarity of the phase diagrams for all three investigated materials. At d z = 0 (which corresponds to a large strain around ∼25%), there is an additional mirror symmetry M z , and all structures are in an topological crystalline insulating (TCI) phase, protected by a mirror Chern number, which we verified by Wilson loop calculations (not shown here). This phase does not have exponentially localized Wannier functions that respect all symmetries of the model. Small buckling breaks the mirror symmetry and the materials then realize an OAL with localized Wannier orbitals centered at the center of the hexagons in the honeycomb lattice (Wyckoff position 1a of the crystal). Upon further decreasing strain, a transition to a Z 2 topological insulator (TI) is observed via a band gap closing around d z = 0.6Å. To confirm this topological phase transition, we compute the Z 2 topological index ∆ TI given by the product of the inversion eigenvalues of the occupied bands at the time-reversal invariant momenta in the BZ 4 , and obtain ∆ TI = 1. As strain decreases further, another band gap closing occurs. Here, the Bi monolayer reenters a TI phase (with different symmetry indicator invariants as shown in Table VII), as confirmed by the Z 2 index remaining nontrivial. In contrast, the almost fully buckled Sb and As monolayers enter once again in an OAL phase, this time with bands induced from the Wyckoff positions 3c (which is located on the bonds of the hexagon, see Fig. 4 d). Hence, our results reveal more details on the previously investigated strain-induced topological phase transitions in these materials [73][74][75] .
Let us consider the systems with open boundary conditions. To establish the presence of corner charges, we perform open flake calculations for distinct OALs using the localized basis DFT method SIESTA 76 . In Fig. 5 e, f, we show results for a fully buckled antimony flake as a representative of the 3c OAL. The most direct indicator of fractional corner charges are corner-localized midgap states. If present, they are expected to appear close to the Fermi level. However, they are not necessarily well-separated from the bulk or edge modes. Therefore, we passivate the structure with tellurium atoms (marked with stars in Fig. 5 f) in order to remove spurious dangling edge states from the bulk gap. The energy spectrum (see Fig. 5 e) then exhibits 12 exactly degenerate corner states at the Fermi level, with only half of them filled. We thus obtain a fractional corner charge of Q c = 1 mod 2 per corner, realizing a filling anomaly, as at the given filling it is not an insulating state that satisfies both charge neutrality and the crystalline symmetries. 32 We confirm this corner charge using the topological indices developed in section III A. In Table VII, we evaluate the symmetry indicators for all discussed phases. We may then compute the corner charge of the 3c OAL on a hexagonal flake using Eq. (49). The relevant unit cell is the hexagonal cell, shown in Fig. 4 I = (4, 0) for the hexagonal cell, from which we obtain Q c = 1 mod 2 by Eq. (49). This is in agreement with the numerical results presented in Fig. 5 e, f. Correspondingly, in the case of the 1a OAL, we obtain χ (3) I = (0, 0) for the hexagonal cell (the primitive cell cannot be used to build a C 3 -symmetric finite geometry). We conclude that there are no fractional charges. This is a case in point: although the 1a atomic limit is obstructed, in the sense that the electrons are localized away from the atomic sites, which are located at the 2b Wyckoff position of the crystal, there are no protected corner charges. (There may however be such charges in C 3 -symmetric geometries that are terminated by cutting through unit cells. We do not consider these geometries here, mainly because there is no bulk-boundary correspondence in this case, and the actual corner charge is dependent on how the boundary unit cells are cut.) phase I corresponding to different regions in the phase diagrams. The symmetry indicators were calculated using the primitive 2-site unit cell of the honeycomb lattice (see Table XII for a decomposition in terms of elementary band representations). The indices χ I allow for a more refined classification even of the strong TIs. We find that the 3c and 1a OALs differ in their inversion indicator [M I 2 ] and thus, as explained in the main text, inversion-symmetric flakes built from their hexagonal unit cells differ by a protected corner charge equal to 1 mod 2.

V. DISCUSSION
As established in Refs. 29, 35, and 36 for insulators with bulk polarization and more recently in Refs. 30 and 32 for second-order topological insulators, the nontrivial bulk topology of OALs can be revealed via charge fractionalization at boundaries. This represents the simplest mechanism for a topological bulk-boundary correspondence that is protected by crystalline symmetries. In this work, we presented theory and material candidates for charge fractionalization at corners in 2D systems with significant spin-orbit coupling, thus providing a broader picture than the one presented in some recent previous treatments of this phenomenon [30][31][32] .
Corner charges in topological insulators are well-defined when there is no edge spectral flow but also only in the absence of an edge-induced filling anomaly due to (timereversal) bulk polarization. Since there is no crystalline symmetry-protected edge spectral flow in 2D (assuming the symmetry acts at least in part non-locally), corner charges are well defined for all 2D systems that are not strong or weak first-order topological insulators, or M z mirror Chern insulators 7,10 .
Diagnosing spinful OALs with time-reversal symmetry in 2D was particularly challenging because the irreducible representations of the occupied bands at HSPs are usually twodimensional, yielding trivial symmetry indicator invariants at C 2 -invariant HSPs. Symmetry indicators were therefore insufficient to identify the Wannier centers in C 6 , C 4 , and C 2 symmetric insulators. This is straightforwardly manifested in the fact that inducing a band representation from spinful maximal Wyckoff positions exhausts the symmetry representations at various HSPs of the BZ. To overcome this difficulty, we considered Wilson loop and nested Wilson loop invariants, which could better "resolve" the positions of the Wannier centers. Wilson loops, however, are essentially one dimensional objects that extract projections of the 2D positions of the Wannier centers along particular directions. Nested Wilson loops are a best-effort attempt to localize the Wannier centers in 2D, but cannot always be interpreted literally due to the possible non-commutation of Wilson loops along different directions. In the presence of crystalline symmetries, however, Wilson and nested Wilson loops have eigenvalues with quantized phases, which clearly distinguish different OALs in C 6 and C 2 symmetric insulators, but are insufficient for insulators which only have C 4 symmetry. We leave the challenge of finding a formula for the corner charge in such systems to future work.
We studied the protection due to only spatial symmetries because corner charge fractionalization is a robust observable that does not require additional spectral symmetries such as chiral or particle-hole symmetry. However, when particlehole symmetry is present, we can additionally predict topologically protected zero-energy corner states. These are characterized by Q c = 1 mod 2: Consider a system with an n-fold symmetry in a phase with 2n degenerate midgap states (the 2 is due to TRS). At half-filling, n midgap states are occupied and there is no gap. To arrive at a gapped system (as required for the corner charge to be well defined), we need to either fill n more states or remove n electrons from the charge-neutral system. When maintaining the crystal symmetry, this implies an excess (or missing) charge of Q c = 1 mod 2 for each of the n corners 32 .
Interestingly, we find that there are obstructed atomic limits, where the electrons are localized away from the atomic sites, which still do not have nontrivial corner charges. These may instead be diagnosed by their response to crystal defects 32,39 . We leave the exploration of the defect response of obstructed atomic limits with significant spin-orbit coupling to future work. r n = −1. From (A1), it follows that h(Rk)r |u n k =rh(k) |u n k = n (k)r |u n k . (A2) Thus,r |u n k is an eigenstate of h(Rk) with energy n (k). This means that we can write the expansion where is the sewing matrix, which is unitary: As before, let us use (A1) to do the following calculation: Thus, it is possible to choose a gauge in which the energy eigenstates are also eigenstates of the rotation operator, This is automatic if there are no degeneracies, but if energy degeneracies exist, one can always choose a gauge such that the above expression is possible. At these invariant points, the sewing matrix is diagonal: Now, we show that the rotation eigenvalues of HSPs that are related by symmetry are equal. Consider the rotation by an angle φ in a crystal with C 2π/φ symmetry. This rotation symmetry relates HSPs that are invariant under rotations by a larger angle θ = nφ, for n integer. Call these HSPs Π θ . Here, we are interested in knowing how the rotation eigenvalues of Π θ and R φ Π θ are related. In particular, this applies to two cases: (i) In C 6 -symmetric crystals, φ = 2π/6. For θ 1 = 2π/3 = 2φ we have K = R φ K , while for θ 2 = π = 3φ we have M = R φ M = R 2 φ M ; (ii) in C 4 -symmetric crystals, φ = π/2, for θ = π = 2φ we have X = R φ X. Let us start by asking what we get from applyingr θ |u n R φ Π θ . Since R φ Π θ is invariant underr θ , we havê Since R φ Π θ and Π θ are related by C 2π/φ symmetry, we can expand where B mn Π θ = u m R φ Π θ |r φ |u n Π θ is the sewing matrix, with the properties shown before. Conversely, we also have that So, replacing this expansion in (B5), we havê Taking a different approach, we calculate directly the rotation eigenvalues in the expansion (B7) to get So, comparing the last two results we conclude that for all n. Furthermore, since the eigenstates form an orthonormal basis, we have for all m and n. Now, the sewing matrix will have non-zero elements for equal energies at the two different HSPs R φ Π θ and Π θ . Thus, for m (R φ Π θ ) = n (Π θ ), we need r m Π θ = r n R φ Π θ , i.e., the rotation spectra at R φ Π θ and Π θ are equal for bands having equal energies. In particular, we have the relations Thus, on one hand we have where V is the sewing matrix for TRS. On the other hand, we havê In the last expression, we have used the fact that R(−k) = −Rk. From these two expressions we conclude that for all l. Since the eigenstates are orthonormal, this relation implies that for all m, l. As noted earlier, of particular interest are the HSPs. At these points, B mn Π = r n Π δ mn in the gauge in which {|u n Π } are rotation eigenstates. Then, at these points, the previous relation results in , we have that V ml Π = 0, which means that there is no restriction on the rotation eigenvalues]. In particular, at TRIMs which are also HSPs, Π = −Π, we have that r l * Π = r m Π for equal energies m (Π) = l (Π). This imposes the following constraints on the rotation eigenvalues: (i) for a non-degenerate state labeled by n, r n * Π = r n Π , i.e., its rotation eigenvalue is real: r n Π = ±1 and (ii) for two degenerate states n = 1, 2 one could have r 1 Π = λ and r 2 Π = λ * , so that r 1 * Π = λ * = r 2 Π and r 2 * Π = λ = r 1 Π , that is, in energy-degenerate states, the rotation eigenvalues can be complex, but have to come in complex conjugate pairs. As said before, these constraints follow for HSPs that are also TRIM. This is the case for all the HSPs except K and K , which map into each other under TRS.

Appendix D: Mapping between spinless and spinful C3 eigenvalues
We start with the spinless indicators whereK i=1,2,3 = {1, e i2π/3 , e −i2π/3 }. Upon introducing spin, each spinless eigenvalue λ contributes two spinful eigenvalues λe ±iπ/3 . From this we obtain the relations In this Section, we explicitly induce the energy band representations at HSPs of the BZ following the prescription in Ref. 43. Given a site symmetry representation, the induced band representation will allow us to identify the symmetry indicator invariants associated with a maximal Wyckoff position. In this section, we induce the band representations and corresponding symmetry indicator invariants for all the allowed site symmetry representations of spinful time-reversal symmetric orbitals at each maximal Wyckoff position. In the following, ρ refers to the representation of the site symmetry group, while ρ k G refers to the band representation at crystal momentum k. We treat each case separately.

C4 symmetry: Representations induced from 2c
Given a site symmetry representation ρ(C 2 ) of the orbitals at 2c, the band representations are Let us consider one the only possible site symmetry representation, ρ(C 2 ) = e i π 2 σz . For C 4 , the band representations at HSPs are Both of these matrices have the four eigenvalues e iπ/4 , e −iπ/4 , e 3iπ/4 , e −3iπ/4 . Therefore, [M 1 ] = 0. For C 2 , the band representations at HSPs are All these matrices have eigenvalues +i,+i, −i, −i, also leading to vanishing symmetry indicators. As we will see, this is also the case when the band representations are induced from 1b: in fact, with spinful time-reversal symmetry, the only possible EBR is given by a pair of states with C 2 eigenvalues (+i, −i). Therefore, no symmetry indicators exist associated with the band representations of C 2 .

C4 symmetry: Representations induced from 1b
Given a site symmetry representation ρ(C 4 ) of the orbitals at 1b, the band representations are Site symm. evals Γ evals K Invariants where ρ(C 2 ) = ρ 2 (C 4 ). Let us consider the site symmetry representation ρ(C 4 ) = e i π 4 σz . For C 4 , the band representations at HSPs are The matrix for the band representation of C 4 at Γ has eigenvalues e iπ/4 , e −iπ/4 , while the one at M has eigenvalues e 3iπ/4 , e −3iπ/4 . Thus, [M For C 2 , the band representations are always of the form ±e i π 2 σz , which has eigenvalues +i, −i, leading to vanishing symmetry indicators.
Let us now consider obstructions arising from the band representation when multiple orbitals localize at 1b. If the two orbitals have the same representation, e.g., ρ(C 4 ) = e i π 4 σz , the overall site symmetry representation, σ 0 e i π 4 σz , induces a band representation with invariant [M where N is the number of Kramers pairs in the site 3c. For C3, the band representation is constant across the C 3 invariant HSPs. Therefore, all invariants are trivial. For C 2 , the band representations at the HSPs are But since the site representation for a Kramers pair, ρ(C 2 ) = e i π 2 σz , has eigenvalues +i, −i, the band representations at Γ and M are the same.
The band representations at the HSPs are Thus, the invariants depend on the site symmetry representation ρ(C 3 ). They are shown in Table IX. Since TRS relates K with K , we only provide the representations at Γ and K. Note that, in order to have a trivial insulator with three movable Kramers pairs at 1b, two of them need to have the representation e i π 3 σz and the third one the representation −σ 0 . The band representations at the HSPs are So, just as for 1b, the invariants depend on the site symmetry representation ρ(C 3 ). They are shown in Table X.
Just as before, in order to have a trivial insulator with three movable Kramers pairs at 1c, two of them need to have the representation e i π 3 σz and the third one the representation −σ 0 .

C2 symmetry: Representations induced from any C2 invariant HSP
Since all C 2 invariants points are also TRIM points, and the C 2 eigenvalues of the site symmetry group of Kramers pairs is always +i, −i, which exhausts the representations, all the invariants due to C 2 are trivial.  Fig. 5 in the main text. '-' indicates that a given band representation cannot be written as a combination of EBRs.

Appendix H: Details of the ab-initio calculations
Fully relativistic DFT calculations were performed via the Vienna ab initio simulation package (VASP) 80,81 by employing the Perdew-Burke-Ernzerhof (PBE) 82,83 exchange-correlation functional and projected augmented-wave pseudopotentials 84,85 . For the self-consistent calculations, we used a 19 × 19 × 1 k-point grid generated for the Monkhorst-Pack method in case of Bi and Sb, and a 17 × 17 × 1 mesh for As. The plane wave basis cutoff was set to 400 eV (Bi and Sb) or 350 eV (As). A finer grid of 30 × 30 × 1 k-points was used later on in order to obtain the energy gaps and band representations. The lattice parameters in the equilibrium configuration, which are in good agreement with previous reports 64,86,87 , are summarized in Table XIII For open flake calculations, we employed the Siesta code 76 . We used pseudo-atomic orbitals (PAO) with a basis of double zeta plus polarization orbitals (DZP) and norm-conserving fully relativistic pseudopotentials from the PseudoDojo library 88 . The bulk crystal structure was terminated to obtain a hexagonal structure of 546 Sb atoms, and 30 Te atoms were added to the edges in order to passivate the edge states (as shown in Fig. 5 f). The distance between Te and edge Sb atoms was set to a value 3.02 A, which was determined from the structure relaxation of an armchair Sb ribbon with Te adatoms at the edge. The DFT data post-processing was performed with the sisl Python package 89 .
The irreducible representations of bands at high-symmetry points were obtained using the irrep code 90 , which relies on the double space group character tables 91 published on the Bilbao Crystallographic server 92 .