Wallpaper Fermions and the Nonsymmorphic Dirac Insulator

TL;DR

This paper introduces a new class of topological insulators called nonsymmorphic Dirac insulators, characterized by unique surface states protected by wallpaper group symmetries, and predicts their realization in specific materials.

Contribution

It identifies new topological crystalline phases protected by wallpaper group symmetries and predicts their realization in materials like Sr₂Pb₃ and Au₂Y₃.

Findings

Discovery of a new topological phase with a single Dirac fermion surface state.

Prediction of materials hosting nonsymmorphic Dirac fermions.

Development of topological invariants for orthorhombic space groups.

Abstract

Recent developments in the relationship between bulk topology and surface crystal symmetry have led to the discovery of materials whose gapless surface states are protected by crystal symmetries. In fact, there exists only a very limited set of possible surface crystal symmetries, captured by the 17 "wallpaper groups." We show that a consideration of symmetry-allowed band degeneracies in the wallpaper groups can be used to understand previous topological crystalline insulators, as well as to predict new examples. In particular, the two wallpaper groups with multiple glide lines, $pgg$ and $p4g$, allow for a new topological insulating phase, whose surface spectrum consists of only a single, fourfold-degenerate, true Dirac fermion. Like the surface state of a conventional topological insulator, the surface Dirac fermion in this "nonsymmorphic Dirac insulator" provides a theoretical exception to a fermion doubling theorem. Unlike the surface state of a conventional topological insulator, it can be gapped into topologically distinct surface regions while keeping time-reversal symmetry, allowing for networks of topological surface quantum spin Hall domain walls. We report the theoretical discovery of new topological crystalline phases in the A$_2$B$_3$ family of materials in SG 127, finding that Sr$_2$Pb$_3$ hosts this new topological surface Dirac fermion. Furthermore, (100)-strained Au$_2$Y$_3$ and Hg$_2$Sr$_3$ host related topological surface hourglass fermions. We also report the presence of this new topological hourglass phase in Ba$_5$In$_2$Sb$_6$ in SG 55. For orthorhombic space groups with two glides, we catalog all possible bulk topological phases by a consideration of the allowed non-abelian Wilson loop connectivities, and we develop topological invariants for these systems. Finally, we show how in a particular limit, these crystalline phases reduce to copies of the SSH model.

Tables

Table 1. Table 1: Materials in SG 55 identified by DFT to be gapless at the Fermi energy or topologically trivial by glide indices.

Glide-Trivial Tested Materials in SG 55 (22 Materials)
Chemical Formula	ICSD Number		Chemical Formula	ICSD Number		Chemical Formula	ICSD Number
Sc2Pt3Si2	247425		Y5Pd2In4	165133		Y2RuB6	427249, 603493,
							615404
Na3Cu4S4	10004		NaIn3S5	426706		AgHg2PO4	2208
Al3Pt5	55579, 58135,		CuHgSeCl	16450		La2InSi2	193223
	656681
La2SnS5	2313, 641853		Ca2PbO4	36629		Ca2SnO4	9011, 173626,
							187742
Cd2SnO4	69296, 9010		Cs2Te2	83351		Mg2Ru5B4	61039
Sr5Sn2P6	63593		Ta2Pd3Se8	73318		Ca5Sn2As6	61037
Rb2Te2	83350		Ta4SiTe4	40207, 659266		Tl2PdSe2	79601
YCrB4	16171

Table 2. Table 2: Materials in SG 125 identified by DFT to be gapless at the Fermi energy or topologically trivial by glide indices.

Glide-Trivial Tested Materials in SG 125 (18 Materials)
Chemical Formula	ICSD Number		Chemical Formula	ICSD Number		Chemical Formula	ICSD Number
SrNa2P4O12	37171		La2NiGa12	161765		K3AgSn3Se8	416330
K3NaSn3Se8	280286		Ge4N6Sr11	170982		La2PdGa12	171486, 183717
K2Sr(VO3)4	155420, 250103		BaAg2Hg2O4	40835		BaAl2Te4	41165
CeMn2Ge4O12	50695		Cs2Sr(VO3)4	250105		KCeSe4	67656
Ag2Ca(VO3)4	161369		Ag2Sr(VO3)4	161371		Ce2CuGa12	161767
La2NiGa12	161765		Na2Sr(VO3)4	155419, 161370,		RbAg5Se3	50738
				250102

Table 3. Table 3: Materials in SG 127 identified by DFT to be gapless at the Fermi energy or topologically trivial by glide indices.

Glide-Trivial Tested Materials in SG 127 (57 Materials)
Chemical Formula	ICSD Number		Chemical Formula	ICSD Number		Chemical Formula	ICSD Number
CaCu9Cd2	424134		YB2C2	427155		Sc2MgGa2	260213
NaNbO3	192406, 280100,		NaMgF3	193088		MgY2Ge2	423457
	23563, 236892
BaPtLa2O5	68794		InNi2Sc2	107333		CsSnI3	69995, 262925
CdY2Ge2	414169		B2Nb3Ru5	423469		AlPt3	107439, 609126,
							609153, 656679
Y2Cu2Mg	411711, 419472		Sc2Ni2Sn	54348		Pb2Br2CO3	250396, 29114
Na2Bi5AuO11	74365, 164986		Mn2Ga5	249632, 634613,		La7Ni2Zn	159116
				634639
La2Ni2Mg	107327		La2Cu2In	411708, 628002		K3Li2(NbO3)5	164890
K2LaTa5O15	421750		Co2Zr2In	107331		CeB4	417745, 24682,
							612705
CoIn3	102501, 623922		Ba4In2Te2Se5	425588		BaHg2O2Cl2	77509
Ba3Nb5O15	69993		Ba3Ta5O15	79810		CeB2C2	88560, 94036
							280180, 88857
Ca2Au2Pb	409531		Cs3GeF7	202917		Hg2PbI2S2	59204
KAlF4	16413, 60524,		K2CsPdF5	72301		KCuF3	21110
	201947
KMo4O6	68533		LaB2C2	94035		La2InGe2	87511
Li2Sn5	26200		Mg2SiIr5B2	69487		NaMo4O6	40962
SnMo4O6	92839		V3B2	88317, 107318,		La2Pd2Pb	99190
				615662
Y2Pd2Pb	99189		Y2BaPdO5	202819		NaTaO3	23322, 88377
							239692, 280101
La2Ni2Mg	107327		La2Cu2Mg	411709		La2Ni5C3	62277, 67376
Ta4	54204		Ti2In5	401730		Au2InY2	658835
Au2SnTb2	658834		Ta3Ga2	107309, 635467,		Tl2GeTe5	69035
				635490

Equations

$${\left[{\mathcal{W}}_{(k_x,k_y,k_{z0})}\right]}_{ij}$$

$$\hat{\Pi }(k_x,k_y,k_{z0})|u^j(k_x,k_y,k_{z0})⟩,$$

$$\displaystyle\hat{\Pi}(k_{x},k_{y},k_{z})\equiv\$$

$$\hat{\mathcal{P}}(k_x,k_y,k_z+$$

$${\chi }_y$$

$${\mathcal{H}}_{\bar{M}}={\tau }^x\left(v_x{\sigma }^x{k}_x+v_y{\sigma }^y{k}_y\right),$$

$$\mathcal{H}=v_x{k}_x{\sigma }^x+v_y{k}_y{\sigma }^y.$$

$$C=\frac{1}{2\pi i}\log \left(P\exp i{\oint }_{FS}\mathcal{A}\cdot d\mathbf{k}\right),$$

$$\mathcal{H}=v_x{k}_x{\tau }^x{\sigma }^y+k_y[v_{y1}{\tau }^y+v_{y2}{\tau }^x{\sigma }^x+v_{y3}{\tau }^x{\sigma }^z],$$

$$g_{x,y}=\{m_{x,y}|\frac{1}{2}\frac{1}{2}\},$$

$$g_x=\{m_x|\frac{1}{2}\frac{1}{2}\},C_{4z}=\{C_{4z}|00\}.$$

$$g_x(\mathcal{T}{g}_y|+⟩)=-\mathcal{T}{g}_y{g}_x{t}_x{t}_{-y}|+⟩=-\mathcal{T}{g}_{y}i{e}^{-i{k}_y/2}{e}^{i{k}_x}|+⟩=-i{e}^{i{k}_y/2}\left(\mathcal{T}{g}_y|+⟩\right).$$

$${\mathcal{H}}_{\bar{M}}={\tau }^x(v_x{\sigma }^x{k}_x+v_y{\sigma }^y{k}_y).$$

$$|{\varphi }_{\mathbf{k},\alpha }⟩=\frac{1}{\sqrt{N}}\sum _{\mathbf{R}}{e}^{i\mathbf{k}\cdot (\mathbf{R}+{\mathbf{r}}_{\alpha })}|{\varphi }_{\mathbf{R},\alpha }⟩.$$

$$H(\mathbf{k}{)}_{\alpha ,\beta }=⟨{\varphi }_{\mathbf{k},\alpha }|\hat{H}|{\varphi }_{\mathbf{k},\beta }⟩,$$

$$\hat{H}(\mathbf{k}+\mathbf{G})=\hat{V}(\mathbf{G}{)}^{-1}\hat{H}(\mathbf{k})\hat{V}(\mathbf{G}),$$

$$|u^n(\mathbf{k}+\mathbf{G})⟩=\hat{V}(\mathbf{G}{)}^{-1}|u^n(\mathbf{k})⟩.$$

$$\hat{g}|{\varphi }_{\mathbf{R},\alpha }⟩=|{\varphi }_{{\mathbf{R}}^{\prime },\beta }⟩{\left[U_g\right]}_{\alpha \beta },$$

$${\hat{g}}_{\mathbf{k}}|{\varphi }_{\mathbf{k},\alpha }⟩$$

$$=\frac{1}{\sqrt{N}}\sum _{\mathbf{R}}{e}^{i\mathbf{k}\cdot (\mathbf{R}+{\mathbf{r}}_{\alpha })}|{\varphi }_{D_g(\mathbf{R}+{\mathbf{r}}_{\alpha })+\delta -{\mathbf{r}}_{\beta },\beta }⟩{\left[U_g\right]}_{\alpha \beta }$$

$$=\frac{1}{\sqrt{N}}{e}^{-i(D_g\mathbf{k})\cdot \delta }\sum _{\mathbf{R}}{e}^{i(D_g\mathbf{k})\cdot (D_g(\mathbf{R}+{\mathbf{r}}_{\alpha })+\delta )}|{\varphi }_{D_g(\mathbf{R}+{\mathbf{r}}_{\alpha })+\delta -{\mathbf{r}}_{\beta },\beta }⟩{\left[U_g\right]}_{\beta \alpha }$$

$$=e^{-i(D_g\mathbf{k})\cdot \delta }|{\varphi }_{D_g\mathbf{k},\beta }⟩{\left[U_g\right]}_{\alpha \beta }.$$

$${\hat{g}}_{\mathbf{k}}=e^{-i(D_g\mathbf{k})\cdot \delta }{U}_g.$$

$$H(\mathbf{k})=g_{\mathbf{k}}^{†}H(D_g\mathbf{k})g_{\mathbf{k}}.$$

$${\mathcal{T}}_{g,\mathbf{k}}\equiv e^{i(D_g\mathbf{k})\cdot \delta }{U}_{\mathcal{T}g}\mathcal{K},$$

$$H(\mathbf{k})={\mathcal{T}}_{g,\mathbf{k}}^{-1}H(-D_g\mathbf{k}){\mathcal{T}}_{g,\mathbf{k}}.$$

$$\hat{\mathcal{P}}(\mathbf{k})=\sum _{n=1}^{n_{\mathrm{occ}}}|u^n(\mathbf{k})⟩⟨u^n(\mathbf{k})|,$$

$$\hat{\mathcal{P}}(\mathbf{k})=\hat{V}(\mathbf{G})\hat{\mathcal{P}}(\mathbf{k}+\mathbf{G})\hat{V}(\mathbf{G}{)}^{†}.$$

$$\hat{\mathcal{P}}(D_g\mathbf{k})=\sum _{n=1}^{n_{\mathrm{occ}}}{\hat{g}}_{\mathbf{k}}|u^n(\mathbf{k})⟩⟨u^n(\mathbf{k})|{\hat{g}}_{\mathbf{k}}^{†}={\hat{g}}_{\mathbf{k}}\hat{\mathcal{P}}(\mathbf{k}){\hat{g}}_{\mathbf{k}}^{†}.$$

$${\hat{g}}_{\mathbf{k}}\hat{V}(\mathbf{G})\hat{H}(\mathbf{k}+\mathbf{G})\hat{V}(\mathbf{G}{)}^{†}{\hat{g}}_{\mathbf{k}}^{-1}=\hat{V}(D_g\mathbf{G})\hat{H}(D_g(\mathbf{k}+\mathbf{G}))\hat{V}(D_g\mathbf{G}{)}^{†}.$$

$${\hat{g}}_{\mathbf{k}+\mathbf{G}}=\hat{V}(D_g\mathbf{G}{)}^{†}{\hat{g}}_{\mathbf{k}}\hat{V}(\mathbf{G}).$$

$$\hat{\mathcal{P}}(-D_g\mathbf{k})={\mathcal{T}}_{g,\mathbf{k}}\sum _{n=1}^{n_{\mathrm{occ}}}|u^n(\mathbf{k})⟩⟨u^n(\mathbf{k})|{\mathcal{T}}_{g,\mathbf{k}}^{†}$$

$${\mathcal{T}}_{g,\mathbf{k}+\mathbf{G}}=\hat{V}(-D_g\mathbf{G}{)}^{†}{\mathcal{T}}_{g,\mathbf{k}}\hat{V}(\mathbf{G}).$$

Full text

I Introduction

Topological phases stabilized by crystal symmetries have already proven to be both a theoretically and an experimentally rich set of systems. The first class of these proposed materials, rotation or mirror topological crystalline insulators (TCIs), host surface fermions protected by the projection of a bulk mirror plane or rotation axis onto a particular surfaceTeo et al. (2008); Fu (2011); Kim et al. (2015). They have been observed in SnTe Hsieh et al. (2012); Tanaka et al. (2012) and related compounds Dziawa et al. (2012); Xu et al. (2012). Recent efforts to expand this analysis to nonsymmorphic systems with surface glide mirrors – operations composed of a mirror and a half-lattice translation – have yielded additional exotic free fermion topological phases, which can exhibit the so-called surface gapless “hourglass fermions,” and the glide spin Hall effect Liu et al. (2014a); Wang et al. (2016); Alexandradinata et al. (2016); Shiozaki et al. (2016); Chang et al. (2017a). The theoretical proposal of Wang et al., 2016; Alexandradinata et al., 2016 has recently also seen incipient experimental support Ma et al. (2017).

In addition, topological insulators (TIs) – crystalline or otherwise – provide exceptions to fermion doubling theorems. These theorems impose fundamental bounds on phenomena in condensed matter physics. For example, in 2D, a single Kramers degeneracy in momentum space must always have another partner Kramers crossing elsewhere in the Brillouin Zone (BZ), otherwise the Berry phase of a loop enclosing the degeneracy suffers from ambiguity Haldane (2004). The discovery of the topological insulator provided the first exception to this theorem: in these systems 2D Kramers pairs are allowed to be isolated on a single 2D surface because they are connected to a 3D topological insulating bulk and have their partners on the opposite surface Fu and Kane (2007); Moore and Balents (2007).

Higher-fold-degenerate bulk fermions, such as Dirac points, which are stabilized by crystal symmetry, may come with their own fermion doubling theorems Young et al. (2012); Steinberg et al. (2014); Wieder et al. (2016); Bradlyn et al. (2016). As noted in Young and Kane, 2015, a single fourfold-degenerate Dirac fermion cannot be stabilized by 2D crystal symmetries as the only nodal feature at a given energy; it must always have at least one partner or accompanying hourglass nodal points. This is because a single Dirac point in 2D represents the quantum critical point (QCP) separating a trivial insulator (NI) from a topological insulator. Shown in more detail in Appendix A.2, stabilizing just one of these Dirac points with crystal symmetries would therefore force the broken-symmetry NI and TI phases to be related by just a unitary transformation, violating their $\mathbb{Z}_{2}$ topological inequivalence. In this manuscript, we report a new class of symmetry-protected topological materials which, like the topological insulator before it, circumvents this restriction by placing a single, stable Dirac point on the surface of a 3D material.

To realize this, the crucial requirement is that the surface preserves multiple nonsymmorphic crystal symmetries (Appendix A.2). Until now, most attention has been paid to crystal systems with surfaces that preserve only a single glide mirror. However, two of the 17 two-dimensional surface symmetry groups, called wallpaper groups, host two intersecting glide lines Conway et al. (2008). As we show in Appendix B, the algebra of the two glides requires that bands appear with fourfold degeneracy at a single time-reversal-invariant momentum (TRIM) at the edge of the BZ.

In this work, we study the non-interacting topological phases allowed in bulk crystals with surfaces invariant under the symmetries of these two wallpaper groups, $pgg$ and $p4g$. We show that, in addition to generalizations of the hourglass fermions introduced in Wang et al., 2016, they host a novel topological phase characterized by a single, symmetry-enforced fourfold Dirac surface fermion, i.e., twice the degeneracy of a traditional topological insulator surface state. This Dirac fermion is nonsymmorphic symmetry-pinned to the QCP between a TI and an NI, allowing for controllable topological phase transitions of the 2D surface under spin-independent glide-breaking strain.

We classify the allowed topological phases for orthorhombic, tetragonal, and cubic crystals with two perpendicular glides that are preserved on a single surface by considering the possible connectivities of the non-abelian Wilson loop eigenvalues Fu and Kane (2006); Ryu et al. (2010); Soluyanov and Vanderbilt (2011); Yu et al. (2011); Taherinejad et al. (2014); Alexandradinata et al. (2014a, 2016). We demonstrate that these systems allow for three classes of topological phases: an hourglass phase with broken $C_{4z}$ symmetry, a previously uncharacterized “double-glide spin Hall” phase, and the novel topological nonsymmorphic Dirac insulating phase mentioned above. We present topological invariants to distinguish these phases and use these invariants to predict material realizations of our new phases. Using density functional theory (DFT) to calculate bulk Wilson loop eigenvalues and surface Green’s functions, we report the existence of the nonsymmorphic Dirac insulating phase in Sr2Pb3 in space group (SG) 127. We also report the discovery of a related topological hourglass fermion phase in the narrow-gap material Ba5In2Sb6 in SG 55. In Appendix F.2, we show that two additional members of the SG 127 A2B3 family of materials, Au2Y3 and Hg2Sr3, can also realize this topological hourglass phase under $(100)$-strain, giving unique promise for strain-engineering topological phase transitions in these materials. Finally, we show in Appendix G how these crystalline phases reduce in a particular limit to copies of the Su-Schrieffer-Heeger (SSH) model Su et al. (1979).

II Wallpaper Groups $pgg$ and $p4g$

The surface of a nonmagnetic crystal is itself a lower-dimensional crystal, which preserves a subset of the bulk crystal symmetries, and all 2D nonmagnetic surfaces are geometrically constrained to be invariant under the action of the 17 wallpaper groups. The set of spatial wallpaper group symmetries is restricted to those 3D space group symmetries which preserve the surface normal vector: rotations about that vector and in-plane lattice translations, mirror reflections, and glide reflections. Surface band features will therefore be constrained by the irreducible corepresentations of these symmetries and their momentum-space compatibility relations Dong and Liu (2016); Bradlyn et al. (2017) (Appendix B). Focusing on the 17 wallpaper groups which describe the surfaces of 3D crystals with strong spin-orbit coupling (SOC) and time-reversal symmetry, we designate all topological 2D surface nodes formed from symmetry-enforced band degeneracies “wallpaper fermions.” In the language of group theory, wallpaper fermions are therefore fully captured by the irreducible corepresentations of the wallpaper groups, such that, for example, the fourfold-rotation-protected quadratic topological node introduced in Fu, 2011 is a “spinless wallpaper fermion,” whereas the magnetic twofold surface degeneracy in Fang and Fu, 2015, for which both bands have the same irreducible representations and are instead prevented from gapping by an additional rotation- and time-reversal-enforced ($C_{2}\times\mathcal{T}$) 1D loop invariant, is not a wallpaper fermion.

Though other works have focused on the mathematical classification of topological phases protected by wallpaper group symmetries Shiozaki et al. (2017); Song et al. (2017a), we here, in addition to providing a further topological classification, search for topological insulating phases with phenomenologically distinct surface states. For the symmorphic wallpaper groups, topological rearrangements of the minimal surface band connectivities, recently enumerated for 3D bulk systems in Bradlyn et al., 2017, allow for only quantum spin Hall (QSH) phases and rotational variants of the mirror TCI phases, which exhibit twofold-degenerate free fermions along high-symmetry lines in the surface BZ Hsieh et al. (2012); Tanaka et al. (2012); Dziawa et al. (2012); Xu et al. (2012).

For the four wallpaper groups with nonsymmorphic glide lines ($pg$, $pmg$, $pgg$, and $p4g$), this picture is enriched. Even in 2D, glide symmetries require that groups of four or more bands be connected, an effect which frequently manifests itself in hourglass-like band structures Young and Kane (2015); Watanabe et al. (2015); Wieder and Kane (2016). For the wallpaper groups with only a single glide line, $pg$ and $pmg$, surface bands can be connected in topologically-nontrivial patterns of interlocking hourglasses Wang et al. (2016); Alexandradinata et al. (2016), or in “Möbius strip” connectivities if time-reversal-symmetry is relaxed Shiozaki et al. (2015); Fang and Fu (2015), both of which exhibit twofold-degenerate surface fermions along some of the momentum-space glide lines. In the remaining two wallpaper groups with multiple glide symmetries, $pgg$ and $p4g$, higher-degenerate wallpaper fermions are uniquely allowed.

We consider a $z$-normal surface with glides $g_{x,y}\equiv\{m_{x,y}|\frac{1}{2}\frac{1}{2}0\}$, i.e., a mirror reflection through the $x,y$-axis followed by a translation of half a lattice vector in the $\hat{x}$ and $\hat{y}$ directions (Fig. 1). When spin-orbit coupling is present, $g_{x,y}^{2}=-e^{ik_{y,x}}$. At the corner point, $\bar{M}$ of the surface BZ ($k_{x}=k_{y}=\pi$), $g_{x}^{2}=g_{y}^{2}=+1$, and $\{g_{x},g_{y}\}=0$. When combined with time-reversal, $\mathcal{T}^{2}=-1$, this symmetry algebra requires that all states at $\bar{M}$ are fourfold-degenerate. Furthermore, wallpaper groups with two glides are the only nonmagnetic surface groups that admit this algebra, and therefore the only surface groups that can host protected fourfold degeneracies on strong-SOC crystals Wieder and Kane (2016). The examination of symmetry-allowed terms reveals that fourfold points in these wallpaper groups will be linearly-dispersing (Appendix B), rendering them true surface Dirac fermions, more closely related by symmetry algebra and quantum criticality to the bulk nodes in nonsymmorphic 3D Dirac semimetals Fu and Kane (2007); Young et al. (2012); Steinberg et al. (2014) than to the surface states of a conventional TI. In Appendix B, we provide proofs relating this algebra to Dirac degeneracy and dispersion.

For bulk insulators, the glide-preserving bulk topological phase and, consequently, $z$-normal surface states, can be determined without imposing a surface by classifying the allowed connectivities of the $z$-projection Wilson loop holonomy matrix Fidkowski et al. (2011); Alexandradinata et al. (2014a, 2016), a bulk quantity defined by:

[TABLE]

where we have defined the product of projectors,

[TABLE]

and $\hat{\mathcal{P}}(\mathbf{k})$ is the projector onto the occupied bands at $\mathbf{k}$. The rows and columns of $\mathcal{W}$ correspond to filled bands, where $|u^{j}(\mathbf{k})\rangle$ is the cell-periodic part of the Bloch wavefunction at momentum $k$ with band index $j$. The eigenvalues of $\mathcal{W}$ are gauge invariant and of the form $e^{i\theta(k_{x},k_{y})}$. As detailed in Appendices C and D, the Wilson bands inherit the symmetries of the $z$-normal wallpaper group, and therefore must also exhibit the required degeneracy multiplets of wallpaper groups $pgg$ and $p4g$. In particular, both surface and Wilson bands are twofold-degenerate along $\bar{X}\bar{M}$ ($\bar{Y}\bar{M}$) by the combination of time-reversal and $g_{y}$ ($g_{x}$) and meet linearly in fourfold degeneracies at $\bar{M}$.

By generalizing the $\mathbb{Z}_{4}$ invariant defined in Shiozaki et al., 2016 for the single-glide wallpaper groups Alexandradinata ; Xiong and Alexandradinata (2018), we define topological invariants for double-glide systems using the (001)-directed Wilson loop eigenvalues. For $g_{y}$ in a surface BZ in wallpaper group $pgg$, the quantized invariant $\chi_{y}$ is defined in Shiozaki et al., 2016 by integrating the Wilson phases, $\theta_{j}(k_{x},k_{y})$, along the path $\bar{M}\bar{Y}\bar{\Gamma}\bar{X}$:

[TABLE]

where $n_{\rm occ}$ is the number of occupied bands, the superscript $\pm$ indicates the glide sector, and the absence of a superscript indicates the line where $g_{y}$ is not a symmetry and the sum is over all bands. In the presence of an additional glide, $g_{x}$, one can obtain $\chi_{x}$ by the transformation $x\leftrightarrow y,\bar{X}\leftrightarrow\bar{Y}$ in Eq. (3). Though Eq. (3) appears complicated, $(\chi_{x},\chi_{y})$ can be easily evaluated by considering the bands within each glide sector which cross an arbitrary horizontal line in the Wilson spectrum (Appendix E.1). Wallpaper group $p4g$ also has $C_{4z}$ symmetry, which requires $\chi_{x}=\chi_{y}$ and implies the existence of the symmorphic mirrors, $\{m_{110}|\frac{1}{2}\frac{1}{2}0\}$ and $\{m_{1\bar{1}0}|\frac{1}{2}\bar{\frac{1}{2}}0\}$. These mirrors yield $\mathbb{Z}$ mirror Chern numbers, $n_{110},n_{1\bar{1}0}$, respectively, with values constrained by the glide invariant $\chi_{x}$, $(-1)^{n_{110}}=(-1)^{n_{1\bar{1}0}}=(-1)^{\chi_{x}}$ (Appendix E.4).

To enumerate the allowed topological phases shown in Fig. 2, we consider possible restrictions on $(\chi_{x},\chi_{y})$. Though $\chi_{x,y}$ can individually take on values $0,1,2,3$; only pairs that satisfy $\chi_{x}+\chi_{y}\mod 2=0$ are permitted in bulk insulators; this can be understood as follows: if $\chi_{x}+\chi_{y}$ is odd, the 2D surface consisting of the four partial planes $(0\leq k_{x}\leq\pi,0(\pi),k_{z})$ and $(0(\pi),0\leq k_{y}\leq\pi,k_{z})$ possesses an overall Chern number, which implies the existence of a gapless point Fang et al. (2016); Alexandradinata and Bernevig (2016), contradicting our original assumption that the system is insulating. We present a rigorous proof in Appendix E.2, and show that the remaining collection of eight insulating phases is indexed by the group $\mathbb{Z}_{4}\times\mathbb{Z}_{2}$.

For $\chi_{x,y}=1,3$, the system is a strong topological insulator (STI). These four “double-glide spin Hall” phases possess the usual twofold-degenerate Kramers pairs at $\bar{\Gamma},\bar{X},$ and $\bar{Y}$ as well as a fourfold-degenerate Dirac point at $\bar{M}$. The four STI phases are topologically distinct, but will appear indistinguishable in glide-unpolarized ARPES experiments. However, if two double-glide spin Hall systems with differing $\chi_{x,y}$ are coupled together, the resulting surface modes will distinguish between $\chi_{x,y}=1,3$ Shiozaki et al. (2016) (Appendix E.2).

When $\chi_{x,y}=0,2$, the system is in a topological crystalline phase. For $(\chi_{x},\chi_{y})=(0,2)$ or $(2,0)$, which is only permitted in a $C_{4z}$-broken surface $pgg$, a variant of the hourglass insulating phase Wang et al. (2016) is present on the surface. For example, when $(\chi_{x},\chi_{y})=(0,2)$, either time-reversed partners of twofold-degenerate free fermions live along $\bar{\Gamma}\bar{X}$ or both twofold-degenerate fermions live along $\bar{\Gamma}\bar{Y}$ and a fourfold-degenerate Dirac fermion exists at $\bar{M}$.

Finally, but most interestingly, for $\chi_{x}=\chi_{y}=2$, we find that the system exists in a previously uncharacterized “nonsymmorphic Dirac insulating” phase, capable of hosting just a single fourfold-degenerate Dirac surface fermion at $\bar{M}$.

III Materials Realizations

We apply DFT including the effects of SOC to predict topological phases stabilized by wallpaper groups $pgg$ and $p4g$ in previously synthesized, stable materials. The details of these large-scale calculations are provided in Appendix F. We find double-glide topological phases on the $(001)$-surface (wallpaper group $p4g$) of three members of the SG 127 ($P4/mbm$) A2B3 family of materials: Sr2Pb3 Merlo (1984); Bruzzone et al. (1981), Au2Y3 Chai and Corbett (2011), and Hg2Sr3 Druska et al. (1996); Gumiński (2005). Shown in Fig. 3, we find that Sr2Pb3 has at each crystal momentum a gap at the Fermi energy, in spite of the presence of electron and hole pockets. A Wilson loop calculation of the bands up to this gap (Fig. 3(d)) indicates that this material possesses the bulk topology of a $(2,2)$ nonsymmorphic Dirac insulator (Appendix F.2). Calculating the surface spectrum through surface Green’s functions (Fig. 4), we find that the $(001)$-surface of Sr2Pb3, while displaying an overall metallic character, develops a gap of 45 meV at the Fermi energy at $\bar{M}$. Inside this gap, we observe a single, well-isolated, fourfold-degenerate surface Dirac fermion.

Unlike Sr2Pb3, Au2Y3 and Hg2Sr3 in SG 127 ($P4/mbm$) are gapless, with bulk $C_{4z}$-protected Dirac nodes Wang et al. (2013) present near the Fermi energy. In Appendix F.2, we show that under weak $(100)$-strain, these Dirac nodes can be gapped to induce the $(0,2)$ topological hourglass phase in these two materials.

We additionally find that the $(001)$-surface (wallpaper group $pgg$) of the narrow-gap insulator Ba5In2Sb6 in SG 55 ($Pbam$) Cordier and Steher (1988) hosts a double-glide topological hourglass fermion. Shown in Fig. 5, we find that Ba5In2Sb6 develops an indirect band gap of 5 meV (direct band gap: 17 meV). The Wilson loop spectrum obtained from the occupied bands, shown in (Fig. 5(d)), demonstrates that this material is a $(2,0)$ double-glide topological hourglass insulator (Appendix F.1). We find that the $(001)$-surface of Ba5In2Sb6 has a projected insulating bulk gap which is spanned along $\bar{Y}\bar{\Gamma}$ by the top and bottom bands of two different topological hourglass fermions, which themselves are degenerate with states in the bulk spectrum (Fig. 6). However, these fermions are topologically connected to a clearly distinguishable hourglass fermion along $\bar{\Gamma}\bar{X}$ and a fourfold-degenerate surface Dirac fermion at $\bar{M}$, both of which could in-principle be observed through ARPES.

IV Discussion

We have demonstrated the existence of a nonsymmorphic Dirac insulator – a topological crystalline material with a single fourfold-degenerate surface Dirac point stabilized by two perpendicular glides. After an exhaustive study of the 17 time-reversal-symmetric, strong-SOC wallpaper groups, only $pgg$ and $p4g$ are revealed to be capable of supporting this fourfold fermion. This phase is one of eight topologically distinct phases that can exist in insulating orthorhombic crystals with surfaces that preserve two perpendicular glides; we have classified all eight phases by topological indices $(\chi_{x},\chi_{y})$ that characterize the connectivity of the $z$-projection Wilson loop spectrum. We report the discovery of the nonsymmorphic Dirac insulating phase in Sr2Pb3 and of related double-glide topological hourglass phases in Ba5In2Sb6, as well as in $(100)$-strained Au2Y3 and Hg2Sr3. We also report the theoretical prediction of a set of novel double-glide spin Hall phases. Though their surface Kramers pairs and fourfold Dirac fermions should be distinctive in ARPES experiments, a characterization of transport in the double-glide spin Hall phases remains an open question.

We also find that there exists a simple intuition for the topological crystalline phases $\chi_{x,y}=0,2$. In Appendix G.1, we present an eight-band tight-binding model which, when half-filled, can be tuned to realize all $\mathbb{Z}_{4}\times\mathbb{Z}_{2}$ double-glide insulating phases. In a particular regime of parameter space, in which SOC is absent at the $\bar{X}$ and $\bar{Y}$ points and bulk inversion symmetry is imposed, the Wilson loop eigenvalues at the edge TRIMs are pinned to $\pm 1$ ($\theta(\bar{M}/\bar{X}/\bar{Y})=0,\pi$) and each TRIM represents the end of a doubly-degenerate SSH model Su et al. (1979). In this limit, when the product of parity eigenvalues at $\bar{\Gamma}$ satisfies $\xi(\bar{\Gamma})=+1$, the bulk topology is fully characterized by the relative SSH polarizations, $\chi_{x,y}=2\{\frac{1}{\pi}[\theta(\bar{M})-\theta(\bar{Y},\bar{X})]\mod 2\}$.

Finally, as the $(2,2)$ topological surface Dirac point is symmetry-pinned to the QCP between a 2D TI and an NI, we examine its potential for hosting strain-engineered topological physics. Consider the two-site surface unit cell in wallpaper group $pgg$ from Fig. 1. In the $(2,2)$ nonsymmorphic Dirac insulating phase, the surface Dirac fermion can be captured by the $k\cdot p$ Hamiltonian near $\bar{M}$:

[TABLE]

where $\tau$ is a sublattice degree of freedom, $\sigma$ is a $\mathcal{T}$-odd orbital degree of freedom, and $g_{x/y}=\tau^{y}\sigma^{x/y}$ ($g_{x/y}^{2}=+1$). There exists a single, $\mathcal{T}$-even mass term, $V_{m}=m\tau^{z}$, which satisfies $\{\mathcal{H}_{\bar{M}},V_{m}\}=0$ and is therefore guaranteed to fully gap $\mathcal{H}_{\bar{M}}$. Therefore, surface regions with differing signs of $m$ will be in topologically distinct gapped phases and must be separated by 1D topological QSH surface domain walls, protected only by time-reversal symmetry Jackiw and Rebbi (1976). As $pgg$ has point group $2mm$, and $\{V_{m},g_{x,y}\}=0$, $V_{m}$ can be considered an $xy$ $A_{2}$ distortion Tinkham (1964), which could be achieved by strain in the $x+y$ direction and compression in the $x-y$ direction. These domain walls would appear qualitatively similar to those proposed in bilayer graphene Zhang et al. (2011); Ohta et al. (2006); Oostinga et al. (2007). However, whereas those domain walls are protected by sublattice symmetry and are therefore quite sensitive to disorder, domain walls originating from nonsymmorphic Dirac insulators are protected by only time-reversal symmetry, and therefore should be robust against surface disorder. Under the right interacting conditions or chemical modifications, a nonsymmorphic Dirac insulator surface may also reconstruct and self-induce regions of randomly distributed $\pm m$, separated by a network of 1D QSH domain walls. These domain walls are closely related to the 1D helical hinge modes of “second-order” topological insulators Schindler et al. (2018); Song et al. (2017b); Schindler et al. (2018). Furthermore, in the presence of electron-electron interactions these domain walls will form helical Luttinger liquids Wu et al. (2006). Though Sr2Pb3 is insufficiently insulating to experimentally isolate these effects, future, bulk-insulating, nonsymmorphic Dirac insulators could provide an experimental platform for probing and manipulating the physics of time-reversal-invariant topological Luttinger liquids.

Acknowledgements.

*†*Corresponding author. Email: bwieder@princeton.edu (BJW), kane@physics.upenn.edu (CLK), bernevig@princeton.edu (BAB). We thank Aris Alexandradinata for a discussion about the invariant in Eq. (3), and Eugene Mele for fruitful discussions. BJW and CLK acknowledge support through a Simons Investigator grant from the Simons Foundation to Charles L. Kane and through Nordita under ERC DM 321031. ZW and BAB acknowledge support from the Department of Energy Grant No. DE-SC0016239, the National Science Foundation EAGER Grant No. DMR-1643312, Simons Investigator Grants No. ONR-N00014-14-1-0330, No. ARO MURI W911NF-12-1-0461, and No. NSF-MRSEC DMR- 1420541, the Packard Foundation, the Schmidt Fund for Innovative Research, and the National Natural Science Foundation of China (No. 11504117). YK and AMR thank the National Science Foundation MRSEC Program for support under DMR-1120901 and acknowledge the HPCMO of the U.S. DOD and the NERSC of the U.S. DOE for computational support. BAB wishes to thank École Normale Supérieure, UPMC Paris, and the Donostia International Physics Center for their generous sabbatical hosting during some of the stages of this work.

Appendix A Fermion Doubling in Quasi-2D Crystals

In this section, we discuss how fermion doubling theorems in 2D, which apply to both truly 2D wallpaper-group and quasi-2D layer-group systems, constrain band structures, and how apparent exceptions to them manifest on the surfaces of bulk 3D topological phases. Specifically, we develop fermion doubling constraints for the larger set of 2D Hamiltonians invariant under layer groups, which consequently also apply to the subset of those Hamiltonians also invariant under wallpaper groups. In this work, we use the phrase “Dirac fermion” only to refer to fermions with point-like fourfold degeneracies and linear dispersion in the two in-plane reciprocal lattice directions of a two-dimensional surface or a quasi-two-dimensional material, or in all three independent directions if it is in a three-dimensional bulk. In this nomenclature, the surface states of a topological insulator are therefore not “Dirac points,” but are twofold-degenerate linearly dispersing fermions by another name (in a topological insulator they are linearly dispersing Kramers pairs).

This is an unfortunate consequence of the competing and contradicting contexts in which the name “Dirac point” has been previously applied. In earlier works, it has been used to describe the twofold degeneracies in 2D in spinless graphene Semenoff (1984); DiVincenzo and Mele (1984), the spinful twofold degeneracies on the 2D surface of a topological insulator Hsieh et al. (2008), and the spinful fourfold degeneracies in the 3D bulks of Na3Bi Liu et al. (2014b) and Cd3As2 Borisenko et al. (2014). Although it differs slightly from the usage in high-energy physics for two dimensions Weinberg et al. (1995), we choose to designate fourfold linear degeneracies in two and three dimensions as condensed matter realizations of “Dirac fermions.” In so doing, we keep the classification of the bulk degeneracies in graphene as Dirac points, but are forced to expand the $2\times 2$ matrices from Semenoff, 1984; DiVincenzo and Mele, 1984 to explicitly include an additional spin degree of freedom. As we show in A.2, this definition has precise symmetry and topology underpinnings; it enforces the relationship of these fourfold degeneracies with crystalline-symmetry-enhanced fermion doubling, and it preserves their role as the symmetry-pinned quantum critical points between trivial and topological insulating phases.

A.1 2D Fermion Doubling for Twofold-Degenerate Linear Fermions

In 2D, a system may host linearly-dispersing twofold-degenerate fermions. Such gapless fermions may exist as fine-tuned points or, if additional symmetries are present to protect them, may exist in pairs in a stable phase (i. e.  spinless graphene). Here, we are interested in stable phases. Thus, we restrict ourselves to only discussing systems for which all possible symmetry-allowed hopping terms have been included. For example, the critical point separating two topologically distinct insulating phases in a two-band model features a single twofold-degenerate gapless fermion. However, without imposing additional symmetries, this fermion can be gapped, and, generically, the system will be insulating at half-filling.

The simplest example of a symmetry protecting a twofold linear fermion occurs in a crystal with time-reversal symmetry, $\mathcal{T}$, satisfying $\mathcal{T}^{2}=-1$. This requires states to be twofold-degenerate at the Time-Reversal-Invariant-Momenta (TRIMs) by Kramers’ theorem. Twofold-degenerate fermions can also appear pairwise along high-symmetry lines and planes, with symmetry stabilization coming from a combination of crystalline and time-reversal symmetries.

In 2D systems, these symmetry-protected gapless points come in pairs, a consequence of the parity anomaly Alvarez-Gaumé and Witten (1984); Redlich (1984); Jackiw (1984). We illustrate this result in the case of a single twofold-degenerate linear fermion at a TRIM protected by time-reversal symmetry, $\mathcal{T}=i\sigma^{y}K$, described by the following 2D $k\cdot p$ theory:

[TABLE]

There is a single remaining Pauli matrix, $\sigma^{z}$, which anticommutes with all of the terms in Eq. (5) and opens a gap. The mass term $V_{m}=m\sigma^{z}$, which could originate from an external magnetic field or mean-field magnetic order, breaks $\mathcal{T}$ symmetry and gaps locally to a $k\cdot p$ theory of a Chern insulator, with two bands of winding number $C$ and $C+1$, respectively, for some $C\in\mathbb{Z}$. If Eq. (5) is a complete description of the low-energy physics, a contradiction arises: because $\{V_{m},\mathcal{T}\}=0$, the two gapped phases that result from choosing opposite signs of $m$ are related to each other by transformation under $\mathcal{T}$. In particular, the band with Chern number $C$ when $m>0$ is related to the band with Chern number $C+1$ when $m<0$. Since $C$ is odd under time-reversal, $-C=C+1$. This condition cannot be satisfied by $C\in\mathbb{Z}$.

Thus, the hypothesis that the stable gapless fermion (Eq. (5)) is a complete description of the low-energy physics cannot be true. For a system with a single twofold fermion, even at $V_{m}=0$, time-reversal symmetry must be anomalously broken by terms beyond the $\mathbf{k}\cdot\mathbf{p}$ level, and the system forms an anomalous Hall state with $|C|=1/2$ (here $C$ need not be integer since the system is gapless). In order for time-reversal to remain unbroken, there must be a compensating second degeneracy point somewhere else in the Brillouin zone, also with $|C|=1/2$. The only details of time-reversal symmetry that entered into the preceding argument are that $\mathcal{T}V_{m}\mathcal{T}^{-1}=-V_{m}$, and so the general result remains true for any twofold symmetry that protects a gapless fermion in two dimensions111This is due to the fact that a Pauli-Villars regulator is a function of $V_{m}$, and so this regularization must break the symmetry, and so we derive the same anomaly-generating functional as in the literature.. As the generic Hamiltonian of a twofold-degenerate linear fermion in 2D only contains two linear terms which exhaust two of the Pauli matrices, the statement that the fermion is symmetry-protected also implies that a mass term proportional to the remaining Pauli matrix will be odd under the symmetry it breaks, and that it will anticommute with the Hamiltonian and necessarily open a gap.

We can gain some further intuition about the parity anomaly by noting that the anomalous Hall conductance $C$ is related to the Berry phase at the Fermi surface by Haldane (2004):

[TABLE]

where $P$ is the path-ordering operator. Let us take a compact Brillouin zone with $N$ gapless fermions, and let the Fermi level be above all bands. By evaluating the Berry phase we find that $2\pi iC=\log((-1)^{N})$. However, with time-reversal symmetry, we also have that $C$ is the total Chern number of all bands, and hence $C=0$. Thus we conclude $N\in 2\mathbb{Z}$. In the context of three-dimensional topological insulators, we note that each surface taken in isolation is a 2D system with a single gapless fermion. From the above discussion we thus recover the well-known result that each surface of a topological insulator has a half-quantized anomalous Hall conductivity, and that only when both surfaces are connected by a bulk can the system be described in an anomaly-free way Fu and Kane (2007); Qi et al. (2008); Mulligan and Burnell (2013). This is the sense in which a topological insulator is sometimes said to “cheat” the doubling of twofold-degenerate linear fermions. Of particular note is that the strongly spin-orbit coupled topological crystalline insulating phases discovered to date, the mirror TCI in SnTe Hsieh et al. (2012); Tanaka et al. (2012); Dziawa et al. (2012); Xu et al. (2012) and the hourglass insulator in KHgSb Wang et al. (2016); Alexandradinata et al. (2016); Ma et al. (2017), exhibit time-reversed pairs of twofold-degenerate surface fermions, and thus do not carry the same relationship with fermion doubling as does a QSH insulator. Rather their individual surfaces exhibit integer-quantized Hall conductivities and conversely do not anomalously violate parity.

The same logic can be applied to a 3D material by replacing the mass term $m\sigma^{z}$ by $v_{z}k_{z}\sigma^{z}$; in this case the gapless point would be a Weyl point of Chern number $+1$ and the two gapped Hamiltonians of opposite-signed mass lie in the planes above and below it. This expresses the so-called “descent relation” between the parity anomaly in two dimensions and the chiral anomaly in three dimensions Alvarez-Gaumé et al. (1985). If $\mathbf{k}_{z}$ is periodic, this would imply that, absent another Weyl point, two systems with different Chern number (above and below the Weyl point) could be adiabatically connected through the BZ boundary, which is impossible. To avoid this contradiction, the doubling theorem then requires that the low-energy physics cannot be described by only a single Weyl point: there must be another Weyl point or other band crossing at the Fermi level. In 3D, this is the celebrated Nielsen-Ninomiya theorem Nielsen and Ninomiya (1981).

A.2 Quasi-2D Fermion Doubling for Dirac Fermions

We now extend these arguments to show why fourfold-degenerate quasi-2D Dirac fermions cannot be stabilized as the only nodal features in a metallic phase with time-reversal symmetry at a given energy. As in the twofold-degenerate case, while many models might display Dirac fermions upon fine-tuning, here we are interested in robust Dirac fermion phases. Thus, we only consider systems that display Dirac fermions when all symmetry-allowed hopping terms are present. The crux of the arguments in this section was originally highlighted in Young and Kane, 2015; Fu et al., 2007.

When Dirac points occur off of the TRIMs, which can only occur in quasi-2D systems invariant under nonsymmorphic layer groups Wieder and Kane (2016), time-reversal requires that they come in pairs Yang and Nagaosa (2014). The remaining Dirac points which occur in 2D or quasi-2D systems are filling-enforced and pinned to the TRIMs by crystalline symmetries and time-reversal. As shown in Wieder and Kane, 2016, these filling-enforced, high-symmetry Dirac points can only occur in layer or wallpaper groups where either inversion anticommutes with a twofold nonsymmorphic symmetry or where two perpendicular nonsymmorphic symmetries anticommute, such as the two glides in $pgg$.

Consider first a quasi-2D spinful system with inversion symmetry $\mathcal{I}$ and a screw rotation, $s_{2y}=t_{y/2}C_{2y}$. At $k_{y}=\pi$, $\mathcal{I}^{2}=s_{2y}^{2}=+1$ and $\{\mathcal{I},s_{2y}\}=0$. In a time-reversal-invariant system, these symmetries require a four-dimensional corepresentation: more generally, for any Hamiltonian invariant under two symmetries, $A$ and $B$, in addition to $\mathcal{T}$, such that $\{A,B\}=[A,\mathcal{T}]=[B,\mathcal{T}]=0$ and $A^{2}=1$, any eigenstate of the Hamiltonian, $\psi$, which is also an eigenstate of $A$, is part of a fourfold-degenerate quartet of orthogonal states, $\psi,\mathcal{T}\psi,B\psi,\mathcal{T}B\psi$. In our example, there is a fourfold degeneracy at each TRIM with $k_{y}=\pi$. In general, the combination of $\mathcal{I}$ and a twofold nonsymmorphic symmetry will always mandate that there are no fewer than two Dirac points for a given filling: since $\mathcal{I}$ anticommutes with the nonsymmorphic symmetry at two of the TRIMs in 2D and $\mathcal{I}^{2}=+1$, there are always at least two TRIM points with fourfold degeneracies.

In the case where there are two perpendicular nonsymmorphic symmetries, but no inversion symmetry, the obstruction to forming an isolated, stable Dirac fermion takes a slightly different form. Consider the trivial phase of wallpaper group $pgg$, for example, which is characterized by $\mathcal{T}$ and two glides, $g_{x,y}=\{m_{x,y}|\frac{1}{2}\frac{1}{2}\}$. At $k_{x}=k_{y}=\pi$, $\{g_{x},g_{y}\}=0$ and $g_{x}^{2}=g_{y}^{2}=+1$, and therefore a Dirac point exists. However, this condition is only met at this corner TRIM, and no other fourfold degeneracies are allowed elsewhere in this system. In this case, the Dirac point is obstructed from being alone by the filling-enforced hourglass structures also required to exist by the presence of singly degenerate eigenstates of $g_{x,y}$. In these systems, the Dirac point occurs at the same filling as four 2D twofold-degenerate linear fermions, and is also prevented from being alone and stable at any filling.

We pause to briefly note here that the expression “nonsymmorphic symmetry,” is a slight abuse of terminology, though one rampant in this field. For precision, and for consistency with the terminology employed in Wang et al., 2016; Alexandradinata et al., 2016, we define for this work a nonsymmorphic symmetry to be a generating element $g$ of the maximal fixed-point-free subgroup(s), or Bieberbach subgroup(s) Watanabe et al. (2015); Wieder and Kane (2016), of a nonsymmorphic wallpaper or space group $\mathcal{G}$, modulo full lattice translations $T$. Consider, for example, wallpaper group $pgg$, generated by $g_{x}=t_{y/2}M_{x}$ and $g_{y}=t_{x/2}M_{y}$. There are two maximal fixed-point-free subgroups of $pgg$: one generated by $g_{x}$ and $t_{x}$ and one generated by $g_{y}$ and $t_{y}$, where $t_{x,y}$ is a full lattice translation. Both of these groups are isomorphic to the Bieberbach wallpaper group Wieder and Kane (2016) $pg$, and when the generating elements are taken modulo the group of full lattice translations respectively give $g_{x}$ and $g_{y}$. For a more complicated example, consider space group 14 $P2_{1}/b$, generated by $s_{2x}=t_{x/2}t_{z/2}C_{2x}$, $t_{y}$, $t_{z}$, and spatial inversion $\mathcal{I}$ about the origin Bradley and Cracknell (1972). This group also has two maximal fixed-point-free subgroups, which when taken modulo the group of full lattice translations respectively give $s_{2x}$ (isomorphic to Bieberbach space group 4 $P2_{1}$ modulo $T$) and $g_{x}=\mathcal{I}\times s_{2x}$ (isomorphic to Bieberbach space group 7 $Pb$ modulo $T$).

We pose an explanation for this obstruction, which is similar to the resolution of the parity anomaly in the previous section. Suppose a combination of symmetries were to allow a single stable Dirac point. The $k\cdot p$ model Hamiltonian at a corner TRIM in our previous system with $\mathcal{I}$ and $s_{2y}$ is described by a linear $k\cdot p$ model:

[TABLE]

where $\mathcal{I}=\tau^{z}$, $s_{2y}=\tau^{y}\sigma^{y}$, and $\mathcal{T}=i\sigma^{y}K$. The four Dirac matrices present in $\mathcal{H}$ span the space of symmetry-allowed matrices. Thus, the system supports a robust fourfold-degenerate gapless fermion stabilized by crystal-symmetries.

As before, we can examine the consequences of locally breaking one of the symmetries. To guarantee that the $k\cdot p$ Hamiltonian is gapped everywhere, we seek a mass term to anticommute with all of the terms in $\mathcal{H}$. Generically, the Clifford algebra of Dirac matrices is spanned by four $\mathcal{T}$-odd matrices which couple to crystal momenta, and one $\mathcal{T}$-even matrix, here $V_{m}=m\tau^{z}$. For either sign of $m$, the resulting phase is 2D, gapped, and $\mathcal{T}$-symmetric with $\mathcal{T}^{2}=-1$; a quantum spin Hall (QSH) index can thus be defined. Noting that at $\mathbf{k}=0$, $\mathcal{H}=V_{m}=m\mathcal{I}$, a generic feature of Hamiltonians restricted to TRIM points, c. f. Bernevig and Hughes, 2013, we see by the Fu-Kane formula Fu and Kane (2007) for the QSH index that occupied bands for $m>0$ and $m<0$ have opposite parity $\mathbb{Z}_{2}$ indices. This is expected, as it is known that such a Dirac point is the boundary between a trivial and a topological insulator Bernevig and Hughes (2013).

To show the need for fermion doubling, we expand to the full BZ of a hypothetical system where this fermion is the only feature at the Fermi energy and show that there is a contradiction. We label the Bloch wavefunctions of the phase when $m>0$ by $|u(k)\rangle$ and those of the phase with $m<0$ by $|\tilde{u}(\tilde{k})\rangle=s_{2y}|u(k)\rangle$, where $\tilde{k}=s_{2y}k$. The integral of the pfaffian of the matrix $w_{ij}=\langle u(-k)_{i}|\mathcal{T}|u(k)_{j}\rangle$ gives the QSH $\mathbb{Z}_{2}$ topological invariant Kane and Mele (2005). Consider relating this matrix for one gapped phase to the other by the operation of the broken symmetry, $s_{2y}$: $\tilde{w}_{ij}=\langle u(-k)_{i}|s_{2y}^{\dagger}\mathcal{T}s_{2y}|u(k)_{j}\rangle=\langle u(-k)_{i}|\mathcal{T}s_{2y}^{\dagger}s_{2y}|u(k)_{j}\rangle=w_{ij}$ because $[\mathcal{T},s_{2y}]=0$. Therefore, having the same $w$ matrix, the two phases have the same QSH invariant, contradicting the earlier Fu-Kane requirement that the two insulating phases are topologically distinct.

The resolution of this is a fermion doubling requirement for 2D fourfold-degenerate Dirac points. Specifically, a closed 2D crystal, invariant under either just a wallpaper group or both a wallpaper group and its layer supergroup, cannot host a single symmetry-stabilized Dirac point at the Fermi level. In the cases where the Fermi surface is gapped except for exactly two Dirac points, each Dirac point can have a single $\mathcal{T}$-symmetric mass term $m_{1,2}$ such that the overall $\mathbb{Z}_{2}$ QSH invariant, $n$, satisfies $n=\frac{1}{2}[\operatorname{sgn}(m_{1})-\operatorname{sgn}(m_{2})]\mod 2$. Under the action of the broken symmetry operation, the signs of both $m_{1,2}$ are flipped and $n$ is preserved Young and Kane (2015).

Like the topological insulator before it, the nonsymmorphic Dirac insulator that we present in this manuscript “cheats” this fermion doubling by placing each of its two Dirac points on opposite surfaces of a 3D bulk. While they each live alone on a surface, and therefore pin that surface to a 2D QSH transition, the Dirac points (or a Dirac point and four twofold degenerate points) in the combined system of two opposing surfaces respect the fermion doubling requirement.

Appendix B Symmetries and Degeneracies of Wallpaper Groups $pgg$ and $p4g$

It has been extensively demonstrated that the bulk electronic degeneracies of a 3D crystal are constrained by the allowed irreducible (co)representations of the space group of that crystal, up to topological features, such as Weyl points Young et al. (2012); Steinberg et al. (2014); Wieder et al. (2016); Bradlyn et al. (2016, 2017); Bradley and Cracknell (1972); Mañes (2012); Chang et al. (2017b). Additional degeneracies may be realized in toy models, but those fermions are unlikely to manifest in real materials, as they require symmetry-allowed terms to be manifestly zero. For example, the threefold symmorphic fermions predicted in the toy model in Fulga and Stern, 2017 appear to be realizable using only wallpaper group symmetries, but in fact they additionally require an unphysical anticommuting mirror to artificially constrain their $k\cdot p$ Hamiltonians. In fact, similar symmorphic realizations of Spin-1 Weyl fermions Bradlyn et al. (2016) can be realized in cubic crystals with point group $23$ ($T$) and weak SOC, such as, for example at the $\Gamma$ point in SG 195 Mañes (2012); Chang et al. (2017b). But these fermions are inherently three-dimensional, and are completely characterized by the single-valued corepresentations of these cubic space groups. Therefore, in this work, to find novel surface degeneracies, we explore the allowed band degeneracies of the time-reversal-symmetric wallpaper groups as captured by their double-valued corepresentations.

The wallpaper, or plane, groups describe the 17 possible configurations of symmetries on the two-dimensional surface of a three-dimensional, time-reversal-symmetric crystal Conway et al. (2008). To be mathematically precise, in this work, our use of the term “wallpaper groups” is shorthand for the Shubnikov Bradley and Cracknell (1972) wallpaper groups $W_{II}$ formed from index-2 supergroups of their type-I counterparts $G_{I}$ by the addition of time-reversal-symmetry $W_{II}=G_{I}+\mathcal{T}G_{I}$. The irreducible corepresentations of $W_{II}$ are formed by considering the action of $\mathcal{T}$ on the irreducible representations Bradlyn et al. (2017); Bradley and Cracknell (1972) of $G_{I}$. Of these groups, only two have multiple nonsymmorphic symmetries Bradley and Cracknell (1972) (as defined previously in A.2). Specifically, wallpaper groups $pgg$ and $p4g$ contain glide lines in the $x$ and $y$ (surface in-plane) directions; $p4g$ has an additional $C_{4z}$ symmetry that is not present in $pgg$. The group $pgg$ is generated entirely by the two glides:

[TABLE]

while wallpaper group $p4g$ is generated by

[TABLE]

The product of the two generators for $p4g$ yields two additional symmorphic mirror symmetries, $\{m_{x+y}|\frac{1}{2}\frac{1}{2}\}$ and $\{m_{x-y}|\frac{1}{2}\bar{\frac{1}{2}}\}$. Though containing translations, these symmetries are not glides, as glides are, by definition, free of fixed points in position space, and $\{m_{x\pm y}|\frac{1}{2},\pm\frac{1}{2}\}$ leaves the line $y=\pm(-x+\frac{1}{2})$ invariant. In Fig. 1, we show the locations of the glides, mirrors, and $C_{4z}$ rotation centers for both wallpaper groups realized with a two-site unit cell.

In crystal momentum space, the symmetry generators enforce required band groupings along lines and at points. In Figure 7, we show one quarter of the surface BZ and identify relevant lines and points with the letters (a), (b), and (c). Lines sharing the same letter obey the same symmetry restrictions, though for wallpaper group $pgg$, which lacks a $C_{4z}$ symmetry, they may individually display different features.

The lines designated (a) in Fig. 7 are $\bar{\Gamma}\bar{X}$ and $\bar{\Gamma}\bar{Y}$ and host singly degenerate bands which are eigenstates of $g_{y}$ and $g_{x}$, respectively. As detailed extensively in Young and Kane, 2015; Wieder and Kane, 2016; Wang et al., 2016; Alexandradinata et al., 2016, singly degenerate bands along lines invariant under a glide form hourglass or quantum spin Hall flows.

Bands along lines designated (b), $\bar{X}\bar{M}$ and $\bar{Y}\bar{M}$, are also eigenstates of $g_{x}$ and $g_{y}$, respectively. However, unlike the bands along (a), they are doubly degenerate, because the symmetry operation $\mathcal{T}g_{y}$ ($\mathcal{T}g_{x}$) leaves the line $\bar{X}\bar{M}\ (\bar{Y}\bar{M})$ invariant and squares to $-1$.

We now show that all potential band crossings along the line $\bar{X}\bar{M}$ are avoided by utilizing the fact that bands along this line are also eigenstates of $g_{x}$. If $g_{x}|(k_{x}=0,\pi;k_{y}),+\rangle=g_{x}|+\rangle=ie^{ik_{y}/2}|+\rangle$, then, using the commutation relation $g_{x}g_{y}=-g_{y}g_{x}t_{x}t_{-y}$, and by evaluating the resulting expression along $\bar{X}\bar{M}$ at $k_{x}=\pi$,

[TABLE]

Thus, the Kramers partners $|+\rangle$ and $\mathcal{T}g_{y}|+\rangle$ have opposite $g_{x}$ eigenvalues. We note that eigenstates of $g_{y}$ along $\bar{Y}\bar{M}$ must behave in the same manner by $x\leftrightarrow y$ exchange symmetry. Since all sets of Kramers pairs along line (b) have the same pair of eigenvalues, they belong to the same corepresentation, and so crossings between these bands will generically be avoided.

Finally, at the point $\bar{M}$, labeled (c), bands are fourfold-degenerate, a unique property of wallpaper groups $pgg$ and $p4g$. This can be easily seen because the $\bar{M}$ point is invariant under $g_{x},g_{y}$, and $\mathcal{T}$. Thus, if $\psi$ is an eigenstate of $g_{x}$, $\psi,\mathcal{T}\psi,g_{y}\psi$, and $\mathcal{T}g_{y}\psi$ form a degenerate quartet of states; clearly the Kramers pairs are orthogonal and since the first two states have the opposite $g_{x}$ eigenvalue as the last two states, they are also orthogonal.

This algebra can only be satisfied in the strong spin-orbit coupled wallpaper groups by two perpendicular glides. Furthermore, an examination of two-dimensional filling constraints in Wieder and Kane, 2016 confirms that no other algebra in 2D can enforce fourfold point-like band degeneracies. Therefore a fourfold point degeneracy can only be hosted on the surfaces of time-reversal-symmetric, strong spin-orbit, three-dimensional crystals invariant under the symmetries of $pgg$ or $p4g$.

We now determine the low-energy band dispersion near a fourfold crossing at (c). Returning to the two-site unit cell depicted in Fig. 1, we construct a $k\cdot p$ theory around $\bar{M}$. Letting $\tau$ represent the sublattice degree of freedom and $\sigma$ the local spin degree of freedom, we choose the representation $g_{x/y}=\tau^{y}\sigma^{x/y}$ and $\mathcal{T}=i\sigma^{y}\mathcal{K}$. Then, to linear order in $k_{x,y}$, the most general allowed Hamiltonian is given by,

[TABLE]

This is the equation of a linearly-dispersing fourfold-degenerate 2D fermion. Therefore bands at $\bar{M}$ generically touch in Dirac points. It is worth noting that unlike the 3D Dirac fermions characterized in Wang et al., 2013, 2012, which are optional degeneracies created by band inversion, all band multiplets at $\bar{M}$ in $pgg$ and $p4g$ must be at least fourfold-degenerate, and therefore the Dirac fermions in these wallpaper groups are instead more closely related to the filling-enforced 2D Dirac points proposed in Young and Kane, 2015; Wieder and Kane, 2016; Young and Wieder, 2017. Like the 2D inversion-broken magnetic Dirac points in Young and Wieder, 2017, these surface Dirac points also have generically non-degenerate cones; the cones here are only degenerate along the glide lines $\bar{X}\bar{M}$ and $\bar{Y}\bar{M}$. In our $k\cdot p$ model, this behavior emerges upon the introduction of the quadratic term $v_{xy}k_{x}k_{y}\tau^{z}$.

Inducing a $\mathcal{T}$-symmetric mass term $V_{m}=m\tau^{z}$ breaks each surface into a 2D trivial or topological insulator. On a single surface, as discussed in more detail in A.2, domain walls between regions with different signs of $m$ will host 1D QSH edge modes, which are topologically protected by time-reversal symmetry alone. These domain walls have been shown to host Luttinger liquid physics Wu et al. (2006), and intersections between them can form effective quantum point contacts with switching behavior characterized by universal scaling functions and critical exponents Teo and Kane (2009). Additionally, these surface domain walls would provide an improvement over recent efforts to test this physics in qualitatively similar domain walls in gapped bilayer graphene Castro Neto et al. (2009); McCann and Fal’ko (2006); Zhang et al. (2011); Ohta et al. (2006); Oostinga et al. (2007), which have also been proposed as Luttinger liquids Killi et al. (2010); Wieder et al. (2015). Notably, as the domain walls in bilayer graphene are protected by valley index, they are quite sensitive to disorder, and thus far Luttinger liquid physics in bilayer graphene has not been observed Ju et al. (2015); Li et al. (2016). Domain wall modes on gapped nonsymmorphic Dirac insulator surfaces, conversely, should remain robust against nonmagnetic disorder.

Appendix C Tight-Binding Notation

Here we provide the tight-binding notation that will be used in subsequent sections, following Alexandradinata et al., 2016. In the unit cell labeled by the Bravais lattice vector, $\mathbf{R}$, the wavefunction corresponding to an orbital labeled by $\alpha$ at position $\mathbf{R}+\mathbf{r}_{\alpha}$ is denoted by $|\phi_{\mathbf{R},\alpha}\rangle$. The Fourier transformed operators are given by,

[TABLE]

The single-particle Hamiltonian, $\hat{H}$, defines a tight-binding Hamiltonian,

[TABLE]

whose eigenstates are denoted $|u^{n}(\mathbf{k})\rangle$.

The Fourier transform in Eq. (12) shows that the Hamiltonian is not necessarily invariant under a shift of a reciprocal lattice vector, $\mathbf{G}$. Instead,

[TABLE]

where $\hat{V}(\mathbf{G})_{\alpha\beta}=\delta_{\alpha\beta}e^{i\mathbf{G}\cdot\mathbf{r}_{\alpha}}$. Thus, we can choose the basis of eigenstates so that,

[TABLE]

C.1 Symmetries

If the lattice is invariant under a spatial symmetry, $\hat{g}$, that acts in real space by taking $\mathbf{r}\rightarrow D_{g}\mathbf{r}+\mathbf{\delta}$, where $D_{g}$ is the matrix that enforces a point group operation and $\mathbf{\delta}$ is a (perhaps fractional) lattice translation, then $\hat{g}$ acts on states by:

[TABLE]

where $\mathbf{R}^{\prime}=D_{g}(\mathbf{R}+\mathbf{r}_{\alpha})+\delta-\mathbf{r}_{\beta}$ is a Bravais lattice vector and $U_{g}$ is a unitary matrix. It is convenient to define the Fourier transformed operator, $\hat{g}_{\mathbf{k}}$, which can have explicit $\mathbf{k}$ dependence when it acts on Bloch states, by:

[TABLE]

Thus, $\hat{g}_{\mathbf{k}}$ separates into a product of a $\mathbf{k}$-dependent phase and a matrix, $U_{g}$, that rotates the orbitals:

[TABLE]

If $\hat{H}$ is invariant under $\hat{g}$, then Eq. (17) shows that:

[TABLE]

We can follow the same procedure for an antiunitary operator, $\mathcal{T}_{g}\equiv\mathcal{T}\hat{g}$, to find:

[TABLE]

where $\mathcal{K}$ is the complex conjugation operator. Similarly to Eq. (19),

[TABLE]

C.2 Projector onto Occupied States

We define the projector onto the $n_{\rm occ}$ occupied states:

[TABLE]

which satisfies, using Eq. (15),

[TABLE]

Given a spatial symmetry, $\hat{g}$, Eq. (19) shows that $\hat{g}_{\mathbf{k}}|u^{n}(\mathbf{k})\rangle$ has the same energy as $|u^{n}(\mathbf{k})\rangle$ and is a state at momentum $D_{g}\mathbf{k}$. Hence, the projector onto the occupied states at momentum $D_{g}\mathbf{k}$ is given by:

[TABLE]

In subsequent sections, we will want to know how $\hat{g}_{\mathbf{k}+\mathbf{G}}$ is related to $\hat{g}_{k}$. Plugging Eq. (14) into Eq. (19):

[TABLE]

Thus,

[TABLE]

For an antiunitary symmetry, $\mathcal{T}_{g}\equiv\mathcal{T}\hat{g}$, the analogous equations are:

[TABLE]

and

[TABLE]

Appendix D Wilson Loops

A precise way to distinguish the distinct surface connectivities is obtained through the eigenvalues of Wilson loops, which, as we will elaborate, also give information about the edge state spectrum Fu and Kane (2006); Ryu et al. (2010); Soluyanov and Vanderbilt (2011); Yu et al. (2011); Taherinejad et al. (2014); Alexandradinata et al. (2014a); Fidkowski et al. (2011); Alexandradinata et al. (2016). Here, we are interested in the Wilson loop matrix:

[TABLE]

where $P$ indicates that the integral is path-ordered, and $A_{z}(\mathbf{k})_{ij}\equiv i\langle u^{i}(\mathbf{k})|\partial_{k_{z}}u^{j}(\mathbf{k})\rangle$ is a matrix whose rows and columns correspond to filled bands. The eigenvalues of $\mathcal{W}$ are gauge invariant and of the form $e^{i\theta(k_{x},k_{y})}$, i.e., they are independent of the ‘base point,’ $k_{z0}$.

In this appendix, we show that a space group with two glides, $g_{x,y}\equiv\{m_{x,y}|\frac{1}{2}\frac{1}{2}0\}$, as well as time-reversal symmetry, $\mathcal{T}$, has the following constraints on its $z$-directed Wilson loop eigenvalues:

• •

Along the lines $(k_{x},\pi)$ and $(\pi,k_{y})$, the Wilson loop eigenvalues are doubly degenerate, due to the antiunitary symmetries $g_{x}\mathcal{T}$ and $g_{y}\mathcal{T}$, respectively (see D.2).

• •

At the $(\pi,\pi)$ point, the doubly degenerate bands meet at a fourfold band crossing (see D.3).

• •

Along the $\bar{X}\bar{M}\ (\bar{Y}\bar{M})$ and $\bar{\Gamma}\bar{Y}\ (\bar{\Gamma}\bar{X})$ lines, bands can be labeled by the eigenvalues of $g_{x}\ (g_{y})$ (see D.3).

The gauge invariance of the loops allows us to write down a topological invariant that distinguishes the states, as derived in E.1.

D.1 Discretized Wilson Loop

For completeness, following Alexandradinata et al., 2016, we derive a discretized version of the Wilson loop (Eq. (29)), which is useful for clarity when deriving symmetry constraints. Using the projector onto occupied states, $\hat{\mathcal{P}}$, defined in Eq. (22), we derive the discretized Wilson loop matrix,

[TABLE]

where $k_{\perp}\equiv(k_{x},k_{y})$ and in the last line we have defined the ordered product of projectors,

[TABLE]

Eq. (30) shows that the discretized Wilson loop can be written in the basis-invariant form,

[TABLE]

By applying Eq. (23) to Eq. (32), for a reciprocal lattice vector $\mathbf{G}$,

[TABLE]

which shows that the Wilson loop eigenvalues are invariant under shifts of $\mathbf{G}$, whether $\mathbf{G}$ is along or perpendicular to the direction of the loop.

D.2 Effect of Time-Reversal-Like Symmetries on the Wilson Loop

Here, we consider an antiunitary symmetry, $\mathcal{T}_{g}\equiv\mathcal{T}\hat{g}$, where $\hat{g}$ rotates $k_{\perp}$ but does not affect $k_{z}$ and $\hat{g}$ does not include a real space translation along $\hat{z}$. If we write the action of $\hat{g}$ in real space by $\mathbf{r}\rightarrow D_{g}\mathbf{r}+\mathbf{\delta}$, these conditions require that $D_{g}(k_{\perp},k_{z})=(D_{g}k_{\perp},k_{z})$ and $\delta=(\delta_{x},\delta_{y},0)$. Using Eq. (20), these properties imply $\mathcal{T}_{g,(k_{\perp},k_{z})}=\mathcal{T}_{g,(k_{\perp},0)}\equiv\mathcal{T}_{g,k_{\perp}}$. We would like to relate $\mathcal{W}_{(-D_{g}k_{\perp},0)}$ to $\mathcal{W}_{(k_{\perp},0)}$, to which end we compute:

[TABLE]

We have shown that, $\mathcal{T}_{g,k_{\perp}}^{-1}\mathcal{W}_{(-D_{g}k_{\perp},0)}\mathcal{T}_{g,k_{\perp}}=\mathcal{W}_{(k_{\perp},0)}^{\dagger}$. Thus, if $\psi(k_{\perp})$ is an eigenvector of $\mathcal{W}_{(k_{\perp},0)}$ with eigenvalue $e^{i\theta(k_{\perp})}$, then $\mathcal{T}_{g,k_{\perp}}\psi(k_{\perp})$ is an eigenvector of $\mathcal{W}_{(-D_{g}k_{\perp},0)}$ with the same eigenvalue. Furthermore, if for some $k_{\perp}$, $\mathcal{T}_{g}$ leaves $k_{\perp}$ invariant up to a reciprocal lattice vector (i.e, $-D_{g}(k_{\perp},k_{z})=(k_{\perp}+G_{\perp},-k_{z})$, where $(G_{\perp},0)$ is a reciprocal lattice vector), then from Eq. (33), $\hat{V}(G_{\perp})\mathcal{T}_{g,k_{\perp}}\psi(k_{\perp})$ is an eigenstate of $\mathcal{W}_{(k_{\perp},0)}$, also with the same eigenvalue as $\psi(k_{\perp})$. Since $\hat{V}(G_{\perp})\mathcal{T}_{g,k_{\perp}}$ is an antiunitary symmetry that squares to $-1$, $\psi(k_{\perp})$ and $\hat{V}(G_{\perp})\mathcal{T}_{g,k_{\perp}}\psi(k_{\perp})$ are orthogonal. Thus, at momenta whose projections onto $(k_{x},k_{y})$ are invariant under $\mathcal{T}_{g}$ up to a reciprocal lattice vector, eigenstates of $\mathcal{W}_{(k_{x},k_{y},0)}$ come in Kramers pairs.

D.3 Effect of Unitary Symmetries that Leave $k_{z}$ Invariant

Here we consider a unitary symmetry, $\hat{g}$, which leaves $k_{z}$ invariant and does not translate in the $\hat{z}$ direction, i.e., $D_{g}(k_{\perp},k_{z})=(D_{g}k_{\perp},k_{z})$ and $\delta=(\delta_{x},\delta_{y},0)$. Using Eq. (17), $g_{(k_{\perp},k_{z})}=g_{(k_{\perp},0)}\equiv g_{k_{\perp}}$. Eq. (24) shows that $\hat{V}(2\pi\hat{z})\hat{\Pi}(D_{g}k_{\perp},k_{z})=\hat{V}(2\pi\hat{z})\hat{g}_{k_{\perp}}\hat{\Pi}(k_{\perp},k_{z})\hat{g}_{k_{\perp}}^{\dagger}=\hat{g}_{k_{\perp}}\hat{V}(2\pi\hat{z})\hat{\Pi}(k_{\perp},k_{z})\hat{g}_{k_{\perp}}^{\dagger}$. Thus, by definition,

[TABLE]

Now specialize to momenta invariant under $\hat{g}$ up to a reciprocal lattice vector, so that $D_{g}(k_{\perp},k_{z})=(D_{g}k_{\perp},k_{z})=(k_{\perp}+G_{\perp},k_{z})$, where $(G_{\perp},0)$ is a reciprocal lattice vector. At these momenta, $\mathcal{W}_{(D_{g}k_{\perp},k_{z})}=\mathcal{W}_{(k_{\perp}+G_{\perp},k_{z})}$, which, combined with Eqs. (35) and (33), shows that at these momenta, $\left[\mathcal{W}_{(k_{\perp},k_{z})},\hat{V}(G_{\perp})\hat{g}_{k_{\perp}}\right]=0$. Thus, at values of $k_{\perp}$ that are invariant under $D_{g}$ up to a reciprocal lattice vector, the Wilson loop, $\mathcal{W}_{(k_{\perp},k_{z})}$, and the operator $\hat{V}(G_{\perp})\hat{g}_{k_{\perp}}$ can be simultaneously diagonalized. Wilson loop bands can then be labeled by their $\hat{V}(G_{\perp})\hat{g}$ eigenvalue in exactly the same way as energy bands can be labeled by their $\hat{g}_{k_{\perp}}$ eigenvalue.

We now apply the results of this section and the previous section to the glide symmetries defined in the main text, $g_{x,y}\equiv\{M_{x,y}|\frac{1}{2}\frac{1}{2}0\}$. At $k_{\perp}=(\pi,\pi)$, $D_{g_{x}}(k_{\perp},k_{z})=(k_{\perp},k_{z})-2\pi\hat{x}$ and $D_{g_{y}}(k_{\perp},k_{z})=(k_{\perp},k_{z})-2\pi\hat{y}$. From the previous paragraph, the Wilson loop operator $\mathcal{W}_{(k_{\perp},k_{z})}$ commutes with both $\hat{V}(-2\pi\hat{x})g_{x,k_{\perp}}$ and $\hat{V}(-2\pi\hat{y})g_{y,k_{\perp}}$, while $\{\hat{V}(-2\pi\hat{x})g_{x,k_{\perp}},\hat{V}(-2\pi\hat{y})g_{y,k_{\perp}}\}=0$ (the last statement follows from Eq. (26) and the fact that $\{g_{x,k_{\perp}},g_{y,k_{\perp}}\}=0$ at $k_{\perp}=(\pi,\pi)$). Furthermore, from D.2, eigenstates of $\mathcal{W}_{(k_{\perp},0)}$ come in Kramers pairs due to time-reversal symmetry, which takes $(k_{\perp},0)$ to $(k_{\perp},0)-2\pi(\hat{x}+\hat{y})$. Thus, if $\psi(k_{\perp})$ is a simultaneous eigenvector of $\mathcal{W}_{(k_{\perp},0)}$ and $\hat{V}(-2\pi\hat{x})g_{x,k_{\perp}}$, it forms a quartet with three other states, $\hat{V}(-2\pi(\hat{x}+\hat{y}))\mathcal{T}_{k_{\perp}}\psi(k_{\perp})$, $\hat{V}(-2\pi\hat{y})g_{y,k_{\perp}}\psi(k_{\perp})$, and $\hat{V}(-2\pi(\hat{x}+\hat{y}))\mathcal{T}_{k_{\perp}}\hat{V}(-2\pi\hat{y})g_{y,k_{\perp}}\psi(k_{\perp})$, which share the same eigenvalue of $\mathcal{W}_{(k_{\perp},0)}$ but are orthonormal (the orthonormality is verified because the last two states have the opposite $\hat{V}(-2\pi\hat{x})g_{x,k_{\perp}}$ eigenvalues as do the first two). Thus, the Wilson loop eigenvalues are fourfold-degenerate at $k_{\perp}=(\pi,\pi)$.

D.4 Effect of Unitary Symmetries that Flip the Sign of $k_{z}$

Here we consider a unitary symmetry, $\hat{g}$, which flips the sign of $k_{z}$ and does not translate in the $\hat{z}$ direction, i.e. $D_{g}(k_{\perp},k_{z})=(D_{g}k_{\perp},-k_{z})$ and $\delta=(\delta_{x},\delta_{y},0)$. As in the previous section, Eq. (17) guarantees $g_{(k_{\perp},k_{z})}=g_{(k_{\perp},0)}\equiv g_{k_{\perp}}$. Using Eq. (24),

[TABLE]

Thus, the Wilson loop eigenvalues at $k_{\perp}$ and $D_{g}k_{\perp}$ come in complex conjugate pairs. Consequently, at momenta whose surface projection is invariant under this symmetry (or invariant up to a shift of a reciprocal lattice vector $(G_{\perp},0)$, according to Eq. (33)), the spectrum of the phase of the Wilson loop eigenvalues is particle-hole symmetric.

Appendix E Topological Invariant

E.1 Single Glide

In Shiozaki et al., 2016, the authors introduce a $\mathbb{Z}_{4}$ invariant to classify strong topological phases of time-reversal-invariant systems with a single glide symmetry; in Xiong and Alexandradinata, 2018, the authors show how to interpret this invariant in terms of Wilson loops. In this section, we refine this invariant for systems with two glide symmetries.

We first summarize the relevant result from Shiozaki et al., 2016. Consider a time-reversal-invariant system with a single glide symmetry, $g_{y}=\{m_{y}|\frac{1}{2},t,0\}$, where $t$ can be a fractional or integer lattice translation. Since $g_{y}^{2}=\{\mathcal{R}|100\}$, where $\mathcal{R}$ indicates a $2\pi$ rotation, Eq. (17) dictates that $g_{y,(k_{x},k_{y},k_{z})}^{2}=-e^{ik_{x}}$, where the minus sign comes because we are considering spinful systems that acquire a minus sign upon a $2\pi$ rotation. Hence, along the glide-invariant planes, $k_{y}=0$ and $k_{y}=\pi$, each energy or Wilson band can be labeled by its $g_{y,(k_{x},k_{y},k_{z})}$ eigenvalue, $\pm ie^{ik_{x}/2}$, and we say it belongs to the $\pm$ glide sector. The $\mathbb{Z}_{4}$ invariant is defined as,

[TABLE]

where $A_{\pm,i}\equiv i\langle u^{\pm}|\partial_{k_{i}}|u^{\pm}\rangle$, $F_{\pm}\equiv\nabla\times A_{\pm}$, the subscript $\pm$ indicates the glide sector (when there is no subscript, there is no glide symmetry and all occupied bands are summed over), and the superscript $I$ indicates one state in a Kramers pair. It is shown in Shiozaki et al., 2016 that $\chi_{y}$ always takes integer values. Because it is a function of Wilson loops, it is also gauge invariant: the first line compensates for any change of gauge implemented in the last two lines. Thus, $\chi_{y}$ at least partially classifies time-reversal-invariant systems with one glide. Using $K$ theory, it is claimed in Shiozaki et al., 2016 that this is a complete classification of strong topological phases for these space groups. In Alexandradinata et al., 2016, the authors show that the completeness of this classification can also be deduced from the possible windings of the glide-allowed Wilson band connectivities.

Before moving to systems with two glide symmetries, it is useful to translate Eq. (37) into a pictorial computation from a plot of the Wilson loop (defined in Eq. (32)) eigenvalues along the line segment $\bar{M}\bar{Y}\bar{\Gamma}\bar{X}$ (notice this is a bent line Alexandradinata and Bernevig (2016) consisting of three segments, not a complete loop):

Draw a horizontal line across the plot. 2. 2.

Count the number of times a positively sloped band in the $+$ sector crosses the line along the $\bar{M}\bar{Y}$ segment and subtract from it the number of times a negatively sloped band in the $+$ sector crosses the line along the same segment. Multiply the total by 2. 3. 3.

Repeat along the $\bar{\Gamma}\bar{X}$ segment. 4. 4.

Along the $\bar{Y}\bar{\Gamma}$ segment, count the number of times any positively sloped band crosses the line and subtract from it the number of times a negatively sloped band crosses the line (since bands along this line are not eigenstates of $g_{y}$, all bands contribute). 5. 5.

Add the numbers from the previous three steps together; taken mod 4, this sum is $\chi_{y}$.

An example is shown in Fig. 8.

E.2 Two Glides

We now consider systems with two glide symmetries, $g_{x}$ and $g_{y}$, which satisfy $g_{x,(k_{x},k_{y},k_{z})}^{2}=-e^{ik_{y}},\ g_{y,(k_{x},k_{y},k_{z})}^{2}=-e^{ik_{x}}$. Such a system can be described by a pair of invariants $(\chi_{x},\chi_{y})$, where $\chi_{x}$ is defined by exchanging $k_{x}$ and $k_{y}$ in Eq. (37). $\chi_{x}$ can be easily computed from a plot of the Wilson loop by following steps 1-4 in the previous section after interchanging $\bar{X}$ and $\bar{Y}$. An example is shown in Fig. 9.

However, not all pairs $(\chi_{x},\chi_{y})$ are compatible with our assumption of an insulating bulk: we now show that a bulk band insulator only permits $\chi_{x}+\chi_{y}=0\mod 2$. The total number of Wilson bands that cross the reference line around the closed loop $\bar{\Gamma}\bar{X}\bar{M}\bar{Y}\bar{\Gamma}$ must be even because the system is gapped and, consequently, the integral of the Berry curvature over any closed surface must be zero (else there would be a Weyl point contained in the bulk region enclosed by the planes $(0\leq k_{x}\leq\pi,0,k_{z}),(\pi,0\leq k_{y}\leq\pi,k_{z}),(0\leq k_{x}\leq\pi,\pi,k_{z}),(0,0\leq k_{y}\leq\pi,k_{z})$, as shown in Fig. 10). Since the bands along the segment $\bar{Y}\bar{M}\bar{X}$ are doubly degenerate (shown in D.2), the parity of the number of bands that cross the reference line around the closed loop $\bar{\Gamma}\bar{X}\bar{M}\bar{Y}\bar{\Gamma}$ is equal to the parity of the number of bands that cross the reference line along the segment $\bar{X}\bar{\Gamma}\bar{Y}$. Since, in the computation of $\chi_{x,y}$, the results of steps 2 and 3 are necessarily even, the parity of the number of bands that cross along $\bar{Y}\bar{\Gamma}$ is exactly $\chi_{y}\mod 2$; similarly, the parity of the number of bands that cross along $\bar{X}\bar{\Gamma}$ is exactly $\chi_{x}\mod 2$. Thus if we disallow Weyl points by mandating a gapped bulk,

[TABLE]

Consequently, there are eight topologically distinct surface phases that describe a gapped system with two glide symmetries. The eight possible Wilson loops corresponding to the pair of invariants are shown in Fig. 2.

Because $\chi_{x}+\chi_{y}=0\mod 2$, we can rewrite each pair of invariants as $(\chi_{x},\chi_{y})=(\chi_{x},\chi_{x}+1-(-1)^{(\chi_{x}-\chi_{y})/2})$, which shows that the eight topological phases are classified by a $\mathbb{Z}_{4}\times\mathbb{Z}_{2}$ index: the $\mathbb{Z}_{4}$ index is given by the $\mathbb{Z}_{4}$ index of a single $g_{x}$ glide, $\chi_{x}$, and the $\mathbb{Z}_{2}$ index is $\eta_{\chi_{x},\chi_{y}}\equiv\frac{1}{2}\left(1-(-1)^{(\chi_{x}-\chi_{y})/2}\right)$. It is straightforward to check that the $\mathbb{Z}_{2}$ index satisfies the desired group addition:

[TABLE]

We now consider what would happen if instead of computing $\chi_{x,y}$ in the first Brillouin zone, we looked at an adjacent Brillouin zone, shifted by $2\pi$ along the $k_{x}$ axis. In this case, the $g_{y}$ eigenvalues change sign, $\pm ie^{ik_{x}/2}\rightarrow\mp ie^{i(k_{x}+2\pi)/2}$, while the $g_{x}$ eigenvalues remain invariant. Consequently, $\chi_{y}\rightarrow-\chi_{y}$, while $\chi_{x}$ remains unchanged. Similarly, if we moved to a Brillouin zone shifted by $2\pi$ in the $k_{y}$ direction, $\chi_{y}$ would be unchanged, but $\chi_{x}\rightarrow-\chi_{x}$. Thus, one might worry that the $\mathbb{Z}_{4}$ invariant is not a robust characterization of the phase. However, the characterization remains robust: while the labels depend on our choice of Brillouin zone, once that choice is made, there are always eight distinct topological phases. Within a given choice of BZ labeling, transitioning between $\chi_{x,y}=1$ and $\chi_{x,y}=3$ requires closing a bulk gap. Because the phases with odd $\chi_{x,y}$ have topological indices that depend on choice of BZ, and Bloch’s theorem requires that all choices of BZ be equivalent under reciprocal lattice translation, physically distinguishing between $\chi_{x,y}=1,3$ for a single, isolated sample requires a glide-polarized measurement scheme. However, one can compare the relative values of $\chi_{x,y}$ for two samples stacked in the $z$-direction by measuring the presence or absence of a boundary mode, which depends on the difference between $\chi_{x,y}$ in each sample. For example, placing two samples with $\chi_{x}=\chi_{y}=1$ next to each other would yield a trivial boundary, whereas placing such a sample next to one with $\chi_{x}=\chi_{y}=3$ would yield a $(2,2)$ topological Dirac point on the interface.

This effect is similar to what occurs in time-dependent adiabatic pumping cycles of topological superconducting Josephson junctions, for which there are two distinct QSH-like phases that in practice are only distinguishable when coupled to other similar systems Zhang and Kane (2014).

E.3 $\mathbb{Z}_{2}$ Topological Invariant in the Presence of Inversion Symmetry

In the presence of inversion symmetry, $\mathcal{I}$, the $\mathbb{Z}_{2}$ strong topological invariant, $\nu$, is given by,

[TABLE]

where the product is over the eight inversion-symmetric points, $\mathbf{k}_{\text{inv}}$, and $\xi_{\mathbf{k}_{\text{inv}}}$ is the product of the inversion eigenvalues of the occupied bands at $\mathbf{k}_{\text{inv}}$ that are not time-reversal partners (since time-reversal partners have the same inversion eigenvalues, there is no ambiguity in choosing just one of the partners). In this section, we show that in the presence of the two glides, $g_{x},\ g_{y}$, and inversion symmetry, Eq. (40) can be simplified to only involve two points:

[TABLE]

Without loss of generality, in this section we choose the crystal origin so that the inversion operator involves no translation. The two glides can be expressed as $g_{x}=\{m_{x}|t_{x},\frac{1}{2},0\},\ g_{y}=\{m_{y}|\frac{1}{2},t_{y},0\}$, where $t_{x,y}=0$ or $\frac{1}{2}$, depending on the space group. We do not consider translations in the $\hat{z}$ direction since these symmetries would not be preserved by a surface parallel to the $xy$-plane. The operators obey the following commutation relation:

[TABLE]

as do the operators under the transformation $x\leftrightarrow y$; $t_{\mathbf{v}}$ indicates a translation by $\mathbf{v}$.

At the two points $(\pi,\pi,0(\pi))$, which project in the surface Brillouin zone to the $\bar{M}$ point, where a Dirac node is located, filled bands come in groups of four that are eigenstates of $g_{y}$: $\psi$, $g_{x}\psi$ and their time-reversed partners, $\mathcal{T}\psi,\mathcal{T}g_{x}\psi$; notice $\psi$ and $g_{x}\psi$ are linearly independent because they have different $g_{y}$ eigenvalues (assume $g_{y}\psi=\pm\psi$; then because $\{g_{x},g_{y}\}=0$ at the $\bar{M}$ point, $g_{y}(g_{x}\psi)=-g_{x}g_{y}\psi=\mp(g_{x}\psi)$); they are not time-reversed partners because their distinct $g_{y}$ eigenvalues are real. At $(\pi,\pi,0)$, for each eigenstate, $\psi$, with inversion eigenvalue $\lambda=\pm 1$, the state $g_{x}\psi$ has inversion eigenvalue $\lambda e^{(-2t_{x}+1)\pi i}$, following Eq. (42). Thus, $\xi_{(\pi,\pi,0)}=\left(-e^{2\pi it_{x}}\right)^{n_{\rm occ}/4}$. Since the same logic holds at the $(\pi,\pi,\pi)$ point, $\xi_{(\pi,\pi,\pi)}=\left(-e^{2\pi it_{x}}\right)^{n_{\rm occ}/4}$, as well. Since $t_{x}=0$ or $\frac{1}{2}$, $\xi_{(\pi,\pi,0)}\xi_{(\pi,\pi,\pi)}=\left(e^{4\pi it_{x}}\right)^{n_{\rm occ}/4}=1$.

We now assume the presence of $C_{4z}$ to show that the points $(0,\pi,0(\pi))$ and $(\pi,0,0(\pi))$ also contribute a factor of $+1$ to $\nu$. Since $\left[C_{4z},\mathcal{I}\right]=0$, an eigenstate at $(0,\pi,0(\pi))$ has a $C_{4z}$ partner at $(\pi,0,0(\pi))$ with the same inversion eigenvalues; hence the product of the inversion eigenvalues at these two points is $+1$. This proves Eq. (41) in the presence of $C_{4z}$ symmetry.

We now prove Eq. (41) without $C_{4z}$ symmetry using the $\mathbb{Z}_{4}$ invariants, $\chi_{x,y}$. We already proved in Eq. (38) and the surrounding text that the two $\mathbb{Z}_{4}$ invariants have the same parity: $(-1)^{\chi_{x}}=(-1)^{\chi_{y}}$. It is proved in Shiozaki et al., 2016 that the parity of the $\mathbb{Z}_{4}$ invariant, $\chi_{x,y}$, is exactly the $\mathbb{Z}_{2}$ invariant of the $k_{y,x}=0$ plane. Writing the $\mathbb{Z}_{2}$ invariant of the $k_{y,x}=0$ plane in terms of its inversion eigenvalues and equating the parity of $\chi_{x}$ with $\chi_{y}$ yields $\xi_{(0,0,0)}\xi_{(0,0,\pi)}\xi_{(\pi,0,0)}\xi_{(\pi,0,\pi)}=(-1)^{\chi_{x}}=(-1)^{\chi_{y}}=\xi_{(0,0,0)}\xi_{(0,0,\pi)}\xi_{(0,\pi,0)}\xi_{(0,\pi,\pi)}$. By equating the left-most and right-most expressions, we find that $\xi_{(\pi,0,0)}\xi_{(\pi,0,\pi)}=\xi_{(0,\pi,0)}\xi_{(0,\pi,\pi)}$. Hence, the inversion eigenvalues at these points contribute a (trivial) $+1$ to Eq. (40). We already showed above that the two points $(\pi,\pi,0(\pi))$ also contribute a factor of $+1$. Together, this proves Eq. (41).

E.4 Mirror Chern Number in Wallpaper Group $p4g$

The wallpaper group $p4g$ has $C_{4z}$ symmetry, in addition to the two glides, $g_{x,y}$. Consequently, it has the two mirror symmetries $m_{1\bar{1}}\equiv\{m_{1\bar{1}0}|\frac{1}{2}\bar{\frac{1}{2}}0\}$ and $m_{110}\equiv\{m_{110}|\frac{1}{2}\frac{1}{2}0\}$, pictured Fig. 1. Though these symmetries contain translations, they are symmorphic mirrors, and not glides, because their associated translations are along the same axes as their reflections. Equivalently, for a different choice of origin, these symmetries could therefore be written without an accompanying translation. One can define mirror Chern numbers Teo et al. (2008), $n_{1\bar{1}0}$ and $n_{110}$, associated with $m_{1\bar{1}0}$ and $m_{110}$, respectively.

We now show that $(-1)^{n_{1\bar{1}0}}=(-1)^{n_{110}}=(-1)^{\chi_{x}}=(-1)^{\chi_{y}}$ (the last equality was proved above Eq. (38)). We focus on the mirror Chern number $n_{1\bar{1}0}$, associated with $m_{1\bar{1}0}$, which leaves the line $(k,k,k_{z})$ invariant; an identical argument holds for $n_{110}$. As shown in Alexandradinata et al., 2014b,

[TABLE]

where $k_{11(1\bar{1})}\equiv\frac{1}{\sqrt{2}}(k_{x}\pm k_{y})$ and here, $F_{\pm,1\bar{1}}\equiv\partial_{k_{11}}A_{\pm,z}-\partial_{k_{z}}A_{\pm,11}$, where $\pm$ indicates that the trace is over bands with $m_{1\bar{1}0}$ eigenvalue $\pm i$. This quantity is gauge invariant, but we can evaluate it in our gauge choice of Eq. (12). In particular, $A_{\pm,11}(k,k,2\pi)\equiv\langle u^{i}(k,k,2\pi)|\partial_{k_{11}}|u^{j}(k,k,2\pi)\rangle=\langle u^{i}(k,k,0)|V(2\pi\hat{z})\partial_{k_{11}}\left(V(2\pi\hat{z})^{-1}|u^{j}(k,k,0)\rangle\right)=\langle u^{i}(k,k,0)|\partial_{k_{11}}|u^{j}(k,k,0)\rangle\equiv A_{\pm,11}(k,k,0)$, using Eq. (15) and the fact that $\partial_{k_{11}}V(\mathbf{G})=0$. Thus, $\int_{0}^{2\pi}\partial_{k_{z}}A_{\pm,11}=0$, and

[TABLE]

where $\mathcal{W}^{\pm}$ is defined as the Wilson loop evaluated on the bands with $m_{1\bar{1}0}$ eigenvalue $\pm i$. Thus, $n_{1\bar{1}0}$ can be evaluated from a plot of the phases of the eigenvalues of the Wilson loop $\mathcal{W}_{(k,k,0)}$ along the line $\bar{\Gamma}\bar{M}$. in a similar manner to the $\mathbb{Z}_{4}$ invariant calculation in E.1:

Draw a horizontal reference line across the plot. 2. 2.

Count the number of times a positively sloped line in the $+$ sector crosses the horizontal reference line along $\bar{\Gamma}\bar{M}$ and subtract from that the number of times a negatively sloped line in the $+$ sector crosses the horizontal reference line; this gives the $F_{+,1\bar{1}}$ part of Eq. (43), up to edge contributions. 3. 3.

Repeat for the $-$ sector. 4. 4.

$n_{1\bar{1}0}$ is equal to the result from step 2 minus the result from step 3.

This is illustrated in Alexandradinata et al., 2014b and in Fig. 11. From steps 1-4, it is evident that the parity of $n_{1\bar{1}0}$ is equal to the parity of the number of lines crossing the horizontal reference line drawn in step 1, along the segment $\bar{\Gamma}\bar{M}$ (the negative/positive slope no longer matters, since we are considering parity). Now consider the closed loop in the surface Brillouin zone, $\bar{\Gamma}\bar{M}\bar{X}\bar{\Gamma}$. Because the system is insulating, when the Wilson loop, $\mathcal{W}_{(k_{x},k_{y},0)}$, is plotted along this path, an even number of bands must cross the reference line (the assumption that the system is insulating actually places a stronger constraint – that there must be the same number of bands crossing the reference line with positive as with negative slope – but again, since we are considering parity, we only need this constraint mod 2). Since the bands come in pairs along the line $\bar{M}\bar{X}$, the parity of the number of bands crossing the reference line along $\bar{\Gamma}\bar{M}$ is equal to the parity of the number of bands crossing the reference line along $\bar{X}\bar{\Gamma}$; the latter is equal to the parity of $\chi_{y}$, as derived at the end of E.1. It follows that $n_{1\bar{1}0}=\chi_{y}\mod 2$, completing the proof.

Since $(-1)^{\chi_{x}}=(-1)^{\chi_{y}}$ gives us the strong $\mathbb{Z}_{2}$ index, as can be seen from the examples in Fig. 2, which show that the parity of $\chi_{x,y}$ is exactly the parity of the number of linearly dispersing surface states (the Dirac cone at $\bar{M}$ counts twice towards this parity since it comes from four intersecting bands), we conclude that $(-1)^{n_{1\bar{1}0}}$ also provides the strong $\mathbb{Z}_{2}$ index, i.e., if $n_{1\bar{1}0}$ is odd, the occupied bands constitute a strong topological insulating phase.

Appendix F Density Functional Theory Methods and Additional Double-Glide Topological Materials

In this section, we provide additional details on the density functional theory (DFT) calculations in the main text and present supplemental materials results. First, we detail our methodology for generating the bulk electronic band structures for both the materials in the main text and for those analyzed in this supplement. We then provide material-specific descriptions of the $(001)$-directed Wilson loop calculations used to confirm the bulk topology. Finally, we present the details of the semi-infinite surface Green’s function calculations used to predict topological surface features for the materials in the main text. All of the materials candidates presented in this work were previously synthesized, and all were sufficiently stable as to have been powderized and examined under X-ray diffraction in an argon environment Merlo (1984); Bruzzone et al. (1981); Chai and Corbett (2011); Druska et al. (1996); Gumiński (2005); Cordier and Steher (1988).

To identify our topological materials candidates, we first searched the previously calculated electronic structures tabulated in the open online database Materials Project Jain et al. (2013) for double-glide materials in space groups 32, 50, 54, 55, 100, 106, 117, 125, and 127 with small or negligible listed band gaps, zero net magnetic moment, and fewer than 50 atoms per unit cell. Using the most promising candidates from this search, as well as additional materials listed in the inorganic crystal structure database (ICSD) Inorganic Crystal Structure Database (Fachinformationszentrum Karlsruhe, Karlsruhe, Germany, 2015) (ICSD) with fewer than 50 atoms per unit cell, we performed DFT examinations of $\sim 100$ materials. In Tables 1, 2, and 3, we list the chemical formulas and ICSD numbers for the materials which our DFT and Wilson loop calculations determined to be gapless at the Fermi energy or gapped with topologically trivial $\mathbb{Z}_{4}\times\mathbb{Z}_{2}$ glide indices.

Though our investigations only yielded four materials with band gaps at the Fermi energy and nontrivial glide indices, we note that expanding consideration to materials with greater unit cell complexity, the ICSD reports over 2000 previously synthesized materials in SGs 55 and 127 alone. By examining these more complicated materials, or by performing strain or symmetry-preserving chemical substitution, improved topological double-glide materials candidates should be readily discoverable.

Of the two materials presented in the main text, we first detail our analysis of Ba5In2Sb6 in SG 55 Cordier and Steher (1988). We then consider members of the the Si2U3-like A2B3 family of materials in SG 127, which include Sr2Pb3, highlighted in the main text, as well as Au2Y3 and Hg2Sr3. All of the numerical data for the first-principles calculations and figures in this work are freely accessible at https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/EUGQDU.

F.1 Ba5In2Sb6 in SG 55

To explore the electronic properties of Ba5In2Sb6 (ICSD No. 62305), we performed first-principles calculations with the projector-augmented wave (PAW) potential method Blöchl (1994) as implemented in the Vienna ab initio simulation package software (VASP) Kresse and Furthmüller (1996), using the Perdew-Burke-Ernzerhof (PBE)-type generalized gradient approximation (GGA) Perdew et al. (1996) for the exchange-correlation functional. The cutoff energy for the wave-function expansion was set to 400 eV, and the $k$-point sampling grid was set at $6\times 4\times 12$. The experimental lattice parameters Cordier and Steher (1988) were used for Ba5In2Sb6 without structural relaxation, as the band structure near the Fermi level and (thus potentially the bulk topology) was found to be sensitive to the structural parameters, possibly due to the material’s relatively small band gap of 5 meV (indirect) and 17 meV (direct). We additionally confirmed convergency using the above settings.

To calculate the Wilson loop spectrum, we used the following subroutine in the mlwf.f90 module in VASP to calculate the overlap of wavefunctions:

CALL PEAD_CALC_OVERLAP(W,KI,KJ,ISP,P,CQIJ,LATT_CUR,T_INFO,S,LQIJB$=$.TRUE.)

which we then inputted into Eq. (29). The Wilson loop spectrum was calculated for the highest 60 valence bands. To separate the Wilson spectrum into glide sectors, we determined the position-space form of the glide matrices $\hat{g}_{x,y}$ in the plane wave basis:

[TABLE]

where, along the glide-invariant lines:

[TABLE]

We then calculated the Wilson loop individually for the energy eigenstates within each sector of $\hat{g}_{x,y}$. We determined the bulk topology of Ba5In2Sb6 by utilizing the rules in E.1. In Fig. 5(d) we draw a dashed green line across the Wilson spectrum and count the glide-polarized Wilson bands that cross it. Calculating the eigenvalues of $g_{y}$, where valid, along the path $\bar{M}\bar{Y}\bar{\Gamma}\bar{X}$, we count no bands in the $(+)$ sector along $\bar{M}\bar{Y}$, one positively and one negatively sloped band along $\bar{Y}\bar{\Gamma}$, and no bands in the $(+)$ sector along $\bar{\Gamma}\bar{X}$. This results in $\chi_{y}=0+1-1+0\mod 4=0$. To obtain $\chi_{x}$, we perform a similar calculation using the Wilson spectrum along $\bar{M}\bar{X}\bar{\Gamma}\bar{Y}$ labeled by the eigenvalues of $g_{x}$ along the two lines where that symmetry is valid. We count no bands in the $(+)$ sector along $\bar{M}\bar{X}$, no bands along $\bar{X}\bar{\Gamma}$, and one negatively sloped band in the $(+)$ sector along $\bar{\Gamma}\bar{Y}$ (the slope here is taken to be negative because we are moving right to left along this line in Fig. 5(d)). This results in $\chi_{x}=0+0+(0-1)\times 2\mod 4=2$, giving an overall bulk topology of $(\chi_{x},\chi_{y})=(2,0)$. The projected surface states were obtained from the surface Green’s function of a semi-infinite system. For this purpose, we constructed maximally localized Wannier functions for the $s$ orbitals of Ba and the $p$ orbitals of Sb from first-principles calculations using Wannier90 Marzari et al. (2012); Marzari and Vanderbilt (1997); Mostofi et al. (2008); Souza et al. (2001).

F.2 Sr2Pb3, Au2Y3, and Hg2Sr3 in SG 127

We find that three previously-synthesized members of the Si2U3 family of materials in SG 127, Sr2Pb3 (ICSD No. 105627) Merlo (1984); Bruzzone et al. (1981), Au2Y3 (ICSD No. 262043) Chai and Corbett (2011), and Hg2Sr3 (ICSD Nos. 107371, 247135) Druska et al. (1996); Gumiński (2005), are capable of hosting double-glide topological crystalline phases. To explore the electronic properties of these materials, we employed first-principles calculations based on DFT. The PBE-type GGA Perdew et al. (1996) was used to describe exchange correlation as implemented in the Quantum Espresso package Giannozzi et al. (2009). Core electrons were treated by norm-conserving, optimized, designed nonlocal pseudopotentials, generated using the Opium package Rappe et al. (1990); Ramer and Rappe (1999). We used the energy threshold of 680 eV for the plane wave basis. Relativistic effects of spin-orbit interaction were fully described using a non–collinear scheme. The atomic structures were initially obtained from the ICSD Inorganic Crystal Structure Database (Fachinformationszentrum Karlsruhe, Karlsruhe, Germany, 2015) (ICSD), and then fully relaxed to achieve consistency between the electronic and structural states using a force threshold of 0.005 eV/Å. The lattice constants of Sr2Pb3 changed negligibly after the full structural relaxation from a (c) = 8.367 (4.883) Å to 8.383 (4.944) Å. We calculated the Wilson bands of Sr2Pb3 within each glide sector following the methodology detailed in F.1. We found that the bulk topology was unaffected by this structural relaxation. The Wilson bands of Sr2Pb3 were calculated using 20 bands from the $N-19$-th to the $N$-th bulk bands, where N is the number of the valence electrons per unit cell. We determined the bulk topology of Sr2Pb3 by utilizing the rules in E.1. In Fig. 3(d) we draw a dashed green line across the Wilson spectrum and count the glide-polarized Wilson bands that cross it. Calculating the eigenvalues of $g_{y}$, where valid, along the path $\bar{M}\bar{X}^{\prime}\bar{\Gamma}\bar{X}$, we count one positively sloped $(+)$ band along $\bar{M}\bar{X}^{\prime}$, one positively and one negatively sloped band along $\bar{X}^{\prime}\bar{\Gamma}$, and no bands in the $(+)$ sector along $\bar{\Gamma}\bar{X}$. This results in $\chi_{y}=(1\times 2)+1-1+0\mod 4=2$. By the $C_{4z}$ symmetry of SG 127, $\chi_{x}=\chi_{y}$, giving an overall bulk topology of $(\chi_{x},\chi_{y})=(2,2)$, that of a nonsymmorphic Dirac insulator. The surface states of Sr2Pb3 were obtained by calculating the surface Green’s function for a semi-infinite system Sancho et al. (1984, 1985) based on the maximally localized Wannier functions for the $d$ orbitals of Sr and the $p$ orbitals of Pb using Wannier90 Marzari et al. (2012); Marzari and Vanderbilt (1997); Mostofi et al. (2008); Souza et al. (2001).

As detailed in the main text and in Fig. 12, all three of these materials share the same unit cell structure, with two orthogonal glide mirrors in the $x$ and $y$ directions related by $C_{4z}$ symmetry and a mirror in the $z$ direction that relates the upper and lower layers of the unit cell. We find that Sr2Pb3 posses a gap at the Fermi energy at each crystal momentum but that, as shown in Fig. 12(d,f), pristine Au2Y3 and Hg2Sr3 are gapless, with $C_{4z}$-protected Dirac points near the Fermi energy. Applying strain in the $(100)$-direction ($x$) breaks $C_{4z}$ while still preserving the two glides that project to form wallpaper group $pgg$ on the $(001)$-surface. As this lower-symmetry bulk system still possess inversion and time-reversal symmetries, Weyl points are forbidden Vafek and Vishwanath (2014), and therefore under $C_{4z}$-breaking strain the Dirac points in Au2Y3 and Hg2Sr3 become fully gapped. We applied incremental strain and subsequent relaxation of the internal atomic coordinates and evaluated the symmetries of the strained systems using FINDSYM Stokes and Hatch (2005) with a tolerance of 0.001 Å. We find that Au2Y3 develops a gap for up to 5% strain, whereas Hg2Sr3 develops a gap for up to around 4% strain, after which it becomes unstable and undergoes a structural transition to the symmorphic SG 75.

For each of these materials, we therefore evaluate the $(001)$-directed Wilson loop under 2% strain, shown in Fig. 13, where the Wilson bands within each glide sector are calculated following the methodology detailed in F.1. The Wilson bands for Au2Y3 are calculated using the 16 bands from the $N-17$-th to the $N-2$-th bands, valid as the band gap between the $N-2$-th and the $N-1$-th bands widely overlaps with the Fermi level. For Hg2Sr3, we use the 12 bands from band indices $N-11$ to $N$. Though the Wilson spectra for these materials appear to be more complicated than those of Ba5In2Sb6 and Sr2Pb3, we can still use the simple rules developed in E.1 to evaluate the $\mathbb{Z}_{4}$ indices.

We determine the bulk topology of these two strained materials by following the rules detailed in E.1 and used previously in F.1 to evaluate the topology of Ba5In2Sb6 in SG 55. For 2%-strained Au2Y3, we draw a dashed green line across the Wilson spectrum in Fig. 13(a) and count the glide-polarized Wilson bands that cross it. Calculating the eigenvalues of $g_{y}$, where valid, along the path $\bar{M}\bar{Y}\bar{\Gamma}\bar{X}$, we count one positively and one negatively sloped band in the $(+)$ sector along $\bar{M}\bar{Y}$, two positively sloped bands along $\bar{Y}\bar{\Gamma}$, and two negatively sloped bands in the $(+)$ sector along $\bar{\Gamma}\bar{X}$. This results in $\chi_{y}=(1-1)\times 2+2+(0-2)\times 2\mod 4=2$. To obtain $\chi_{x}$, we perform a similar calculation using the Wilson spectrum along $\bar{M}\bar{X}\bar{\Gamma}\bar{Y}$ labeled by the eigenvalues of $g_{x}$ where applicable. We note that as we are now taking a right-to-left path in Fig. 13(a), the positive or negative Wilson band slope labeling will here be the opposite of what it was in the calculation of $\chi_{y}$. We count one positively and one negatively sloped band in the $(+)$ sector along $\bar{M}\bar{X}$, two positively sloped bands along $\bar{X}\bar{\Gamma}$, and one negatively sloped band in the $(+)$ sector along $\bar{\Gamma}\bar{Y}$. This results in $\chi_{x}=(1-1)\times 2+2+(0-1)\times 2\mod 4=0$, giving an overall bulk topology of $(\chi_{x},\chi_{y})=(0,2)$.

For 2%-strained Hg2Sr3, we draw a dashed green line across the Wilson spectrum in Fig. 13(b) and count the glide-polarized Wilson bands that cross it. Calculating the eigenvalues of $g_{y}$, where valid, along the path $\bar{M}\bar{Y}\bar{\Gamma}\bar{X}$, we count no bands in the $(+)$ sector along $\bar{M}\bar{Y}$, two positively sloped bands along $\bar{Y}\bar{\Gamma}$, and two positively and two negatively sloped bands in the $(+)$ sector along $\bar{\Gamma}\bar{X}$ (one of the $(+)$ bands has a very sharp slope, and its continuity across the green line can be inferred by tracing where it appears again at $\theta=-\pi$). This results in $\chi_{y}=0+2+(2-2)\times 2\mod 4=2$. To obtain $\chi_{x}$, we perform a similar calculation using the Wilson spectrum along $\bar{M}\bar{X}\bar{\Gamma}\bar{Y}$ labeled by the eigenvalues of $g_{x}$ where applicable. We note that as we are now taking a right-to-left path in Fig. 13(b), the positive or negative Wilson band slope labeling will here be the opposite of what it was in the calculation of $\chi_{y}$. We count no bands in the $(+)$ sector along $\bar{M}\bar{X}$, four positively and two negatively sloped bands along $\bar{X}\bar{\Gamma}$, and one negatively sloped band in the $(+)$ sector along $\bar{\Gamma}\bar{Y}$. This results in $\chi_{x}=0+4-2+(0-1)\times 2\mod 4=0$, giving an overall bulk topology of $(\chi_{x},\chi_{y})=(0,2)$.

To summarize, we find that both 2% strained Au2Y3 and 2% strained Hg2Sr3 host the $(\chi_{x},\chi_{y})=(0,2)$ double-glide topological hourglass connectivity.

Appendix G Tight-Binding Model and the SSH Limit

In this section, we present a simplified tight-binding model that can realize all of the $\mathbb{Z}_{4}\times\mathbb{Z}_{2}$ insulating phases allowed for wallpaper groups $pgg$ and $p4g$, which have perpendicular glides, $g_{x,y}$, in the $x$ and $y$ directions. For simplicity, we further specialize to systems with inversion symmetry, $\mathcal{I}$; as shown in E.3, this simplifies the computation of the $z$-projected Wilson loop. Inversion symmetry also implies the presence of a mirror in the $z$ direction: $\mathcal{I}=g_{x}g_{y}m_{z}$. We find a fine-tuned limit in which the $\mathbb{Z}_{4}$ topological invariants defined in E.1 can be computed by comparing the relative values of three Su-Schrieffer-Heeger (SSH) $\mathbb{Z}_{2}$ invariants Su et al. (1979) defined by effective 1D chains at the corners of the BZ which project to $\bar{X}$, $\bar{Y}$, and $\bar{M}$ on the $z$-normal surface. This limit provides an alternative, intuitive way to understand the four topological phases where $\chi_{x}=\chi_{y}=0\mod 2$. We note that the Wilson spectrum pinning in this tight-binding model can alternatively be understood in terms of the bulk inversion eigenvalues as detailed in Alexandradinata et al., 2014a.

G.1 Tight-Binding Model for Space Groups 55 and 127

As a starting point, consider a single layer of the two-site unit cell shown in Fig. 1. The sites, designated $A$ and $B$, are related to each other by glide reflections, $g_{x,y}=\{m_{x,y}|\frac{1}{2}\frac{1}{2}0\}$; the origin is defined to sit at an $A$ site. An orthorhombic stack of this single layer, with no other symmetries, is in space group 32, $Pba2$, and its $z$-normal $(001)$ surface is characterized by the two-dimensional wallpaper group, $pgg$.

We enforce an extra mirror symmetry, $m_{z}=\{m_{z}|00\frac{1}{2}\}$, by adding two more sublattices, $C$ and $D$, sitting $t_{z/2}$ above the $A$ and $B$ sites respectively (Fig. 14(a)); this mirror also adds an inversion symmetry $\mathcal{I}=g_{x}g_{y}m_{z}$. The resulting system is in space group 55, $Pbam$. In this space group, when the two layers are decoupled, each sheet is equivalent to a single layer of the two-dimensional Dirac semimetal model in Young and Kane, 2015, which possesses fourfold degeneracies at $X$, $Y$, and $M$ due to the algebraic relations between $\mathcal{I}$, $g_{x}$, and $g_{y}$ at those TRIMs. For our 3D orthorhombic system, at all six bulk TRIMs for which $k_{x}$ or $k_{y}$ is equal to $\pi$, states are fourfold-degenerate by the algebra from A.2 and Eq. (10). At the four TRIMs for which $(k_{x},k_{y})=(0,\pi)$ or $(\pi,0)$, at least one of the glides anticommutes with inversion, which, combined with $\mathcal{I}^{2}=-\mathcal{T}^{2}=+1$, requires states to be fourfold-degenerate. At the two TRIMs for which $(k_{x},k_{y})=(\pi,\pi)$, the two glides anticommute and square to $+1$, which, combined with $\mathcal{T}^{2}=-1$, also enforces a fourfold degeneracy.

The following Hamiltonian includes all allowed in-plane hopping terms up to second-nearest-neighbor:

[TABLE]

where $\tau^{x}$ corresponds to hopping between the $A$ and $B$ (or $C$ and $D$) orbitals, $\mu^{x}$ corresponds to hopping between the $A$ and $C$ (or $B$ and $D$), orbitals and the $\sigma$ Pauli matrices correspond to an on-site spin.

To further simplify, we impose $C_{4z}$ symmetry, which is implemented by the operator, $C_{4z}=\sqrt{i\sigma^{z}}f_{4z}(\vec{k})$, where $f_{4z}(\vec{k})$ acts on the crystal momenta by enforcing the cyclical mapping:

[TABLE]

With this additional symmetry, $v_{s1}=-v_{s1}^{\prime},\ v_{s2}=v_{s2}^{\prime\prime\prime},\ v_{s2}^{\prime}=-v_{s2}^{\prime\prime}$ in Eq. (47) and the system is now in the higher-symmetry space group 127, $P4/mbm$, with an in-plane Hamiltonian up to second-nearest neighbor hopping given by Eq. (47) with these restrictions.

We now introduce hopping in the $z$-direction to couple the layers. We begin in a relatively artificial limit by only turning on specific SOC-free dimerizing terms, (Fig. 14(b)):

[TABLE]

and

[TABLE]

The terms proportional to $u_{1,2}$ correspond to hopping between nearest-neighbor $A$ ($B$) and $C$ ($D$) sites. Terms proportional to $v_{1,2}$ originate from longer-range versions of the same type of hopping, and connect $A$ ($B$) and $C$ ($D$) sites separated by $\vec{d}=\{10\frac{1}{2}\}$ and $\vec{d}=\{01\frac{1}{2}\}$.

We now consider the fine-tuned limit where $v_{r1}=v_{s1}=0$, or one in which there is no spin-orbit coupling at $(k_{x},k_{y})=(\pi,0)$ and $(0,\pi)$. In this limit, we can write down the bulk Hamiltonian on the line that projects to $\bar{X}$ on the $z$-normal surface:

[TABLE]

where we have shifted the energy to make $t_{2}=0$ without loss of generality. We note that the $\tau$ sublattice and $\sigma$ spin degrees of freedom no longer play a role, they merely impose additional degeneracies.

This limit places strong constraints on the Wilson loop bands: since the spins are decoupled, we can consider the spinless time-reversal operator, $\tilde{\mathcal{T}}$, which satisfies $\tilde{\mathcal{T}}^{2}=+1$ and the spinless glide, $\tilde{g}_{y}$, which, at the $\bar{X}$ point, satisfies $\tilde{g}_{y}^{2}=-1$. Thus, at $\bar{X}$, $\tilde{\mathcal{T}}\tilde{g}_{y}$ is an antiunitary operator that squares to $-1$, enforcing that all eigenstates are doubly degenerate, within each spin sector. Since there is no SOC, the two spin sectors are also degenerate, resulting in a fourfold degeneracy at the $\bar{X}$ point, though one which is broken into twofold degeneracies when SOC is realistically reintroduced. In addition, inversion requires that the Wilson bands are particle-hole symmetric Alexandradinata et al. (2016). Thus, our four-band model has all four Wilson bands degenerate at the $\bar{X}$ point and, because of inversion symmetry, they are pinned to either [math] or $\pi$.

It was shown in D.3 that all four Wilson bands are degenerate at $\bar{M}$, as well, and by inversion are also pinned to either [math] or $\pi$. Since the Wilson bands must continuously connect the bands at $\bar{X}$ to the bands at $\bar{M}$, there is a $\mathbb{Z}_{2}$ invariant that characterizes the possible connectivities for fixed band inversion at $\bar{\Gamma}$ (Fig. 15):

[TABLE]

One of these connectivities is a trivial phase and the other is the nonsymmorphic Dirac insulating phase. In Fig. 16, we plot the bulk and Wilson bands for our model in this limit and demonstrate both trivial and nonsymmorphic Dirac insulating phases. As long as the system is QSH-trivial and diagonal mirror Chern-trivial, then the Wilson loop eigenvalues are generically unpinned at $\bar{\Gamma}$ and along all low symmetry lines (though, as seen in Fig. 16(e), in the non-generic case that all of bulk inversion eigenvalues are the same at the $\bar{\Gamma}$-projecting TRIMs, the Wilson loop will additionally be pinned at $\bar{\Gamma}$ by the mechanism detailed in Alexandradinata et al., 2014a). Therefore, in this particular model of SG 127, all of the possible crystalline connectivities for the path $\bar{\Gamma}\bar{X}\bar{M}\bar{X}^{\prime}$ are entirely determined by how the eigenvalues are pinned at $\bar{X}$ and $\bar{M}$.

This limit can also be described by two, spin-degenerate Su-Schrieffer-Heeger (SSH) chains: the layer degree of freedom $\mu$ acts like the sublattice degree of freedom in the original, two-band SSH model, and each chain in this limit has an additional degeneracy in the spin subspace $\sigma$. The edge state of one chain is the projection of the Hamiltonian onto $\bar{X}$, and the edge state of the other is the projection onto $\bar{M}$. The relative values of $u_{1,2}$ and $v_{1,2}$, which parameterize hopping in the $z$-direction, correspond to the two choices of dimerization for each SSH chain. For each of these two chains, there is a $\mathbb{Z}_{2}$ polarization, $\theta=0,\pi$, which directly corresponds to the Wilson phases at those surface TRIMs, and the overall bulk topology for these crystalline phases is given by the $\mathbb{Z}_{2}$ relative polarization between the two SSH models.

For the chain which projects to $\bar{M}$, the fourfold surface degeneracy prevents the two Hamiltonians from coupling. However, for the chain projecting to $\bar{X}$, the two, spin-degenerate SSH Hamiltonians are only prevented from coupling in the limit that there is no SOC at any bulk $k$-point projecting to $\bar{X}$. As symmetry-allowed SOC terms are turned on, a real system will escape this limit, and the two copies of the SSH Hamiltonian at $\bar{X}$ will couple and their surface states will gap into pairs. Nevertheless, as discussed in further detail in the next section, if introducing SOC does not result in a bulk gap closure, then the value of $\chi$ cannot change, and the bulk topology will remain the same as it was in the SSH limit.

G.2 Beyond the SSH Limit

Generically, all symmetry-allowed hopping terms will be present in a real material, including $v_{r1}$, $v_{s1}$, and the other $z$-direction hopping terms. In this section, we examine how this affects the SSH-model definition of crystalline invariants.

The line that projects to $\bar{M}$ continues to be described by a doubly degenerate SSH chain even after generic terms are added, as it still hosts a fourfold-degenerate surface state with a $\mathbb{Z}_{2}$ polarization. As states at $\bar{M}$ are fourfold-degenerate due to surface wallpaper symmetries, the Wilson phases $\theta(\bar{M})$ must still be either [math] or $\pi$ as long as there is bulk inversion symmetry. Furthermore, in the limit of weak spin-orbit interaction, a band inversion at an $\bar{M}$-projecting TRIM will result in the presence of bulk Dirac point or line nodes. However, the properties of these line nodes and their relationship to nonsymmorphic symmetries are space-group dependent, and in general go beyond the focus of this manuscript.

At the TRIMs which project onto $\bar{X}$, the immediate consequences of allowing nonzero values of $v_{r1}$ and $v_{s1}$ are to couple the two copies of the SSH Hamiltonians under $\bar{X}$ and gap out their surface states. As states at $\bar{X}$ are only generically twofold-degenerate, the presence of additional SOC terms breaks the artificial symmetries $\tilde{g}_{y}$ and $\tilde{\mathcal{T}}$ and gaps the states into the two ends of an hourglass. Nevertheless, as pictured in Fig. 17, this hourglass along $\bar{\Gamma}\bar{X}$ is still characterized by a $\mathbb{Z}_{2}$ invariant that characterizes whether it is centered about a Wilson phase of [math] or $\pi$. Therefore, as long as there are no additional band inversions, then the Wilson connectivity remains unchanged and, away from the SSH limit there still is an overall $\mathbb{Z}_{2}$ quantity that characterizes the topological crystalline phases. However, as the phase of the hourglass center can be moved by any band inversion on the plane which projects to $\bar{\Gamma}\bar{X}$, then we find that this $\mathbb{Z}_{2}$ invariant is no longer just a property of the TRIMs. Therefore, for a real double-glide system, the Wilson loop remains the only generic method for evaluating the bulk topology.

G.3 Broken $C_{4z}$ Phases in and beyond the SSH Limit

Additional Wilson band topologies are possible in systems without $C_{4z}$ symmetry, whose surfaces are described by the wallpaper group $pgg$. As described in the main text and in F, Ba5In2Sb6, as well as strained Au2Y3 and Hg2Sr3, are expected to exhibit such phases. We break $C_{4z}$ symmetry in the tight-binding model by adding another interlayer hopping term:

[TABLE]

The resulting system is now in SG 55. Without $C_{4z}$ symmetry, the SSH polarizations at $\bar{X}$ and $\bar{Y}$ can now differ, leading us to define a second crystalline invariant for the independent $y$ direction:

[TABLE]

The reason that $\chi_{x,y}$ is determined by $\theta(\bar{Y},\bar{X})$, as opposed to $\theta(\bar{X},\bar{Y})$, is because $g_{y}$ enforces a fourfold degeneracy at $\bar{X}$, as explained in G.1. Thus, $\chi_{x,y}$ is determined by $g_{x,y}$, consistent with the notation in E.2.

If the polarization at $\bar{X}$ differs from that at $\bar{Y}$ and $\bar{M}$, the system can display a single fourfold point at $\bar{X}$. As more symmetry-allowed SOC terms are added, the states at $\bar{X}$ will couple and the four-band crossing will open up into an hourglass along $\bar{\Gamma}\bar{X}$. Nevertheless, as long as the bulk symmetries remain unchanged and the introduction of SOC preserves the bulk gap and band ordering, the Wilson connectivity will be preserved and the $C_{4z}$-broken SSH phase will evolve into a crystalline hourglass phase. An observation of the allowed Wilson band connectivities with broken $C_{4z}$ symmetry confirms in fact that these two phases (hourglasses along $\bar{\Gamma}\bar{X}$ or along $\bar{\Gamma}\bar{Y}$) are the only two allowed topological crystalline connectivities which fundamentally violate $C_{4z}$. The $(\chi_{x},\chi_{y})=(0,2)$ phase is demonstrated in our tight-binding model in Figure 18.

If symmetry-allowed terms are made larger, other band inversions may occur and $\chi_{x/y}$ may eventually change. However, the resulting phases, if QSH-trivial, will still ultimately be tunable back to an SSH limit. Therefore the SSH limit can be considered as a parent phase to all four possible $\chi_{x,y}=0,2$ crystalline insulating phases.

G.4 Strong Topological Insulating Phases

Away from the SSH limit, but preserving $C_{4z}$ symmetry such that the system is in SG 127, as shown in E.3, a band inversion about a TRIM which projects to $\bar{\Gamma}$ can flip the product of parity eigenvalues and induce a strong topological insulating phase. In practice, this band inversion can be accomplished in our tight-binding model by adding the term:

[TABLE]

to the Hamiltonian in Eq. (50). Counting the product of the parity eigenvalues at the bulk TRIMs (Fig. 19(c)) and grouping them into their projections to $z$-normal TRIMs:

[TABLE]

we confirm that this system is a strong topological insulator with $\xi(\bar{\Gamma})=-1$ and $\xi(\bar{X})=\xi(\bar{X}^{\prime})=\xi(\bar{M})=+1$. The Wilson loop in Fig. 19(b) further shows that this system is in a double-glide spin Hall phase characterized by $\chi_{x}=\chi_{y}=1,3$, depending on the choice of an odd- or even-numbered surface BZ.