Glide-Mirror Protected First- and Second-Order Topological Crystalline Insulator

Xiaoting Zhou, ∗ Chuang-Han Hsu, Cheng-Yi Huang, 1 Mikel Iraola, 5 Juan L. Mañes, Maia G. Vergniory, 6 Hsin Lin, and Nicholas Kioussis † Department of Physics and Astronomy, California State University, Northridge, CA 91330, USA Department of Electrical and Computer Engineering, Faculty of Engineering, National University of Singapore, Singapore 117583 Institute of Physics, Academia Sinica, Taipei 11529, Taiwan Donostia International Physics Center, 20018 Donostia-San Sebastian, Spain Department of Condensed Matter Physics, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3, 48013 Bilbao, Spain (Dated: May 14, 2020)

Most topological insulators (TIs) discovered today in spinful systems can be transformed to topological semimetals (TSMs) with vanishing bulk gap via reduction of the spin-orbit coupling (SOC), which manifests the intrinsic links between the gapped TI phases and the gapless TSMs. Recently, we have proposed a novel family of TSMs in time-reversal invariant spinless systems, which host butterfly-like nodal-lines (NLs) consisting of a pair of identical concentric intersecting coplanar ellipses (CICE), protected by spatial symmetries, including a pair of glide symmetries. In addition, we have developed a minimal tight-binding model exhibiting CICE in space group (SG) P bam (No. 55), one of the nine identified feasible SGs. Here, we generalize the model by including the effect of SOC, and demonstrate that, with substantial SOC, the system undergoes a transition from the CICE TSM to a Z4 = 2 topological crystalline insulator (TCI). This TCI supports in turn a double-hourglass fermion on the (001) surface protected by two glide mirror symmetries, which originates from the intertwined drumhead surface states (DSSs) of the CICE NLs. Hence, the origin of the glide-symmetry-protected TCI is the CICE TSM induced by double-band-inversion (DBI). The higher order topology is further demonstrated for the first time by the emergence of one-dimensional (1D) helical hinge states protected by a glide symmetry.
Alongside the traditional Landau's approach [1] which generally characterizes the quantum states by the principle of spontaneous symmetry breaking, topological order has been, over the past two decades, established as an alternative principle in the classification of quantum states. However, the discovery of the quantum spin Hall effect (QSHE) [2][3][4] and topological insulators (TIs) [5][6][7] which are protected by intrinsic symmetries such as time-reversal symmetry (TRS), has indicated that symmetry plays a crucial role in classifying the topology of free fermion states. As it has been tabulated in a tenfold periodic table [8,9], a complete classification was achieved for the ten Altland-Zirnbauer symmetry classes [10]. This table specifies the possible topologically distinct ground states in free fermion systems in any dimensions under certain symmetries, such as, TRS, particlehole and chiral symmetries.
Subsequently, the concept has been generalized to spatial symmetry in crystalline systems. For instance, the weak topological insulators [5,6] rely on the discrete translation symmetry of the crystal lattice, and the topological crystalline insulators (TCIs) [11] are protected by other space-group symmetries (Q), such as mirror [12] and rotational symmetries [13][14][15]. Without breaking the corresponding symmetries, the unique properties owned by the specific systems remain intact in the presence of any perturbations. Consequently, such type of systems are known to harbor symmetry-protected topological (SPT) phase. [16].
In free fermion systems, SPT insulators harbor a central paradigm referred to as the bulk-boundary correspondence [17]. A bulk with gapped excitations hosts anomalous gapless, topologically nontrivial boundary states in lower dimensions. For instance, a d-dimensional topological insulators (TI) exhibits (d − 1)-dimensional Dirac boundary states, irrespective of the orientation of the boundary or the lattice termination [5,18]. However, for a TCI, the emergence of gapless states on the (d − 1)-dimensional boundary is guaranteed only if the boundary is invariant under the designated crystalline symmetry Q [11,19].
In contrast to the gapped topological phase, a topological semimetal (TSM) has gapless bulk band struc-tures, which are characterized by the topologically robust band-crossings manifolds between occupied and unoccupied bands in momentum space. Among them, nodalline semimetals (NLSMs) [23,24], which harbor onedimensional (1D) nodal lines (NLs), possess the highest variability. NLSMs with NLs integrated in various configurations have been studied under the assumption that the spin-orbit coupling (SOC) is negligible or absent, e.g., a chain link [25][26][27], a Hopf link [26], and a knot [28]. However, from another perspective, it is intuitive to raise the question whether additional topology could be unearthed when these intricate degenerate links are gapped out by substantial SOC. It has been known that a Z 2 strong TI can be realized from NLSMs with a single nodal ring when SOC is included, which represents the intrinsic link between a gapped topological phase and a gapless TSM. Nevertheless, similar studies for NLSMs with complex NL configurations remain deficient.
Recently, we have proposed [29] a new type of NLSM in spinless systems, which hosts a butterfly-like NL consisting of a pair of concentric intersecting coplanar ellipses (CICE) residing on a plane in k space. The CICE are protected by time-reversal and various spatial symmetries, including two glide symmetries, and only 9 out of 230 space groups (SGs) are feasible to host them. Furthermore, for one of the SG, P bam (No.55), we have constructed a minimal tight-binding model which exhibits CICE. The objective of this work is to include the (P bam) symmetry-invariant SOC in the model Hamiltonian and to study the topological properties. We demonstrate that the CICE act as the origin of a TCI protected by two glide symmetries. With substantial SOC, the CICE become anticrossing, thus driving a phase transition from NLSM to a Z 4 = 2 TCI [19,[30][31][32], due to the fact that the CICE are essentially sets of NLs stemming from the DBI. We further uncover the coexistence of the wallpaper fermion [33] with the HOTI phase, which is demonstrated by the emergence of double-hourglass fermion on the (001) surface and 1D helical hinge states when the sample geometries are distinctively and properly selected. Moreover, its higher order feature is demonstrated for the first time by 1D hinge states protected by glide symmetry.
Lattice Model-The CICE can be sustained by two glide mirror symmetries and only nine space groups (SGs) are feasible to host it [29] . The minimal 4-band tight-binding model for the CICE NLSM in SG P bam (No. 55) is a bipartite lattice, consisting of two sublattices denoted by A (gray) and B (blue) which occupy the 2a Wyckoff position at r A = (0, 0, 0) and r B = ( 1 2 , 1 2 , 0) in the unit cell (see Fig. 1 (a) for the structure). There are two orbitals, p z and d xy , for each sublattice, described by the Pauli matrix σ and τ for the A and B sublattices, respectively. The SG P bam can be generated by the mirror symmetry, M z = {m 001 |000}, and the two glide-mirror symmetries, G x = {m 100 | 1 where α, β, γ, δ 0 , and λ ij are constants. As discussed in detail in Ref. [29], the CICE emerge on the mirror plane (gray shaded area in Fig. 1(b)) centered at the high symmetry k point S = (π, π, 0) [R = (π, π, π)] ( Fig. 1(b)), under the conditions The various terms in Eq.
(1) describe the pair of concentric elliptic NLs, the NL anisotropy, and the angle between the NLs (see details in [29]). Since the CICE are composed of two NLs, it is anticipated to observe a pair of DSS [23] intertwined on the (001) surface. Topological Crystalline Insulator-In this section we consider the effect of SOC. The minimal SOC Hamiltonian, H SOC (k), to gap out the CICE-NL is of the form, where Γ ijk = τ i σ j s k (i, j, k ∈ {0, 1, 2, 3}), s 1,2,3 are Pauli matrices operating in spin space, and ζ ijk are constants. The band structure of H 0 and H = H 0 +H SOC are shown in Fig. 1(c) by the red and blue curves, respectively. In the presence of SOC, the CICE-NL TSM evolves into an insulating phase. The parameters are tuned to allow the system to host a single CICE centered at the S point in the absence of SOC and to have no additional band inversions at other time-reversal-invariant momentum points (TRIM) including SOC. In addition to {αδ S < 0 ∩ αβ > 0 ∩ α = β}, either condition below should be satisfied, where λ ± = λ 10 ± λ 13 and δ Γ,Z = δ 0 + (α + β ± γ) In order to determine its band topology, we implement the symmetry-indicator theory [19,30,32]. Crystals in the SG P bam are characterized by four symmetry indicators (SIs) [19,30,32], three Z 2 weak TI indices and one Z 4 index. The Z 4 index is defined as, where n + K (n − K ) is the number of occupied bands with parity + (−) at the TRIM points K. Due to the nonsymmorphic symmetries, bands at the TRIM points X, U, Y, T, S and R are four-fold degenerate. Besides, since inversion I anticommutes with the glide, G, or screw, S, symmetry operations (here S includes S y = 2 , at X, U, Y, and T, the parity of each four-fold degenerate state must be (+, +, −, −), which does not contribute to Z 4 . By enumerating the parity of the states at other TRIM points, we obtain Z 2,2,2,4 = (0, 0, 0, 2), corresponding to eight possible topological states [32]. To further narrow down the possible phases, we have calculated the mirror Chern numbers of M z and C m 001 0,π (m 001 0,π denotes the k z = 0, π mirror planes), following the method implemented in Ref. [15]. We find that C m 001 0,π = (2, 0). The corresponding Dirac surface states on the (010) surface are shown in Fig. 2(d), where the relevant k points and the schematic locations of the Dirac cones are illustrated in Fig. 2(a). Therefore, given that C m 001 0,π = (2, 0), there are two possible topological phases, which are listed in Table I. The first one is S-protected TCI with ν 2 [010] ) . Note that the nontrivial characteristics of the bands agree with the analysis from the elementary band representations (EBRs) [34]. The physical EBRs for the 2a Wyckoff position [35][36][37] require that the parity of Γ, Z, S and R has the same sign. However, because of the double band inversion at S guaranteed by the CICE-NL, both valence and conduction bands violate the physical EBRs, suggesting the emergence of nontrivial topology.
The (001) surface bands are shown in Fig. 2(c), where the (001) surface BZ and the corresponding high symmetry k points are displayed Fig. 2(a). Interestingly, the calculations reveal the emergence of nontrivial surface states aroundS which are composed of two intertwined surface Dirac cones, shown schematically in Fig. 2(b). We refer to these topological surface sates (TSSs) as doublehourglass fermions, which are of a particular type of wallpaper fermions [33]. Meanwhile, one can also obtain the double-hourglass fermions by doubling either the type-I or type-II intertwined DSSs of CICE TSM [29] with degenerate dispersions alongS −X (Ȳ ) guaranteed by G. Consequently, the CICE-NL induced topological phase belongs to the G-protected TCI. In the following, we provide more physics insights on this TCI phase and the double-hourglass fermions. In general, for time reversal symmetric systems, strong topological insulators (STIs) with a single band inversion at one TRIM point can be regarded as an elementary building block of the nontrivial insulating phase [14,19]. For each STI, the topological surface states of the (001) surface can be described by the Hamiltonian, h k = k x s 2 − k y s 1 [19]. The gapless feature of h k is protected by the TRS operator, T = −is 2 K, where K is the complex conjugation operator. For the current case, since there is a DBI at theS point, the induced TCI phase can be viewed as two copies of STIs. Accordingly, for the TSSs of the TCI, the only allowed TR invariant mass term takes the form, where µ 1,2,3 are the Pauli matrices acting on the two copies of h k , and m is constant. If M can be prohibited by any spatial symmetry Q, the anomalous gapless surface states will persist, indicating that the existing [011] [100] [010] [001] [011] [011] topology is protected by Q, which can be G x and G y , as derived below. AtS, the eigenvalues of G x and G y are ±1. To preserve TRS, the only available representations are G x = µ 2 ⊗ s 1 and G y = µ 2 ⊗ s 2 , which in turn lead to the rotational symmetry about the z-axis C 2z (= G x × G y ) = −iµ 0 ⊗ s 3 . Obviously, M cannot survive with G x and G y , but is allowed by C 2z . Consequently, there exist representations for G x and G y to support the (001) TSSs atS, which is the double-hourglass fermion shown in Fig. 2(b) and (c). The double-hourglass TSSs, shown in 2(b), can be described by the k · p Hamiltonian H T SS (q) = g 0 (q x s 2 − q y s 1 ) where all g's and a's are real constants.
Higher Order Topological Insulator-In addition to the topological surface states belonging to the (d − 1) bulkedge correspondence, i.e., the double-hourglass fermion and the M z -protected surface states, the G x , G y and M z symmetries can give rise to higher order (d−2) bulk-edge correspondence.
We have considered the nanorod geometry, shown in Fig. 3(a), with open boundary conditions along the [011] and [011] directions, respectively, and periodic boundary conditions along [100]. We find that the nanorod can support two pairs of hinges modes along the intersection lines between the (011) and (011) surfaces and the (011) and (011) surfaces, respectively. None of the above surfaces hosts gapless surface states, since there is no TCI phase supporting them. As discussed in Ref. [22], the hinge modes formed by the intersection of the (011) and (011) facets are generated via bending the (001) surface along [001] direction, which, however, preserves the G y symmetry for the hinges and the entire crystal. The original pair of Dirac cones forming the double hourglass fermion on the (001) surface become gapped with opposite mass terms and reside on the (011) and (011) surfaces, respectively (red and black massive Dirac cones in Fig. 3(a)). Hence, these two facets can be can be regarded as a combination of a 2D TI and an atomic insulator, where a pair of helical states, namely the hinge modes, emerge at the boundary. Similarly, the hinge modes along the intersection of the (011) and (011) facets are generated via bending the (010) surface along the [001] direction, which preserve the M z symmetry. Fig. 3(b) shows the band structure of the nanorod along the k x symmetry direction where the two pair of hinge modes are denoted in red and cyan, respectively. The distribution in real space of the two types of hinge modes along the [100] direction are displayed in Fig. 3(d) with the corresponding colors. For the purpose of clarity, on-site potentials V = 0.2 eV are added on the hinges formed by the {(011), (011)} and {(011), (011)} facets to better differentiate the two hinge states in energy. Notice that due to the symmetry constraint imposed by G y , there are two types of band connectivities in 1D [38], the hourglass type and the analogue of the quantum spin Hall (see Fig. 3(c)). It is interesting that the hourglass type, top panel of Fig. 3(c), is forbidden for hinge states; otherwise, the gapless feature of the 1D hinge states cannot be maintained. Therefore, the band connectivity naturally chooses the analogue of the quantum spin Hall [38], displayed in the bottom panel of Fig. 3(c).
Conclusion-In summary, inclusion of SOC in the model Hamiltonian describing our recently proposed novel family of butterfly-like CICE NDLs in SG P bam, [29] unveils intriguing topological properties. The SOC drives the TSM to a Z 4 = 2 topological crystalline insulator (TCI) with higher order topology, supporting in turn a double-hourglass fermion on the (001) surface protected by two coexisting glide mirrors. Moreover, its higher order topology is corroborated for the first time by the emergence of 1D hinge states protected by glide symmetry. The hamiltonian given by Eqs. (1)(2) in the Main Text can be obtained as a tight-binding model based on the Wyckoff position 2a with site-symmetry group 2/m. Two orbitals (p z , d xy ) are placed at the origin of the cell r A = (0, 0, 0) and at r B = (1/2, 1/2, 0). The Hamiltonian H(k) = H 0 (k) + H SOC (k) is invariant under the generators of SG55 given by where τ 0 , σ 0 and s 0 are 2×2 identity matrices. The matrices τ i act on the p z and d xy orbitals, σ i on the two sublattices and s i on electron spin. Note that the invariance of the hamiltonian under a space-group operation g = {R|v} means that Similarly, the invariance of the hamiltonian under the anti-unitary time-reversal operation T = T K can be expressed as where T = −iτ 0 ⊗ σ 0 ⊗ s 2 and K is ordinary complex conjugation.

FROM ATOMIC ORBITALS TO EBRS
Roughly speaking, a band representation [S1, S2] can be understood as the collection of electronic bands generated by placing a set orbitals at different positions in the unit cell, in such a way that the set is closed under all the symmetries of the crystal. A band representation obtained by placing orbitals at a maximal Wickoff position w that transform according to an irreducible representation τ of the site symmetry group G w is an elementary band representation (EBR) and is denoted by τ ↑ G| w . Any band representation that is not elementary can be written as a sum of EBRs [S3].
As our system respects time reversal symmetry (TRS), we are interested only in real and physically irreducible representations [S4]. Following the Bilbao Crystallographic Server [S5, S6] (BCS) conventions, a physically irreducible representation is denoted as the sum of two irreducible representations, which can be the same or different depending on their reality type [S4], with the + sign between the components of the pair omitted. For instance, Γ + 1 in Eq. (S5) below is a real irreducible representation whereas R + 1 and R + 2 are complex conjugate irreducible representations that together constitute the real, physically irreducible representation R + 1 R + 2 . In order to find the EBR contents of the model it is enough to know the transformation properties of the atomic orbitals s and p z . The site-symmery group 2/m has four single-valued real irreducible representations (irreps), A g , A u , B g , B u , and two double-valued physically irreducible irreps, 1Ē u . The spinless orbitals d xy and p z belong the single-valued irreps A g and A u respectively. As a consequence, in the absence of SOC the spectrum of the TB hamiltonian (5) is described by the sum of two EBRs induced from the 2a Wyckoff position, namely Then the application BANDREP at the BCS gives all the irreps at the different points in the Brillouin Zone (BZ). For the high symmetry TRIM points the result is where the numbers in parentheses give the dimensions of the irreps, which coincide with the degeneracies of the corresponding bands. This means, for instance, that the hamiltonian for spinless electrons in Eq. (5) of the paper necessarily has four non-degenerate bands at the Γ-point labelled by the 1-dim irreps {Γ + 1 (1), Γ + 2 (1), Γ − 1 (1), Γ − 2 (1)}, two doubly-degenerate bands at the R-point labelled by the 2-dim physically irreps R + 1 R + 2 (2) and R − 1 R − 2 (2), etc.
Note that in order to obtain this information one does not have to diagonalize the hamiltonian H 0 (k), as the result depends only on the orbital contents of the model. Note also that the degeneracies of the bands at the high symmetry points are independent of the parameters of the model, barring accidental degeneracies if some parameters are set to zero or otherwise fine-tuned. Electron spin and SOC can be easily incorporated by noting that spin transforms according to the physically irrep g of the site-symmetry group 2/m [S4]. Then, the double-valued irreps for s and p orbitals are obtained by taking the products of the respective single-valued irreps with the spin representation and this implies that the spectrum of the total hamiltonian H(k) = H 0 (k) + H SOC (k) is described by the sum of two EBRs Then BANDREP immediately gives the irrep contents and the structure of the spectrum for the total hamiltonian H(k). For the high symmetry TRIM points the result is

FROM EBRS TO SYMMETRY-BASED INDICATORS
The irreps at all the high symmetry points in the Brillouin zone (BZ) are tabulated at the BCS. For crystals with SOC and time reversal symmetry (TRS), we are interested in double-valued (because of the SOC), physically irreducible (because of TRS) representations. For SG55 there are 12 double-valued, physically irreducible representations Then any band representation can be characterized by a vector that tells how many times each physically irreducible representation appears at every high symmetry point. For instance, for the two EBRs in Eq. (S8) the vectors are and for the bands of the total hamiltonian H(k) The topology of an isolated subset of bands can be characterized by a set of symmetry-based indicators [S7]. In order to determine the indicators, one tries to write the vector giving the irreps for the subset of bands as a linear combination of all the vectors for the EBRs of the space group. There are several possibilities: 1. If the vector for the subset of bands can be written as a linear combination with positive integer coefficients, then all symmetry-based indicators (SI) for the subset vanish. The subset may still be topologically non-trivial, but all the eigenvalue-based topological invariants will be zero. In order to ascertain that the subset is non-trivial, one should use other tools, such as Wilson loops.
2. If the vector for the subset of bands can be written as a linear combination with integer coefficients, but some of the coefficients are necessarily negative, then the subset exhibits non-trivial fragile topology.
3. If the vector for the subset of bands can not be written as a linear combination with integer coefficients, then there are SI of the form Z n1 × Z n2 × . . ., where the Z n -factors reflect the existence of fractional coefficients in the linear combination of EBRs. These subsets of bands are usually listed as NLC, meaning non-linear combination of EBRs with integer coefficients, and are characterized by the existence of different topological invariants, such as Z 2 , Z 4 , etc.
4. If the vector for the subset of bands can not be written as a linear combination at all, then the subset is not really isolated, i.e., the gap between the subset and the other bands in the crystal closes at some point in the Brillouin zone.
The group of symmetry-based indicators for a given space group can be determined by diagonalization of the matrix EBR. This matrix is constructed by taking as columns the vectors of all the EBRs for the group. A look at BANDREP shows that, with TRS, there are eight double-valued EBRs for SG55. The first two have been given in Eq. (S10). Taking all 8 vectors as columns gives the EBR matrix for SG55  1 1 0 1 0  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  1 1 1 1 1 1 1 1  2 0 0 2 2 0 0 2  0 2 2 0 0 2 (S12) Note that EBR is a 12 × 8 non-square matrix, where 8 is the number of EBRs and 12 the number of irreps. Then, for any isolated subset of bands characterized by a 12-component vector B = (B 1 , . . . , B 12 ) giving the numbers of irreps, we will have to solve the linear set of equations where the solution X = (X 1 , . . . , X 8 ) gives the coefficients of the linear combination of the 8 EBRs. The presence of non-integer the coefficients in {X 1 , . . . , X 8 } implies the existence of symmetry-based indicators.
The possible existence of non-integer solutions and hence the existence of SI and non-trivial topological indicators may be predicted by diagonalizing the EBR matrix. As EBR is a non-square matrix with integer coefficients, one has to use a special diagonalization procedure known as Smith decomposition: where L and R are unitary integer square matrices of dimensions 12 × 12 and 8 × 8 respectively, and D is the following The number of non-zero diagonal elements is the rank of the group, equal to five in the case of SG55 with SOC and TRS. One can show that all the eigenvalues must be positive integers. If some of the eigenvalues{n 1 , n 2 , · · · } are greater than one, then the group of SI is Z n1 × Z n2 × . . .. In this case the SI group is Z 2 × Z 4 , in agreement with Table III of Ref. S7. Note that, upon inversion, the existence of greater than one eigenvalues may give rise to non-integer values for some of the coefficients of the solution X, showing that we are in case 3 of the previous discussion.

COMPUTING THE SYMMETRY INDICATORS FOR THE TB MODEL
According to Table IV of Ref. S7, SG55 has no SI in the absence of SOC and, as a consequence, can not have eigenvalue-based topological indices. For that reason we will consider the complete hamiltonian H(k) = H 0 (k) + H SOC (k). On the other hand, Eq. (S7) shows that the total eight-band spectrum can be written as the sum of two EBRs. Thus, in order to find non-trivial topology, we have to consider the possibility of splitting the eight bands into disconnected subsets, in such a way that the subsets have non-trivial SI. Put differently, if the two EBRs in Eq. (S8) represent the valence and conduction bands, both will be topologically trivial. We may get non-trivial topology only through inversions of energy levels between the conduction and valence bands.
A look at Eq. (S8) shows that inversions can take place only at the Γ, R, S and Z points. There are obviously three ways to split the irreps among the conduction and valence bands at the Γ and Z points, two at R and S and just one at the remaining high symmetry points. This gives a total of 3 × 3 × 2 × 2 = 36 patterns. The possible irrep contents for the different subsets of four bands, together with the computed topological invariants, are shown in Table S1 using the vector notation explained below Eq. (S9). The conduction band for phase considered in the Main Text with a single double band inversion at the S-point corresponds to case 4 in the Table S1.
Note that only half the possible subsets have been listed. The remaining subsets are complementary to the ones given there, and so are their SI. Considering for instance the third entry in the table, given by B = (2, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1), the complementary band is obtained by subtracting it from the complete set of bands of the model, as given in Eq. (S11) (2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2) − (2, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1) = (0, 2, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1) (S16) and its SI are Adding the topological invariants of the two complementary bands and remembering that Z n indicators are defined mod n, gives for the total 8-band system as expected. Note that the conduction and valence bands are always complementary. We now comment briefly on the method used to compute the SI in Table S1. The strong Z 2 index is given by the well known Fu-Kane formula [S8] where n − K is the number of bands with parity − at the TRIM point K. For SG55, the set of TRIM vectors coincides with the eight high symmetry points. This formula is easily evaluated with the help of Table S2, which gives the number of parity eigenstates for the twelve physically irreducible representations. The weak topological index Z 2w,1 (Z 2w,2 , Z 2w,2 ) is obtained by restricting the sum in Eq. (S19) to TRIMs with k x = π (k y = π, k z = π). Finally, the Z 4 index is given by the Fu-Kane-like formula [S9] Note that 14 out of the 18 sets of bands in Table S1 have SI and, as a consequence, their topology is necessarily non-trivial. This gives a very rich phase space to explore. On the other hand, the irrep content at the high symmetry points in cases 1, 5, 8 and 10 coincides with the EBR indicated in the last column. The corresponding subsets of bands might therefore be trivial according to our previous discussion. The EBRs for the complementary bands are obtained by changing the irrep label from g to u.

PHASE DIAGRAM FOR THE TB MODEL
The 18 cases in Table S1 reflect the possible orderings of the band energies at the high symmetry points in the BZ. The hamiltonian can be diagonalized analytically at all the high symmetry points. Note that, regarding possible orderings, points T , U , X and Y are irrelevant, as there is only one type of irrep in each case. It turns out that only 7 out of the 15 parameters in the Hamiltonian given by Eqs. (1)(2) in the Main Text enter the formulae for the energies at the relevant high symmetry point. More concretely, the spectra at the relevant high symmetry points can be written in terms of the five independent combinations δ 0 , δ ± = α + β ± γ , ξ ± = (λ 10 ± λ 13 ) 2 + ζ 2 233 (S21) according to Table S3. Then, as long as Z 4 is even, it is easy to choose the parameters in such a way that the desired ordering is obtained and the valence and conduction bands are separated by a gap throughout the BZ. For Z = 1 or 3, the bands must cross at a high symmetry line or plane and the phase is semimetallic [S9].
Three qualitatively different examples of bands for the TB model in the Main Text are presented in Fig. S1. The top bands correspond to the phase studied in the main text, where there is an inversion only at point S. This is a TCI, with the irrep contents and SI given in case 4 of Table S1. The middle bands correspond to case 12 in Table S1 and Z 4 = 3; there is a band crossing at the T Z line and the phase is semimetallic. In the last example, at the bottom of Fig. S1, all the SI vanish and the irrep contents of the valence and conduction bands coincide with the EBRs u ↑ G| 2a respectively. There are no band inversions and the system is in an atomic limit.             Table S1. Middle: Band inversions at the R, S and Z points, Z4 = 3, case 12. This is a semimetallic phase, note the band crossing at the T Z line. Bottom: This phase has no SI and the irrep contents of the valence and conduction bands are those of EBRs, case 1 in Table S1.