Second-order topological modes in two-dimensional continuous media

We present a symmetry-based scheme to create 0D second-order topological modes in continuous 2D systems. We show that a metamaterial with a \textit{p6m}-symmetric pattern exhibits two Dirac cones, which can be gapped in two distinct ways by deforming the pattern. Combining the deformations in a single system then emulates the 2D Jackiw-Rossi model of a topological vortex, where 0D in-gap bound modes are guaranteed to exist. We exemplify our approach with simple hexagonal, Kagome and honeycomb lattices. We furthermore formulate a quantitative method to extract the topological properties from finite-element simulations, which facilitates further optimization of the bound mode characteristics. Our scheme enables the realization of second-order topology in a wide range of experimental systems.

The development of topological insulators (TIs), while originating in electronic systems, has made a profound impact in the field of classical metamaterials. Their signature topological boundary modes have been successfully demonstrated in photonic [1,2], electrical [3,4], phononic [5][6][7], acoustic [8,9], atomic [10] and polaritonic [11] systems. This remarkable universality stems from the origin of the boundary phenomena; these depend only on topological invariants derived from the bulk spectral bands, and not on the specific medium. Early realizations of topology in classical systems relied on engineering synthetic gauge fields in tight-binding (TB) models, obtained by coupling resonator modes in space [12,13] and/or time [14] to create topological boundary modes. Later it was discovered that in patterned continuous media, topologically distinct phases can be formed purely by breaking the spatial symmetries of the pattern [15]; this has opened the field of TIs to a wider range of platforms, and brought it closer to potential applications. Furthermore, the topological phases of matter paradigm has been recently extended to higherorder TIs (HOTIs) [16][17][18][19][20][21][22][23][24]. In a d-dimensional HOTI, topological modes of dimension d − 2 and lower may appear. HOTIs have been experimentally demonstrated by emulating TB models in a variety of platforms [25][26][27][28][29][30][31][32][33][34].
The concepts of TIs and HOTIs are unified in the framework of ten topological classes, each characterized by time reversal, chiral and particle-hole symmetries [35]. For a d-dimensional system with a defect defined on a Ddimensional surface, a topological invariant exists whose possible values depend only on the bulk symmetry class and d − D [36]. At such a defect, modes with dimension d − 1 − D (codimension D + 1) are formed. TIs, with their topology defined purely in their bulk, correspond to D = 0; for D > 0, the invariant generally includes integrals over curvatures combining real and momentum space, leading to the appearance of high-order topological modes [16-20, 22-24, 36] with topological boundary states of codimension > 1. An instance of a HOTI is found in the d = 2 Jackiw-Rossi Hamiltonian of a topological vortex, which can be implemented by deforming the TB model of graphene [37]. Emulating this model has successfully demonstrated the creation of 0D modes at topological vortices in photonic [38][39][40] and elastic [41] devices. Compared with standard 0D modes formed at non-topological lattice defects, topological vortex modes offer several advantages [39]. (i) Their frequency is near mid-gap, resulting in better spatial confinement and quality factors. (ii) The number of modes formed is fixed by topology. (iii) The modal area is scalable, as the modes form independently of the vortex size. These aspects make topological vortex modes promising candidates for the construction of single-mode semiconductor lasers.
In this Letter, we show how to obtain the Jackiw-Rossi Hamiltonian in continuous 2D structures. Crucially, our approach does not require a recourse to a TB model or threading by magnetic fluxes. Instead, we utilize solely the breaking of spatial symmetries of a linear medium. We focus on the p6m space group, pick-  Fig. 1(III)(a)] lattices. Central to our work, the group p6m has a 4D irreducible representation that manifests as two Dirac cones at the K/K points in reciprocal space. We consider two distinct perturbations of the primitive cell, namely the breaking of inversion and translation symmetries, and find their matrix representations in the 4D eigenspace [42,43]. The perturbation matrices are shown to anticommute and thus correspond to different gap-opening terms; these constitute the Jackiw-Rossi model. By spatially varying the perturbations in a single system, a topological vortex defect is formed with in-gap 0D modes bound within. For concreteness, we assume a dielectric implementation of the scheme; we stress however that our analysis holds for any patterned 2D linear medium and is thus applicable to a wide range of systems.
The 2D Jackiw-Rossi Hamiltonian reads where k = (k x , k y ) are quasimomenta, and m = arXiv:2103.05403v2 [cond-mat.other] 31 Jul 2021 (m 1 , m 2 ) are mass terms that can change in space r = (x, y). The four 4 × 4 matrices, γ = (γ 1 , γ 2 ) and Γ = (Γ 1 , Γ 2 ), anticommute with one another. If either of the mass terms m 1 or m 2 is nonzero, the spectrum of H is gapped around zero energy. We first consider the case where one mass term (e.g., m 2 ) vanishes and the other switches signs in space [e.g., m 1 = |m| sgn(x)], thus defining a domain wall at x = 0. The distinct band topology of the two bulk phases then defines a topological defect (D = 0) and a band inversion occurs across x = 0. The domain wall can then be described by two copies of the Jackiw-Rebbi model [15,44], and the system admits a Z 2 invariant [45,46]. Correspondingly, two counter-propagating 1D edge states appear at the domain wall. Depending on the eigenbasis of the gapping term, we thus obtain the states associated with either the quantum spin Hall effect (QSHE) or the quantum valley Hall effect (QVHE). Allowing instead both masses to vary over a spacedependent closed path defines a defect with D = 1 [36]. The Hamiltonian is then characterized by a Z topological invariant associated with the winding number n of m(r). Taking m 1 +im 2 = |m|e iθ , the winding number is defined as where S 1 denotes the closed path. At such a topological vortex, n zero-dimensional bound modes are found. Systems with |n| > 0 can be constructed within the TB model of graphene, where m 1 and m 2 are generated by spatial perturbations of the six-site unit cell [37]. We also emphasize that this model is topologically equivalent to other proposed TB models for HOTIs [23] and the invariant definition involves integrating over both real and momentum space, such that the winding encircles a quantized 4D Dirac cone [24]. Symmetry analysis.-We focus on the p6m space group, exemplified by the three unperturbed structures shown in Fig. 1(a). All three are symmetric under sixfold rotation C 6 , reflection σ x , and lattice translation T . The operations that leave the structures invariant form a group G that is isomorphic to p6m. The group G admits a 4D irreducible representation (irrep) derived from the little group of the K/K points [47]. This irrep manifests itself as the familiar pair of Dirac cones, appearing at the valleys K and K in reciprocal space. Choosing a basis for this irrep defines a function ρ(g), which maps each symmetry element g ∈ G to a 4 × 4 representation matrix ρ(g). For the three group generators, we obtain [48] where τ i are Pauli matrices, R 6 denotes a 2 × 2 rotation matrix of angle π/3, and w = exp{2πi/3}. The basis is chosen such that each ρ(g) is a direct product of a 2 × 2 matrix acting on the two valley degrees of freedom with another 2 × 2 matrix acting on the two degenerate bands within each valley. Opening band gaps in the degenerate cones to create topologically distinct bulks can be achieved in two different ways: breaking inversion or translation symmetry, as shown in Fig. 1(b) and 1(c), respectively. Starting with the former, we break inversion symmetry while preserving C 3 ; this entails different deformations in each of the three example structures, see Fig. 1(b). These deformations leave the symmetry subgroup generated by {C 3 , σ x , T } intact; the resulting space group is isomorphic to p3m, and the corresponding symmetry-breaking perturbation takes the form [48] where Γ 1 = τ 3 ⊗ τ 2 , and c 1 and m 1 are real coefficients that quantify the strength of the perturbation. Breaking translation symmetry, instead, we consider an enlarged real-space unit cell [ Fig. 1(c)], which leads to band folding into a smaller Brillouin zone (BZ) in reciprocal space. The larger unit cell maps both valleys K/K to the new Γ point. Applying the appropriate deformation for each lattice then breaks the original translation invariance so that the system is no longer self-coincident under T . The remaining symmetry subgroup is generated by {C 6 , σ x , T C 6 T }, where T C 6 T is the expanded lattice translation. The corresponding translation symmetrybreaking perturbation matrix is [48] with Γ 2 = τ 1 ⊗ 1 2 , and real coefficients c 2 , m 2 . The perturbation matrices Γ 1 and Γ 2 anticommute. Hence, when the band gaps of the two perturbed structures coincide in energy (i.e., when c 1 = c 2 ), the identity terms in Eqs. (4) and (5) can be absorbed into the unperturbed Hamiltonian, leaving us with two anticommuting mass terms gapping the K/K valleys. Note that while the enlarged unit cell is not primitive in the case of inversion-breaking M 1 , choosing such a cell allows both perturbations to be combined in a single lattice. To complete the construction, we study the dispersion away from the two valleys; quasimomenta away from K/K similarly result in perturbative terms, entering Eq. (1) as [42,48,49] The two matrices γ 1 and γ 2 anticommute with each other and with both Γ 1 and Γ 2 . As the symmetry-breaking perturbations may be position-dependent, we have effectively arrived at the term Γ·m(r) and fully identified the ingredients for realizing the Jackiw-Rossi Hamiltonian in a continuous 2D medium. We highlight that our analysis is applicable to a broad range of geometries and systems, as it relies solely on the symmetries of the perturbations and on time reversal. Moreover, our construction does not rely on introducing pseudo-spin or pseudo-time reversal symmetries [15]. Instead, the initial 4-fold degeneracy is directly obtained by starting with the original unit cell and its associated space group. We can thus also directly verify that our construction falls into the BDI class: in our chosen 4D basis, the time-reversal operator reads T = IK, where I = τ 1 ⊗ 1 2 corresponds to the spatial inversion operator and K is the complex conjugation operator; chiral symmetry is given by C = τ 2 ⊗ 1 2 . The Hamiltonian obeys these symmetries due to the spatial symmetries required by our construction. Since C 2 = T 2 = 1, our system belongs to the BDI class [35,36] and thus admits a topological winding number [cf. Eq. (2)].
The signs of m 1 and m 2 distinguish four topologically distinct phases [48]. We offer an intuitive view of this distinction by examining the induced band-inversions in the expanded unit cell of the simple hexagonal lattice. We focus on the symmetries at the Γ point and consider the point group C 3v . We also choose bases for two instances of the C 3v twofold degenerate representation, e.g., the pairs of orbitals {p x , p y } and {d xy , d x 2 −y 2 }. These can mix freely under C 3v symmetry, and are guaranteed to form a fourfold degeneracy when the original unit cell is intact. Under the perturbation M 1 (Fig. 2, upper panel), the two distinct orientations of the enlarged unit cell are mapped into each other by inversion, fixing the eigenspaces {p x + d xy , p y + id x 2 −y 2 } and {p x − d xy , p y − id x 2 −y 2 }; these match realizations of the QVHE [50,51]. Under the perturbation M 2 (Fig. 2, lower panel), the enlarged unit cell remains inversionsymmetric, so that its eigenspaces must consist of even and odd functions, here {p x , p y } and {d xy , d x 2 −y 2 }, matching the crystalline realizations of the QSHE [15].
Bulk signatures.-Our discussion has so far relied on general symmetry arguments. At the same time, we show in Fig. 1 examples of lattices and their deformations. In each structure, the magnitudes of the mass terms m 1 and m 2 are evaluated numerically as they depend on the specifics of the symmetry-breaking deformations and material/platform realization. To extract the mass terms from the numerics, we focus on the Γ point of the enlarged unit cell and obtain the four numerical solutions H |ψ i = E i |ψ i , with energies E i and corresponding deviations from mid-gap ∆E i = E i −Ē forĒ = i E i /4. In the effective 4 × 4 model, the mass terms are easily evaluated, since m i = ψ j |Γ i |ψ j . However, as seen in Fig. 2, the numerical solutions are obtained in a realspace basis. We therefore need to convert the matrices Γ 1 and Γ 2 to real space operations, i.e., invert the function ρ(g) in Eq.  then yields where we find using Eq. (3) In Fig. 3, we plot the dependence of the effective masses m i on the structural features of the simple hexagonal lattice example, cf. Fig. 1(I). Here, inversion symmetry is broken by deforming the hexagonal incision into a triangle, parameterized by t. Translation symmetry is broken by introducing two different radii r 1 and r 2 for the central and outer hexagonal incisions, parameterized by r 1 /r 2 . Using such a construction, the deformations can be combined in a single unit cell, cf. Figs. 3(1)-(4). The calculations are performed using COMSOL [52], and the mass terms are evaluated using Eq. (7). Crucially, the mass terms are varied independently by the different symmetry-breaking deformations, see Figs. 3(a) and (b). This confirms that we can readily realize a spatiallydependent mass term m(r), as required by the Jackiw-Rossi model. We reiterate that to create a bound mode, the structures used for the winding must have coincident band gaps, analogously to c 1 = c 2 in Eqs. (4) and (5).
We move now to construct a Jackiw-Rossi vortex in a large supercell composed of the different bulk phases. To obtain a topological bound mode, the mass vector must wind [cf. Eq. (2)], i.e., we need to spatially pass through all four 'pure' perturbations shown in Fig. 2. Constructing a vortex of n = 1 in a dielectric medium, we numerically confirm the existence of a bound mode in Fig. 4. To assure an overall band gap, we use only the four unit cells shown in Fig. 2. Since we use periodic boundary conditions, the supercell features four distinct vortex sites, all with the same |n| but with different terminations. Correspondingly, we find four bound modes, each localized at one of the four vortices. By construction, the topological vortex modes should appear mid-gap. One of the vortex modes indeed lies at the bulk gap center, while the remaining three deviate slightly away [see Fig. 4(a)]. We attribute this discrepancy to the inequivalent vortex terminations, where the spatial symmetries are locally broken [39]. Note that gapped topological edge states exist at the domain walls between different mass regions, as expected in second-order topological insulators [17,18,24], see Figs. 4(a) and (c).
Discussion.-The presented symmetry principles and formulas for numerical calculations enable a systematic exploration of a plethora of structures, which feature high-order topological defects. The range of compatible structures is far broader than with the TB approach, e.g., the primitive cell of the simple hexagonal (Kagome) lattice has 1 (3) incisions, which precludes the formulation of a 4-band TB model required by a 2D HOTI. At the same time, we emphasize that the symmetry arguments describe the behavior of the bands near the high-symmetry points of the BZ, but do not guarantee a fully gapped spectrum. In the Kagome lattice, for example, a dielectric patterned with air incisions shows the correct symmetry breaking at the Γ point, but the spectrum is gapless due to the band dispersion at other regions of the BZ, yielding only an indirect gap at the Γ point. Such limitations are peculiar to specific experimental platforms. From a technological standpoint, the bulk band gap size is the key characteristic to keeping the vortex modes spectrally isolated and spatially confined. In Fig. 4(a), its size relative to the gap center is 16.6%, which is a factor of 3 higher compared to existing works based on the graphene-inspired TB model [39]. Further improvement is expected with the use of numerical optimization enabled by our work. We also note that the spatial confinement may be strengthened by nonlinear interactions [53,54]. Finally, we highlight that our scheme can be extended to three dimensions by appropriately extruding the geometry in the out-of-plane direction. This yields one-dimensional, counter-propagating topological modes that can be used to construct optical fibers [55,56].
This work was supported by the Swiss National Science Foundation through grants CRSII5 177198/1 and PP00P2 163818. We also thank I. Petrides, M. Rechtsman, T. Wolf, A. Eichler, and R. Chitra for helpful discussions.
Supplemental Material for "Second-order topological modes in two-dimensional continuous media" In this section, we provide a minimalist overview of the formalism needed to construct the representation matrices introduced in Eq. (??) of the main text. We follow the nomenclature used in Ref. [? ], where our exposition is simplified by focusing on symmorphic space groups.

A. Glossary
Symmorphic space group: The group G generated by a set composed of point group operations and real space translations is called symmorphic [? ].
Isogonal group: The point group operations of a symmorphic space group G form its isogonal group F.
Equivalent k-points: Two points in k-space are equivalent if they differ by a combination of reciprocal lattice vectors.
Little co-group: For a space group G, the operations which take a k-space vector k into an equivalent point form the little co-group of k, denoted G k ⊂ F.
Star: The set of inequivalent k-points generated from k by all operations of a group G is called the star of k.
Representation: The representation ρ of a group G is a homomorphism of G onto a group T of non-singular linear operators acting in a finite-dimensional vector space V .
Irreducible representation (irrep): If V has no proper subspace invariant under the operations of T except V itself, then ρ is an irrep.
Basis of an irrep: A set of functions B that spans the domain V of an irrep ρ is a basis of ρ. Each basis belongs to one irrep only.
Representation matrix: Having chosen a basis B, every symmetry element g ∈ G is mapped by ρ onto a matrix ρ(g) ∈ T . This is the representation matrix of g in the basis B.

B. Constructing an irrep and its basis
We now consider a space group G with the isogonal group F, and seek a basis B = {φ j } for a yet-unknown irrep ρ. As we deal with periodic metamaterials, Bloch's theorem sets the form for the individual basis functions, where k is a k-space vector, and u j (r) is invariant under lattice translations. In such a basis, the representation matrix of translation by a unit cell, ρ(T ), is diagonal. Choosing a vector k, we now proceed in three steps.
(i) Determine G k . For every g ∈ G k , we have g u j e ik·r = e ik·r gu j by definition, up to a periodic function which can be absorbed into the definition of u j . To be closed under G k , the set {u j } must therefore form a basis for its representation. Let us then choose a bdimensional irrep ρ k of G k ., and a basis B ρ k of ρ k such that every u j ∈ B ρ k is invariant under lattice translations. Each u j ∈ B ρ k can be used to define u j e ik·r ∈ B. We have thus constructed b basis functions of the irrep ρ.
(ii) Take an element h ∈ F that is not in G k . As such, h maps a Bloch function labeled by k into another member of its star, h e ik·r = e ik ·r . The transformed basis h (B ρ k ) is still a basis of ρ k , although generally B ρ k = h (B ρ k ). From each function h (u j ) ∈ h (B ρ k ), we again generate a h (u j ) e ik ·r ∈ B.
(iii) Repeat step (ii) for each inequivalent k generated by h.
For a k with a q-membered star and a b-dimensional irrep of G k , we have thus generated a (qb)-dimensional basis for an irrep ρ of G. In the band diagram, this appears as q degeneracy points, each b-fold degenerate. The representation matrix of each operation in G is found explicitly by its action in the basis B.
The application to the K point of the space group p6m is now straightforward. The isogonal point group of p6m is C 6v . The little co-group G K is C 3v and the star of K consists of two points, {K, K }. Seeking a fourfold degeneracy, we focus on the 2-dimensional irrep of C 3v . A possible basis for this irrep is (x, y); we therefore define the functions (p x , p y ), which transform as (x, y) under point group operations, but are invariant under lattice translations. Applying steps (i) and (ii) with h = C 6 gives four basis functions of p6m, {p x e iK·r , p y e iK·r , p x e iK ·r , p y e iK ·r } . (S2) Let us now find the representation matrices of the symmetry elements. The translation T is diagonal, containing phase shifts due to the K/K labels, ρ(T ) = diag{e 4π/3 , e 2π/3 } ⊗ 1 2 . (S3) The mirror plane σ x / ∈ G K maps (p x , p y ) → (p x , −p y ) and (K, K ) → (K , K), and thus gives ρ(σ x ) = σ 1 ⊗ σ 3 .