Dislocation as a bulk probe of higher-order topological insulators

Topological materials occupy the central stage in the modern condensed matter physics because of their robust metallic edge or surface states protected by the topological invariant, characterizing the electronic band structure in the bulk. Higher-order topological (HOT) states extend this usual bulk-boundary correspondence, so they host the modes localized at lower-dimensional boundaries, such as corners and hinges. Here we theoretically demonstrate that dislocations, ubiquitous defects in crystalline materials, can probe higher-order topology, recently realized in various platforms. We uncover that HOT insulators respond to dislocations through symmetry protected finite-energy in-gap electronic modes, localized at the defect core, which originate from an interplay between the orientation of the HOT mass domain wall and the Burgers vector of the dislocation. As such, these modes become gapless only when the Burgers vector points toward lower-dimensional gapless boundaries. Our findings are consequential for the systematic probing of the extended bulk-boundary correspondence in a broad range of HOT crystals, and photonic and phononic or mechanical metamaterials through the bulk topological lattice defects.

Topological materials occupy the central stage in the modern condensed matter physics because of their robust metallic edge or surface states protected by the topological invariant, characterizing the electronic band structure in the bulk. Higher-order topological (HOT) states extend this usual bulkboundary correspondence, so they host the modes localized at lower-dimensional boundaries, such as corners and hinges. Here we theoretically demonstrate that dislocations, ubiquitous defects in crystalline materials, can probe higher-order topology, recently realized in various platforms. We uncover that HOT insulators respond to dislocations through symmetry protected finite-energy in-gap electronic modes, localized at the defect core, which originate from an interplay between the orientation of the HOT mass domain wall and the Burgers vector of the dislocation. As such, these modes become gapless only when the Burgers vector points toward lower-dimensional gapless boundaries. Our findings are consequential for the systematic probing of the extended bulk-boundary correspondence in a broad range of HOT crystals, and photonic and phononic or mechanical metamaterials through the bulk topological lattice defects.

I. INTRODUCTION
The nontrivial topological invariant characterizing the bulk electronic band structure gives rise to robust edge or surface modes, manifesting the hallmark of a topological material -the bulk-boundary correspondence [1,2]. As such, these boundary modes have been so far almost exclusively used to experimentally detect nontrivial electronic topology, both in gapped [3][4][5][6][7][8] and gapless [9,10] systems. Equally important, but much less explored, is the direct probing of topological states in the bulk without invoking the boundary modes, through their response to topological lattice defects, such as dislocations [11][12][13][14][15][16][17][18][19][20]. Moreover, the topological defect modes are more pristine, being immune to contamination by the interfaces and independent of the surface termination. In fact, in the context of experimental probing of topology in the quantum materials this aspect has started to gain prominence only recently [21,22].
In D−dimensional nth order topological states [23][24][25][26][27], bulk probing of the electronic band topology should play an important role, because the extended bulk-boundary correspondence is realized through gapless modes on the lower, (D −n)−dimensional boundaries, characterized by codimension d c = n, such as hinges (d c = D − 1) and corners (d c = D) [28][29][30][31][32][33][34][35]. Their robustness originates from the combination of spatial symmetries, such as discrete rotations, and non-spatial ones, such as the reversal of time. Importantly, these protected modes on lowerdimensional boundaries may be thought of as inherited from the parent, first-order topological state (with n = 1) upon partially gapping out its edge or surface modes * bitan.roy@lehigh.edu † vladimir.juricic@su.se (d c = 1), which can, in turn, yield a hierarchical ladder of HOT states [36,37]. This is accomplished by a suitable domain wall mass which changes sign across corners [see Fig. 1(a)] or hinges, thus localizing topological modes at these lower-dimensional boundaries [see Figs. 1(b), 2, 3(a) and 3(c)]. The reduced dimensionality of the boundary may, however, hinder the experimental detection of the gapless modes, and therefore HOT states require other means to directly probe the bulk electronic topology.
As we demonstrate here, dislocations can serve as bulk probes of HOT insulators through the binding of special topologically protected electronic modes, see Figs. 1 and 3. To formulate the mechanism, we recall the Volterra construction: in a two-dimensional (2D) lattice a dislocation can be created by removing a line ending at the dislocation center (Volterra cut), and reconnecting the sites across this cut so that the translational symmetry is restored away from the defect center (core), see Fig. 1(a). Therefore, any closed loop around the dislocation center features a missing translation by the Burgers vector b, which topologically characterizes the defect. As such, a dislocation provides global frustration to the underlying crystalline order, which translates into a nontrivial effect on the electrons hopping on the lattice. Namely, an electron with a momentum K when encircling the dislocation picks up a phase equal to exp[iΦ dis ], with Φ dis = K · b (mod 2π). In particular, for topological states with the band inversion momentum at K inv , the hopping phase is Φ dis = K inv · b (mod 2π) [11].
To set the stage, recall that a translationally-active first-order topological insulator features at least one band inversion at a finite (non-Γ) momentum in the Brillouin zone (BZ) [15] yielding gapless edge states. When Φ dis = π in a translationally-active topological insulator, after encircling a dislocation, the electrons pick up a hop- The defect is obtained through the Volterra cut-and-glue procedure by removing a line of atoms ending at the center of the lattice (orange) and reconnecting the edges across this Volterra cut, which are denoted by + and −, right and left from the center, respectively. The corresponding Burgers vector is b = −aex. The HOT mass domain walls along (across) which it vanishes (changes sign) are represented by the red and blue dashed lines for θ = 0 and π/2, respectively, as given by Eq. (1). (b) Local density of states (LDoS) for the dislocation mode localized at the defect core together with the four corner modes in a second-order translationally-active HOT insulator with the band inversion at the M point of the BZ. Here, we set θ = 0 in the HOT mass, so that the defect modes are at finite energies, while the corner modes are at zero energy. (c) LDoS for the zero-energy dislocation modes in the periodic system with a dislocation-antidislocation pair for θ = π/2. Here, we set t = 2B = 1, m = 3 and ∆ = 0.20 [see Eq. (2)]. Any site with LDoS less than 10 −3 is left empty. See also Appendix D.
ping phase equal to exp(iπ) = −1 across the Volterra cut, see Fig. 1(a). In turn, to resolve the frustration in the hopping introduced by the defect through the nontrivial phase factor, a Kramers pair of zero-energy states gets localized at the dislocation core [13].
In contrast, a second-order topological insulator features gapped edge states that stem from the mass domain wall in the bulk [Figs. (1)(a) and 4]. The domain wall mass gaps out the edge states but only partially, in turn producing the topological corner modes through the Jackiw-Rebbi mechanism [38]. Now, when a dislocation is inserted, the defect modes, as we show, still survive [Figs. 1(b),(c)], but, are moved away from zero energy, since the edge states across the Volterra cut are gapped [ Fig. 2(a)]. Importantly, when the orientation of the Burgers vector (b) is parallel to the direction of the domain wall of the HOT mass, the defect modes are pinned at zero energy. Hierarchy of the HOT states in this way directly translates into the spectral flow of the dislocation modes, detectable in the tunneling spectroscopy measurements, for instance. The same mechanism is analogously operative for three-dimensional (3D) second-order translationally-active insulators: an edge dislocation hosts gapped modes which become gapless only when the Burgers vector is parallel to the HOT mass domain wall. Furthermore, a screw dislocation hosts gapless propagating modes only when it is orthogonal to a gapless surface or equivalently parallel to the corresponding surface normal [ Fig. 3 and Fig. 5]. Finally, we emphasize that the composite C 4 T , PT and C 4 P symmetries protecting the zero-energy corner (hinge) modes, also protect both finite-and zero-energy dislocation modes, displayed in Figs. 1(b), 3(a), and 3(c). Here C 4 , T and P represent discrete four-fold rotational, time-reversal and parity symmetries, respectively. For details consult Appendix E, where it is also shown that this protection mechanism for dislocation modes extends to C 4n rotational symmetry breaking HOT insulators in two and three dimensions, with n > 1. Furthermore, in Appendix F we show that the defect modes are protected also in the case of C 4n+2 rotational symmetry breaking HOT insulators, where n ≥ 1.
The rest of the paper is organized as follows. In Sec. II, we discuss the universal tight-binding model for the HOT insulators in d = 2 and d = 3. Section III is devoted to the numerical results for the dislocation modes on a square lattice. In Sec. IV we present a general argument for the existence of the dislocation modes in a HOT insulator, and in Sec. V we show the numerical results for the dislocation modes in 3D HOT insulators. We discuss the results in Sec. VI and highlight their possible realizations in HOT crystals and metamaterials. Additional technical details are relegated to the Appendices.

II. TIGHT-BINDING MODEL
To show the outlined general mechanism, we take the minimal, but the universal tight-binding model describ- ing a second-order topological insulator in d = 2 and d = 3 [36,39], with the Hamiltonian H = k Ψ † kĥ Ψ k , whereĥ =ĥ 0 +ĥ ∆ , and Here (Γ, Γ d+1 , Γ d+2 ) are the mutually anticommuting four-component Γ matrices, k is the momentum, and a is the lattice spacing. The above Hamiltonian breaks C 4 rotational symmetry about the z−axis, generated by R 4 = exp(iπΓ 12 /4), where Γ 12 = [Γ 1 , Γ 2 ]/(2i), the timereversal and the parity symmetries with (representationdependent) operators T and P, respectively, but preserves their products C 4 T , T P and C 4 P. Under fourfold rotation k x → −k y and k y → k x . It should be noted that with the above form of the generator of rotations, following the Lie group, the HOT mass always breaks discrete rotational symmetry. In both d = 2 and d = 3, we take the form factors d i (k) = t sin k i a, while the first-order mass is given by We consider translationally-active M phase in d = 2 with a band inversion at M = (π/a, π/a) point in the BZ, which is realized in the parameter range 4 < m/B < 8. In d = 3 we take translationally-active R phase with the band inversion at R point in the BZ, R = (π/a, π/a, π/a), for 8 < m/B < 12. The resulting first-order topological insulator supports edge and surface states, both with d c = 1, respectively, in d = 2 and d = 3. Furthermore, h ∆ acts as a mass term for the topological edge (surface) states, and leaves only the corners (hinges) gapless, yielding a second-order topological insulator with corner (1) we fix this Wilson-Dirac mass term so that it changes sign under the C 4 rotation, transforming k x → −k y , k y → k x . As such, this mass term for each value of the parameter 0 ≤ θ ≤ π/2 necessarily features a line across which it changes sign. In particular, for θ = π/2 the domain wall lies along the principal axes, k x = 0, k y = 0, while for θ = 0 it is located along the diagonals k y = ±k x . See Fig. 1(a) for the HOT mass domain walls in the real space.

III. 2D LATTICE DISLOCATION MODES
We perform numerical analysis of the translationallyactive HOT M phase, hosting a band inversion at the M = (π/a, π/a) point in the BZ. The implementation of the model, given by Eq. (1), was carried out in the real space on a square lattice hosting a dislocation defect with the Burgers vector b = ae x , oriented in the lattice x−direction, as shown in Fig. 1(a). See also Appendix D for details. In an open system, the dislocation modes and the corner states coexist in the HOT M phase, explicitly showing that the defect can probe the extended bulkboundary correspondence, see Fig. 1(b). Furthermore, we find that dislocations bind the modes in a lattice without boundaries (periodic system), further corroborating their role as a pure bulk probe of higher-order electronic topology, see Fig. 1(c). The hybridization effects in both cases can be neglected, as the defect modes are localized within a few lattice sites around its center, which is much shorter than both the system size and the separation between the defects.
Most importantly, for any choice of the domain wall orientation (θ), we find that the defects feature mid-gap bound states at finite energies, see Fig. 2(a). As the domain wall orientation approaches the direction of the Burgers vector (θ → π/2), the spectral gap (δE) between the dislocation modes decreases. Eventually, when the two directions coincide (θ = π/2), the modes become degenerate zero energy states. For the spectral flow of the dislocation modes, see Fig. 6. Also notice that finite energy dislocation modes are particle-hole partners, while they become eigenmodes of the particle-hole operator when pinned at zero energy. See Appendix C for details. The scaling of the gap with the domain wall orientation for various amplitudes of the Wilson-Dirac mass (∆) is displayed in Fig. 2(a), showing that δE → 0 as θ → π/2. We next present a general argument supporting this observation.

IV. DISLOCATION MODES: A GENERAL ARGUMENT
To this end, we recall that in the parent first-order topological insulator before reconnecting the edges across the Volterra cut, each edge (1) features a Kramers pair of zero energy (due to a unitary particle-hole symmetry, see Appendix C) helical modes, and (2) is perpendicular to the Burgers vector b. The zero energy states at each of the edges are then the eigenstates of the matrix and Γ 3 is the mass matrix for the first-order topological insulator in d = 2 [see Eq. (1) and (2)], as explicitly shown in the Appendix A. The dislocation defect, which is created by the Volterra construction, with the associated hopping π phase factor, gives rise to the level repulsion among four zero modes at the pasted edges across the Volterra cut. However, a Kramers pair of modes |Ψ 0 still remains pinned at zero energy and gets localized in the defect core [11].
The crucial observation is that the HOT mass matrix Γ 4 commutes with the dislocation or edge-mode matrix A b , [Γ 4 , A b ] = 0. The HOT mass matrix therefore reduces in the eigen-subspaces of A b , introducing the level repulsion between the two zero modes and thus symmetrically splits them about the zero energy. Furthermore, this implies that the modes do not change the form, i.e. they remain localized around at the defect, after introducing the HOT mass. The energy splitting reads as Unless ∆[K inv − ib(b · ∇)] = 0, i.e. when the HOT mass vanishes in the direction of the Burgers vector, the energy splitting is non-zero (δE = 0). Therefore, the interplay between the orientation of the HOT mass domain wall and the Burgers vector pins the dislocation modes precisely at zero energy when the two directions are parallel [see Fig. 1(a)]. The mechanism for the splitting of the dislocation modes is also operative for gapping out the edges, ultimately yielding the corner modes, implying that a dislocation can directly probe higher-order bulkboundary correspondence. This mechanism captures the existence of the localized pair of dislocation modes in the M phase, split by the energy gap whereẼ

and Ψ
(1,2) 0 (x) are the zero-energy dislocation states in the first-order phase (see Appendix A). The obtained energy gap (δE) implies that the modes are pinned at zero energy only when the Burgers vector is parallel with the mass domain wall (θ = π/2), see Fig. 2(a). Furthermore, for small ∆, δE scales linearly with cos θ [see Fig. 2 Fig. 2(c)]. Finally, a dislocation does not feature any bound states either when Φ dis = 0 (as in the Γ phase) or in the trivial phase.

V. DISLOCATIONS IN 3D HOT INSULATOR
The above arguments can be straightforwardly extended to 3D second-order topological insulators. Notice first that a 3D edge dislocation can be obtained from its 2D analogue by translating it along an out-of-plane lattice vector [40]. Therefore, an edge dislocation in a 3D second-order topological insulator should in general feature gapped modes, which, however, become gapless for the Burgers vector parallel to the HOT mass domain wall direction, analogously to the 2D case [see Fig. 5

(a)].
A screw dislocation, being a true 3D defect, features the Burgers vector b parallel to its orientation. An electron encircling the dislocation defect once, skips a lattice distance |b| along the defect relative to the perfect crystal [40]. The "K · b" rule implies that an electron then picks up a phase Φ dis = K inv ·b. When this phase is nontrivial in a translationally-active first-order insulator, the screw dislocation hosts gapless propagating modes [11]. On the other hand, in a translationally-active HOT insulator some of the surfaces are gapped, and the screw dislocation hosts gapless propagating modes only when it is oriented perpendicular to gapless surfaces. See Appendix B. In particular, in a 3D second-order topological insulator, with a single mass domain wall [see Eq. (1)], a screw dislocation parallel to it hosts gapless modes, since the defect then pierces gapless xy surfaces. Otherwise, a dislocation perpendicular to gapped xz or yz surfaces features gapped one-dimensional modes. See Fig. 5 We numerically confirm this scenario in a 3D secondorder topological insulator described by the tight-binding model, exemplifying the R−phase on the cubic lattice (see Eq. (1) and Appendix D). Both single edge dislocation with Burgers vectors b = ae x extending in the z−direction in an open system, which is also the propagation direction of the zero-energy hinge modes, and a edge dislocation-antidislocation pair extending in the same direction with Burgers vectors b = ±ae x in a periodic system indeed yield finite energy states when θ = π/2, as shown in Fig. 3(a) and Fig. 3(b), respectively. For θ = π/2 the dislocation modes become gapless (same as in 2D). See Figs. 5(a),(b),(c). On the other hand, for a single screw dislocation with b = ae z in a open system (coexisting with the zero-energy hinge modes) and a screw dislocation-antidislocation pair with b = ±ae z in a periodic system, we obtain gapless dislocation modes for any HOT mass domain wall orientation θ in Eq. (1), as displayed in Figs. 3(c) and 3(d). See also Fig. 5(f) where vanishing of the gap is explicitly shown. The propagating modes are gapless in this case because the defect pierces the gapless xy surfaces, for any θ. The dislocation modes are localized within a few lattice sites at the defect core, as can be seen from their LDoS in a plane perpendicular to the dislocation direction displayed in Figs. 5(d) and 5(e). The screw dislocation modes also inherit the C 4 symmetry preserved by the defect.

VI. DISCUSSION AND CONCLUSION
Our findings are experimentally consequential for probing HOT insulators, apart from the crystalline systems, also in metamaterials.
The paradigmatic model of the 2D second-order topological insulator, the Benalcazar-Bernevig-Hughes (BBH) model [23], equivalent to the minimal lattice model in Eq. (1) [41], has been realized in the lattice of microwave resonators [30]. A dislocation defect in this setup should be created by a local hopping modification through π phase factors across a line of missing sites ending at the dislocation center, analogously to the case of a translationally-active first-order topological insulator [42], and a disclination [43]. In the BBH photonic lattice, where the sign of the hopping also can be locally manipulated [33], it should be therefore possible to introduce the dislocation defects and observe the defect modes, as in first-order 2D topological photonic crystals [44,45], and for a disclination defect [46]. Finally, the artificial lattices can host HOT phases, as recently shown for Kagome lattice [35], and we expect that because of their tunability, our theoretical predictions can also be directly tested in these platforms.
Most of the proposed and experimentally studied 3D HOT crystalline materials turn out to be of the translationally-active type. For example, elemental Bi exhibits a double band inversion at the C 3 -symmetric Tpoint in the BZ and supports mixed electronic topology manifesting through coexisting gapless hinge and Dirac surface modes [21,28,47]. Our general mechanism thus implies that an edge dislocation with the Burgers vector in the (111) direction parallel to the HOT mass domain wall (see, Fig. 1c in Ref. [28]), so that also Φ dis = π, features gapless modes, protected by C 3 and time-reversal symmetries. A screw dislocation oriented in the same direction should host one-dimensional gapless (gapped) states if the (111) surface is gapless (gapped). Similarly, for a recently proposed HOT insulator in Zr(TiH 2 ) 2 , with band-inversion away from the Γ point and the gapless modes along all the edges in the cubic geometry [48], we predict gapless (gapped) modes in the core of an edge (a screw) dislocation with the Burgers vector along a principal crystal axis. Finally, the candidate HOT insulators Bi 4 X 4 , with X=Br,I, [48][49][50][51] feature band inversions at R and M = (π/a, π/a, 0) points in the BZ, and hence the dislocations should host the one-dimensional modes, following the above general rule.
Here we demonstrated that dislocations can be instrumental in probing higher-order electronic topology in insulators as a consequence of the subtle interplay between the geometry of the HOT Wilson-Dirac mass in the momentum space and real-space lattice topological defects. We furthermore demonstrate the protection of the dislocation modes in the case of C 2n rotational symmetry breaking HOT insulators in both two and three dimensions (see Appendices E and F), which pertains to physically relevant C 4 and C 6 symmetric crystals. On the other hand, when the order of the rotation is odd, i.e. for C 2n+1 rotations, there is no higher-order Wilson-Dirac mass term that changes sign under such rotation. Therefore, we cannot find second-order topological mass. Further analysis of this case is left for future investigation. Recently, it has been shown that also partial dislocations with a Burgers vector which is a fraction of a primitive lattice vector can host gapless propagating modes in 3D HOT insulators [52]. In addition, disclina-tion can also host topological modes in 3D HOT insulators [53]. Our findings motivate future investigation of the response to the dislocations in HOT insulators on different crystalline lattices, such as, for instance, Kagome lattice [54]. Finally, we expect that our mechanism will be a useful guide for the experimental detection of the HOT phases in diverse platforms, and consequential also for HOT semimetals [36,55] and superconductors.

Data and code availability
The data that support the plots within this paper and other findings of this study are available from the authors upon reasonable request.
Here, Ψ k is a four-component spinor, the exact form of which does not affect the following discussion, while k is the momentum, and a is the lattice spacing. The mutually anticommuting four-component Γ matrices satisfy the Clifford algebra {Γ j , Γ k } = 2δ jk for j, k = 1, · · · , 5. The following discussion only rests on this anicommuting Clifford algebra, not on the exact representations of the Γ matrices.
This model for the regime of parameters 0 < m/B < 8 describes a 2D first-order topological insulator and a HOT insulator (second-order) for ∆ = 0 and finite ∆, respectively. Furthermore, when 4 < m/B < 8, the model features the band inversion at the M = (π/a, π/a) point (the M phase), while for 0 < m/B < 4, the band inversion is at the Γ = (0, 0) point (the Γ phase) of the BZ. Notice that {ĥ 0 ,ĥ ∆ } = 0, and thereforeĥ ∆ acts as a mass term for the topological edge states ofĥ 0 . This mass term changes sign under the C 4 rotation and, as such, assumes the profile of a discrete symmetry breaking Wilson-Dirac mass, the exact form of which depends on the parameter θ ∈ [0, π/2]. In particular, for θ = π/2 the domain wall lies along the diagonals k y = ±k x , while for θ = 0 it is located along the principal axes, k x = 0, k y = 0.
We now comment on the structure of the corner modes for θ = 0 and π/2. Let us consider a 2D square lattice of linear dimension L in each direction, such that four corners are at (±L/2, ±L/2). Four corner modes are then sharply localized around these corners for θ = 0, and with increasing θ they become more delocalized. By contrast, if the crystal is cut in such a way that four corners are located at (±L/2, 0) and (0, ±L/2), the corner modes are most prominently localized when θ = π/2 and gradually delocalize as θ is ramped down to zero. However, irrespective of the sharpness of the corner modes, the system always describes a 2D second-order topological insulator for any θ. These outcomes are shown in Fig. 4. A similar structure also appears for the hinge modes for 3D second-order topological insulator, discussed in Appendix B.
The continuum Hamiltonian is obtained by expanding the above lattice Hamiltonian close to the bandgap closing at the M point, k x = π + q x , k y = π + q y , with ∆ = 0, which in the real space (q → −i∇) reads aŝ Here,m = 8B − m > 0,B = B > 0, and we set a = 1.
The two edges along the lines x ± = ±a [ Fig. 1(a)], before the dislocation is introduced through the Volterra construction, feature topological gapless edge states, resulting from the corresponding zero modes of the edge where the vector of Pauli matrices µ acts in the space of the two edges. When dislocation introduces a π hopping phase, the reconnection of the edges across the Volterra cut is modeled by a hopping Hamiltonian between them in the form H D = t sgn(x)Γ 1 ⊗ µ 1 , where sgn(x) is the "sign" function. This term, through the sign change of the hopping across the cut, takes into account that the (low-energy) electrons acquire a π phase when encircling the defect. Its form ensures that when the phase factor is trivial, the connected edge modes are trivially gapped out, as in the Γ phase.
We then look for the zero energy modes of the Hamiltonian H edge + H D , or explicitly After multiplying this equation from the left by −iΓ 1 ⊗ µ 3 , we obtain The form of the above equation implies that we can take the following ansatz for the solution where which for an exponentially localized solution at both edges f (x) ∼ exp(−λ|x|), yields − tλ + (m +Bλ 2 )σsgn(x) + tρ = 0.
Choosing σsgn(x) = +1, ρ = +1, and considering the regime t 2 > 4B(m + t), we obtain two characteristic inverse localization lengths Notice that the spinors χ σ are doubly degenerate, because of the anticommuting property of the four-component Hermitian Γ matrices. This form of the localization length together with the continuity condition across the reconnected edges, Ψ 0 (x → 0) = 0, yields a pair of the zero energy dislocation modes given by mentioned in the main text.
The rest of the analysis straightforwardly follows from the fact that [Γ 4 , iΓ 1 Γ 3 ] = 0 implying that Γ 4 reduces in the eigen-subspaces of iΓ 1 Γ 3 . Therefore, a perturbation proportional to Γ 4 symmetrically splits the dislocation modes about the zero energy. Explicitly, the HOT perturbation after the expansion about the M point at K M = (π/a, π/a), k x = π/a + q x , k y = π/a, with q x = −i∂ x , takes the form Therefore the splitting of the modes is given by as announced in the main text. Finally, we note that in the above derivation of the dislocation modes in d = 2 and also in d = 3 (see Appendix B), we only used the Clifford algebra of the Γ matrices, {Γ i , Γ j } = 2δ ij , for i, j = 1, ..., 5, and the existence of the modes is therefore independent of the Γ−matrix representation. Since [Γ 4 , Γ 13 ] = 0, even though introduction of the HOT mass generically yields finite energy dislocation modes unless θ = π/2, they cannot be mixed with rest of the bulk states, and therefore remain robust and protected. We further substantiate this claim by demonstrating the symmetry protection of the dislocation modes in Appendix E in the case of C 4n rotational symmetry, while an analogous analysis for C 4n+2 rotational symmetry is presented in Appendix F.

Appendix B: 3D edge and screw dislocation modes: The continuum model
In this Appendix we start with the lattice model for a three-dimensional second-order topological insulator on a cubic lattice, analogous to Eq. (A1), but with [Eqs. (1) and (2) We now consider the R phase, with the band inversion at the R = (π/a, π/a, π/a) point in the BZ, realized for the values of the parameters 8 < m < 12, t = 1 and B = 1, and expand the above Hamiltonian about the R point for ∆ = 0 to obtain wherem = 12B − m > 0 in the topological R phase. We take B = 1, and k i = π/a + q i , for i = x, y, z.
A 3D edge dislocation is obtained by translating its 2D counterpart along a particular lattice direction representing the defect line. Therefore, the conclusions obtained above in the 2D case for a Burgers vector along the C 4 symmetry breaking x or y direction directly apply to the 3D case: an edge dislocation should in general feature gapped modes which become gapless when the Burgers vector is parallel to the domain wall (θ = π/2), consistent with our numerical findings, see Fig. 5(a)-(c). Also, when the Burgers vector is oriented along the C 4 symmetry axis, therefore piercing (f) Scaling of the spectral gap (δE) among the states localized at the dislocation core in a system of linear dimension L = 16 in each direction, in the presence of a screw dislocation-antidislocation pair, when the Burgers vector is b = ±aez (red) and b = ±aey (blue). Respectively in these two cases the screw dislocation pierces surfaces hosting gapless and gapped modes, and concomitantly the dislocation modes are also gapless and gapped. a gapless surface (in this case the xy surface), the defect should feature gapless modes.
A screw dislocation with the Burgers vector b = ae z can be introduced by the Volterra construction as follows. We first choose a slip half-plane for y = 0, x > 0, whose neighbor half-plane y = ae y , x > 0 is displaced by the Burgers vector b = ae z relative to the slip plane. Translational symmetry is then restored by reconnecting the bonds between the slip and the neighbor plane everywhere, except close to the dislocation line along the zdirection. An electron sliding down through the dislocation picks up a phase Φ dis = K inv ·b = π upon encircling it once.
The Hamiltonian for the slip plane and its neighbor before the reconnection reads where the vector of Pauli matrices µ acts in the space of the two surfaces. The π phase factor is then introduced by modifying the sign of the hopping in the y-direction across the domain wall at x = 0, because the slip plane is orthogonal to the y-axis, yielding The form of this Hamiltonian is chosen so that when the phase factor is trivial, the surface zero modes are gapped out, as in the Γ phase, realized for t = B = 1 and 0 < m < 4 in Eq. (B1). We now look for the zero modes of the surface Hamiltonian using that the translational symmetry is preserved in the z-direction and is explicitly broken in the x-direction. Hence, k x = π/a − i∂ x , k z = π/a, and which is identical to the Eq. (A4), up to the change Γ 3 → Γ 4 . The form of the pair of zero modes is therefore given by Eq. (A9), with the difference that here iΓ 1 Γ 4 χ σ = σχ σ . Gapless propagating modes along the dislocation line are obtained by "translating" the zero modes along the dislocation direction Ψ Notice that this solution explicitly breaks C 4 symmetry which is a consequence of the specific Volterra cut not preserving this symmetry. Namely, there are two equivalent slip planes, x − z and y − z, each of them individually breaking C 4 symmetry. On the lattice, this symmetry is certainly preserved away from the dislocation line and therefore the zero mode solution in the continuum inherits it, so that, after a proper regularization the spatially dependent part should obey Ψ 3D (x, y) = Ψ 3D (−y, x). This feature, as we show below is inherited by the HOT dislocation modes, and as we also find in our numerical analysis, see Fig. 5(e) for the plot of the local density of states (LDoS) of the closest to zero energy modes on a particular x − y plane.
We now introduce the HOT mass term proportional to Γ 5 , as given by the Hamiltonian h ∆ in Eq. (B1), that may gap out otherwise gapless propagating modes along the dislocation since [Γ 5 , iΓ 1 Γ 4 ] = 0, analogous as in the case in two dimensions. The splitting depends also on the form factor ∆(k, θ), which after expanding about the band inversion R point becomes which then yields for the energy gap because of the C 4 symmetry [both (∂ 2 x − ∂ 2 y ) and ∂ x ∂ y are odd, while Ψ(x, y) is even under the C 4 ]. The x − y surface therefore yields gapless modes for the screw dislocation piercing it, as also expected from the fact that it is gapless. This is precisely what we find in the numerical analysis, see Fig. 5(f) (red dots). The plot of the LDoS of the closest to zero energy modes shows explicitly the C 4 symmetry, see Fig. 5(e). The result can also be appreciated from the fact that in the continuum limit, featuring the full U (1) rotational symmetry in the x − y plane, the dislocation represents a π magnetic flux tube, which is known to carry the gapless mode when piercing a gapless surface [11,17]. On the other hand, when the screw dislocation is oriented in the x or y direction, piercing therefore a gapped surface, it carries gapped modes. This can be directly seen by recalling that the dislocation zero modes have to break the C 4 symmetry, and thus the spatial part of the zero energy solution for the dislocation with Burgers vector in the y direction is Ψ(x), as given in Eq. (A9). The energy splitting is then for any θ = π/2. This is consistent with the numerically extracted scaling of δE versus the HOT mass parameter ∆, shown in Fig. 5(f) (blue dots). Therefore, as long as the screw dislocation pierces the gapped surfaces, the modes along the dislocation line are also gapped with the gap scaling with the size of the surface gap (∆). The dislocation modes, however, remain within the bulk bandgap, determined by the first-order topological mass proportional to the Γ 4 matrix in Eq. (B1).

Appendix C: Unitary and anti-unitary particle-hole symmetry
In this Appendix we discuss the unitary and anti-unitary particle-hole symmetry of both first-order and secondorder topological insulators that respectively protects zero-energy edge, surface, corner and hinge modes, and both zero and finite energy dislocation modes. For this analysis we choose the following representation of the Γ matrices satisfying the anticommuting Clifford algebra      Fig. 1(a) of the main text). The bulk states are shown in black. As the HOT mass domain wall becomes parallel to the Burgers vector (θ = π/2), the dislocation modes become zero-energy states. Note that for any value of θ the dislocation modes (at zero energy or finite energies) are always well separated from the bulk states, and they are particle-hole partners of each other, due to the unitary particle-hole symmetry generated by Θ = Γ5 (see Appendix C), and thus do not mix with the bulk states. Hence, they are protected by the bulk topology and cannot be removed from the system (by mixing with the bulk states, for example). Here, n is the index for the energy eigenvalues (En). In Appendix E, we also show that the composite symmetries of HOT insulators forbid addition of any new perturbation in the system, leaving the dislocation modes at finite and zero energy symmetry protected.
in the Hamiltonian for the second-order HOT insulator in d = 2 and d = 3, given by Eq. (1) in the main text.
Next we discuss the particle-hole or the spectral symmetry of both the first-order and second-order topological insulators in d = 2 and d = 3. It is defined in terms of an operator, say Θ, which anticommutes with a generic Hamitlonianĥ gen , i.e., {ĥ gen , Θ} = 0. If such operator (Θ) exists, then all the eigenstates of theĥ gen at positive and negative energies ±E n , denoted by | ± E n , are related to each other according to Θ| ± E n = | ∓ E n . This statement follows from the relation h gen | ± E n = ±E n | ± E n ⇒ Θ ĥ gen | ± E n = ±E n (Θ| ± E n ) ⇒ −ĥ gen (Θ| ± E n ) = ±E n (Θ| ± E n ) and hence Θ| ± E n = | ∓ E n . Otherwise, the particle-hole symmetry generator Θ can be either unitary or antiunitary [41]. The above definition of the particle-hole or spectral symmetry then shows that if there exist  (1) in the main text. The representation of the Γ matrices is given in Eq. (C1). Time-reversal (T ), parity (P) and C4 rotational (R4) symmetries are, respectively, represented by the operators T = τ2 ⊗ σ1K, P = τ1 ⊗ σ3, and R4 = exp i π 4 τ0 ⊗ σ3 . For a scalar (pseudoscalar) M under an operation X it holds XM X † = X (XM X † = −X). Scalar and pseudoscalar operators are denoted by and × in the table. Here (×) also indicates whether an operator is even (odd) under a specific symmetry.
any state at precise zero energy, which can only be achieved in the true thermodynamic limit, then it must be an eigenstate of Θ with eigenvalue +1 or −1. However, in any finite system, all the states are at finite energies, and the states which reside at zero energy in the thermodynamic limit are also placed at finite but extremely small (typically ∼ 10 −6 − 10 −8 ) energies. Next we show that such particle-hole or spectral symmetry operator always exists for the universal models of both first-order and second-order topological insulators in d = 2 and d = 3 [see Eq. (1)], and thereby protects both zero energy topological boundary modes (such as the edge, surface, corner and hinge modes), as well as the dislocation modes (both zero and finite energy ones).
We begin with the first-order topological insulators. Notice that for such systems in both d = 2 and d = 3, we can always find a unitary particle-hole symmetry operator, namely Θ = Γ 5 , irrespective of the representation of the Γ matrices. We note that with the representation of the Γ matrices, shown in Eq. (C1), we also also find one antiunitary particle-hole symmetry generator, namely Θ = (τ 2 σ 2 )K, for 3D first-order topological insulator, where K is the complex conjugation and Θ 2 = +1. Such antiunitary particle-hole operator will play an important role for 3D HOT insulator. For now, one can immediately appreciate that the unitary particle-hole operator Γ 5 , pins the boundary edge and surface modes at zero-energy and prevents them from mixing with the bulk states. In addition, it also forbids the gapless dislocation modes in both 2D and 3D first-order topological insulators from mixing with the bulk mode. Next, we proceed to HOT insulators.
For a 2D HOT insulator, we can still find Θ = Γ 5 . It protects the four zero-energy corner states [see Fig. 1(b)]. Furthermore, given that the in-gap, but finite-or zero-energy dislocation modes are always well separated from the bulk states, see Fig. 6, and appear in pairs, such unitary particle-hole operator also prevents the finite-energy (for any θ = π/2) and the zero-energy (for θ = π/2) dislocation modes from mixing with the bulk states. Thus neither the zero nor the finite energy dislocation modes can be removed from the system, and they are protected by the particle-hole symmetry. Finally, for a 3D HOT insulator, we cannot find any unitary particle-hole operator, since the corresponding Hamiltonian operator exhausts all five mutually anticommuting four-component Γ matrices. Nevertheless, the particle-hole symmetry is then generated by an antiunitary operator Θ = (τ 2 σ 2 )K. Then the above discussion for the 2D HOT insulator immediately generalizes to 3D HOT insulator, and one can immediately conclude that (1) the four zero-energy hinge modes, and (2) both finite-and zero-energy dislocation modes are protected by the antiuntary particle-hole symmetry. Finally, we note that the existence of such an antiunitary particle-hole symmetry generator is also independent of the explicit representation of the Γ matrices [41]. In Appendix E, we show that composite symmetries of both 2D and 3D HOT insulators forbid addition of any new term to their corresponding universal Hamiltonian, shown in Eqs. (A1) and (B1), respectively. Therefore, both finite and zero energy dislocation modes are symmetry protected.

Appendix D: Details of the numerical analysis and additional results
This Appendix is devoted to highlight the parameter values used in Figs. 1 and 2 of the main text, and discuss additional numerical results for three-dimensional second-order topological insulators, hosting one-dimensional hinge modes along the z direction, in the presence of edge and screw dislocation-antidislocation pair, shown in Fig. 5.
For two-dimensional topological insulators, we always choose t = 1, B = 1/2, and m = 3 [see Eq. (A1)], so that the system is in the M phase. To realize a second-order topological insulator, we set ∆ = 0.20 in Figs. 1(b) and 1(c). For Figs. 2(a)-(c), we consider a dislocation-antidislocation pair in a periodic system with linear dimension L = 28 in each direction. The locations of the dislocations in these figures are identical to the ones in Fig. 1(c). For three-dimensional topological insulators, we always set t = B = 1 and m = 10 [see Eq. (B1)], so that the system is in the R phase. To realize a second-order topological insulator, we set ∆ = 0.40 in Fig. 3.
We now discuss additional numerical results for 3D second-order topological insulators in the presence of dislocations, displayed in Fig. 5. In Figs. 5(a)-(c) we show the scaling of the spectral gap (δE) among four states localized in the dislocation core in the presence of an edge dislocation-antidislocation pair with Burgers vectors b = ±ae x . It shows that only when the Burgers vector is parallel to one of the C 4 symmetry breaking axis (θ = π/2) these modes become gapless [ Fig. 5(a)]. Otherwise, the spectral gap scales linearly with cos θ and ∆, as shown in Figs. 5(b) and (c), respectively, when ∆ is small, in agreement with our scaling argument presented in Appendix B.
Figs. 5(d) and (e) show the localization (through LDoS) of the dislocation modes around the core of edge and screw dislocation-antidislocation pair, respectively, in the xy plane for a specific z = 6. These two figures are projections of Figs. 3(b) and (d) of the main text, respectively, for a specific z.

Fig
. 5(f) shows the scaling of the spectral gap (δE) between the states localized in the dislocation core in the presence of a screw dislocation-antidislocation pair for two specific orientations of the corresponding Burgers vectors, namely, when b = ±ae z (red) and b = ±ae y (blue). For b = ±ae z the screw dislocation pair pierces the xy surfaces that host gapless modes, while for b = ±ae y it pierces the xz planes which, on the other hand, accommodate gapped modes. Concomitantly, for these two orientations of the screw dislocation-antidislocation pair the modes localized in its core are respectively gapless and gapped, in agreement with our analytical arguments from Appendix B.
Finally, we show the real space version of the tight-binding model, reported in Eq. (1) of the main text in the momentum space, in the presence of a single edge dislocation in open system on a 2D square lattice for θ = 0. By performing the Fourier transform, we arrive at In this construction, the dislocation center is located at ( x , y ) [see Fig. 1 of the main text]. Now this construction can be generalized to introduce a pair of edge dislocation-antidislocation in a 2D periodic system, as well as edge and screw dislocations in 3D.

Appendix E: Symmetry protection of dislocation modes
In this Appendix we take the representation of the Γ matrices as given by Eq. (C1) in the Hamiltonian for the second-order HOT insulator in d = 2 and d = 3, given by Eq. (1) in the main text. Antiunitary time-reversal (T ), unitary parity (P) and C 4 rotational (R 4 ) symmetries are respectively represented by the operators T = τ 2 ⊗ σ 1 K, P = τ 1 ⊗ σ 3 , and R 4 = exp i π 4 τ 0 ⊗ σ 3 , where K is the complex conjugation. Under time-reversal and parity, momentum changes sign, k → −k, while it is a vector under the C 4 rotation: (k x , k y ) → (−k y , k x ). The first-order topological Hamiltonian (ĥ 0 ) and the second-order mass term (ĥ ∆ ) are both invariant only under the composite PT , C 4 T , C 4 P symmetries, see Table I. We show that there are no additional symmetry allowed terms that can be added to the HOT Hamiltonian besides the ones already present and the trivial one, proportional to unity matrix, which is, of course, always symmetry allowed in an insulating system, and it only causes an overall shift of energy eigenvalues. Finally, we point out that this proof pertains to any second-order topological insulator with C 4n rotational symmetry breaking second-order Wilson-Dirac mass because in this case for any n the HOT mass is even under k → −k.
As a first step, we classify all the sixteen four-dimensional Hermitian matrices, constructed from five mutually anticommuting Hermitian Γ matrices given by Eq. (C1), according to the above three composite symmetries, see Table II. We see that only the matrix Γ 3 , which represents the first-order topological mass, is invariant under all the three composite or product symmetry operations. One can also form scalars under the C 4 rotation in terms of the four vectors in the form, where j = 1, 2, 3, 4 corresponds to four vectors Y j from the last four rows in Table II, and the angle χ → χ + π/2 under the C 4 . However, as explicitly shown in Table III, none of these Hamiltonians is invariant under both parity and time-reversal, and therefore they violate either C 4 T , C 4 P, or PT . Notice that the PT symmetry assures the double degeneracy of the electronic bands in the model, since it is antiunitary and squares to −1 [41]. Therefore, there are no new symmetry allowed terms in the 2D HOT Hamiltonian given by Eq. (1) in the main text besides the ones already present. This in turn implies that both both finite and zero energy dislocation modes in a 2D second-order topological insulator are symmetry protected by the same combinations of the spatial and non-spatial symmetries as the bulk Hamiltonian.
In three spatial dimensions, we choose the first order Hamiltonian so that the matrix Γ 5 multiplies the factor t sin(k z a) (see Table II), while we keep the rest of the terms identical as in the Hamiltonian for the 2D second-order topological insulator [a slightly different notation than the one in Eq. (1)]. Therefore, requiring the same symmetries, the above results imply that the same composite symmetries protect the dislocation modes in 3D as in two spatial dimensions.
In Fig. 6, we show that in the presence of dislocations the energy spectra possess two gaps. One is set by the firstorder topological mass (among the black colored bulk states) and the other one is determined by the HOT mass ∆(k, θ) among the dislocation modes (red states). The latter one in addition also depends on the parameter θ, measuring the relative orientation between the HOT mass domain walls and the Burgers vector. Only when θ = π/2, for which the HOT mass domain wall is parallel to the Burgens vector, the dislocation modes become gapless. However, as we argued above, these two distinct energy scales are symmetry protected and we cannot add any symmetry allowed perturbation that can mix these two scales. The same conclusion holds analogously in 3D. Thus both finite and zero energy dislocation modes are symmetry protected.
Appendix F: Symmetry protection for C4n+2 symmetry breaking second-order topological insulator In this Appendix, we discuss the symmetry protection of the dislocation modes when the rotational symmetry is of the form C 4n+2 , where n is an integer, as for instance is the case for C 6 rotations. Specifically for C 6 symmetry Operator Symmetry in the M phase [compare with Eq. (1)] [56]. In this case, the second-order Wilson-Dirac mass is odd under k → −k, and the protection of the dislocation modes can be realized only through the parity (P) and the time-reversal (T ) symmetries (no rotational or a composite symmetry is required now). Therefore, the inspection of the Table II, shows that the only Γ matrices preserving both parity and time-reversal symmetries are Γ 0 and Γ 3 . Then our argument from Appendix E can be extended straightforwardly to show that dislocation modes (at both finite and zero energies) cannot be mixed with the other bulk states and hence they are symmetry protected. Furthermore, this argument straightforwardly generalizes to all 2D second-order topological insulators for which the second order Dirac-Wilson mass is odd under the transformation k → −k. Finally, due to the form of the Hamiltonian for 3D second-order topological insulators, this argument straightforwardly applies also to this case.