Möbius edge band and Weyl-like semimetal flat-band in topological photonic waveguide array by synthetic gauge flux

2.1 Theoretical model

After setting up the gauge field configurations in the 2D coupling lattices, one immediately observes that for the type-I configuration shown in Figure 1e, the primitive translation L _y is preserved, whereas primitive translation L _x is not preserved. While for the type-II configuration sketched in Figure 1f, both L _x and L _y are not preserved. To recover the original gauge pattern, an additional gauge transformation G must be incorporated [16]. Namely, under a gauge field, the proper primitive translation operator along the x/y direction is changed to L_x/y = G_x/yL_x/y. Here, for the type-I gauge configuration, the transformation G _x is associated with G₁ where sites at odd (even) rows are multiplied with a π (0) phase, and the transformation G _y remains a unit transformation that leaves the original primitive translation unchanged. Whereas for the type-II configuration, G _y corresponds to G₁, and G _x corresponds to G₂, whereby sites located at diagonal (off-diagonal) positions undergo multiplication with a π (0) phase. Then, under the Z2 gauge field, the anti-commutation relation between the rebuilt primitive translations is satisfied, i.e.,

A consequence for the type-I Z2 PRSG [see Figure 1e] is that each energy band is two-fold degenerated for a generic momentum k and a four-fold degenerate Dirac point appears at the high symmetry point M = (π, π) of the Brillouin zone (BZ) [19, 20]. In sharp contrast, for the type-II Z2 PRSG [see Figure 1f], the eigenvalues of the translation operators, L _x and L _y , are ±ieikx/2 and ±ieiky/2, respectively (see more in Supplementary Material Section 1). This is a significant hallmark of the Möbius insulator, wherein 4π periodicity is implied. Together with the time reversal symmetry (T) at the BZ center Γ = (0, 0), the dual eigenvalues are complex conjugated, contributing to Möbius twisted bands crossing at k_x/y = 0. Note that the twofold degeneracy at every k can also be understood by the projective PT symmetry, i.e., PT2=−1 with P=G2P (see Supplementary Material Section 2) [17, 18, 21].

To further elucidate the effect of the gauge field and the Z2 PRSG, the general model Hamiltonian can be written as

where tx1±=(h11±h21)/2, tx2±=(h12±h22)/2, ty1±=(v11±v21)/2, ty2±=(v12±v22)/2, Γ _μ with μ = 1, 2, …, 5 are the 4 × 4 Dirac matrices: Γ₁ = τ₃ ⊗ σ₂, Γ₂ = τ₂ ⊗ σ₀, Γ₃ = τ₁ ⊗ σ₀, Γ₄ = τ₃ ⊗ σ₁, Γ₅ = τ₃ ⊗ σ₃, with τ and σ being two sets of the Pauli matrices. The real hopping amplitudes h(v)₁₁, h(v)₁₂, h(v)₂₁, and h(v)₂₂ are indicated in Figure 1e and f. Under the Z2 PRSG, the Hamiltonian for the type-I gauge configuration schematically shown in Figure 1e is reduced to HIk=t1+coskxiΓ1Γ5−t1+coskyiΓ2Γ5−t⁡sinkyiΓ3Γ5+t⁡sinkxiΓ4Γ5, with tthe strength of all the equal hopping amplitudes, i.e., |h11=h12=|h21=h22=|v11|=v12=|v21|=v22=t.

For the type-II gauge configuration [Figure 1f], the Hamiltonian becomes different to the Hamiltonian HIk whose four-fold Dirac point locates at (k _x , k _y ) = (π, π) [19, 20], a four-fold degenerate point appears at (k _x , k _y ) = (0, 0) for HIIk as shown in Figure 2a. This agrees with the above theoretical analysis and represents a case of Dirac semimetal. The projected band for the corresponding finite-sized structure (e.g., truncated in the y direction) is shown in Figure 2b which clearly characterizes the gapless nodal point.

Figure 2:

Band structure and edge bands. The band structure of the Hamiltonian (a) HIIk, (c) HIIIk, (e) HIVk and (g) HVk. The corresponding bulk bands are represented by blue curves while the edge bands are shown in red. The open edges in the y direction are illustrated in (b), (d), (f), and (h), respectively. The color map in (d) shows the phase evolution of the eigenvalues in the Möbius-twist edge bands. In (g), the bulk bands are for a specific case θ = π/4 of the dimerization (|h11|=h22=t+t′⁡sin⁡θ, |h12|=h21=t−t′⁡sin⁡θ, |v11|=v22=t+t′⁡cos⁡θ, |v12|=v21=t−t′⁡cos⁡θ). The colored ring indicates the trajectory of the θ-dependent nodal point splitting.

Initiated with the four-fold Dirac semimetal, it is noted that selectively breaking the primitive translation L _x and/or L _y leads to a series of topological phase transition. For example, simultaneously breaking L _x and L _y yields the nontrivial quadrupole topological insulator. This has been demonstrated by setting |v₁₁| < |v₁₂|, |v₂₁| < |v₂₂|, |h₁₁| < |h₁₂| and h21<h22 to dimerize the intra- and inter-hopping intensity along both x and y directions [26–28]. On the other hand, breaking the primitive translation by dimerization of the hopping along only one direction would also enforce the bulk gap opening and give rise to novel topological physics. In the following, we focus on the type-II gauge configuration which could be realized in topological photonic Möbius waveguides as we show in the proceding section. The dimerization along y-axis is one of the simplest ways to achieve the desired symmetry breaking. In this way, the symmetry-constrained lattice model reads

where the strength of v₁₁ and v₂₁ (v₂₁ and v₂₂) is denoted as t₁ (t₂). The dimerization is introduced by setting t₁ ≠ t₂ which breaks the primitive symmetry L _y , while L _x is still preserved in this case. Its band structure is shown in Figure 2c which becomes gapped. In addition, the model has a sublattice symmetry S characterized by Γ₅, which anti-commutes with L _x , i.e., {Ŝ,L̂x}=0. In the eigenspace of L _x , the Hamiltonian HIIIk can be block diagonalized on account of its sublattice symmetry

where

with qky=t1−t2eiky and mkx=2t⁡sin(kx/2). Here, the first term in Equation (6) is nothing but the standard Su–Schrieffer–Heeger (SSH) model and the second term is a mass term. Note that for t₁ < t₂ the SSH model have topological nontrivial phase [29, 30]. Particularly, for k _x = 0, the system is decoupled to a pair of SSH models, which feature a pair of zero modes for an edge perpendicular to the y axis. Furthermore, there appear two Möbius-twist edge bands possessing opposite L _x eigenvalues ℓ±=±ieikx/2 which have a period of 4π and are inverted after wrapping around the BZ once, as shown in Figure 2d [corresponding to the case in Figure 1(b)]. Interestingly, we note that this case corresponds to the topological Möbius insulator that has been discovered in acoustics [19, 20]. The Z2 topological invariant is given by [16]

where F=∇kA is the Berry curvature with Berry connection A(k)=ψ−|i∂k|ψ− for the valence bands of h₁ featured by eigenstate |ψ−, and γkx is the Berry phase for subsystem h₁(k _x , k _y ) with fixed k _x . A nontrivial v indicates a Möbius topological insulator, with Möbius-twist edge bands at any L _x -invariant edge [26] (see more in Supplementary Material Section 3).

An alternating dimerization along the y-axis (i.e., v11=v22=t1<t2=v12=v21) would simultaneously breaks the primitive translation L _x and L _y , yielding another distinct topological phase transition. For that, the Hamiltonian of the lattice model is then changed to

In this case, the four-fold degenerate Dirac point splits into two nodal points with a two-fold degeneracy along the line k _y = 0, as shown in Figure 2e. This is in similar to the creation of a Weyl point pair from a four-fold degenerate Dirac point [31]. The topological charge can be determined by calculating the winding number around a circle C surrounding the nodal point: w=(1/4πi)∮Cdk⋅trΓ5H−1k∇Hk [19].

Figure 2e shows the complete band structure and Figure 2f shows that a flat edge band connecting the projections of the two nodal points, and the 2D counterpart of the Fermi-arc connecting the projection of the Weyl points, can be observed due to the nontrivial topological charge [32, 33]. To further break the hopping terms along the x direction, the position of the paired nodal points can be modulated in the full quadrants of the momentum space, instead of just being along the high-symmetric lines. To exemplify that, the dimerization along both the x direction (|h11|=h22=t+t′⁡sin⁡θ, |h12|=h21=t−t′⁡sin⁡θ) and the y direction (|v11|=v22=t+t′⁡cos⁡θ, |v12|=v21=t−t′⁡cos⁡θ) is applied and the Hamiltonian becomes

In this case, the paired nodal points are θ-dependent and behave as the colored ring shown in Figure 2g. Particularly for θ = π/4, Figure 2g shows that the bulk bands exhibit the linear crossing in the vicinity of the paired nodal points which fall on the momentum line ΓM. Due to the nontrivial topological charge, an edge band appears for an open edge either in the x or in the y direction, as shown in Figure 2h.

2.2 Full wave simulation and coupled model theory

In the previous section, we demonstrated the richness of topological phases and the accompanying edge modes. In this section, we move on to showcase the implementation and potential realization of these theoretically predicted phenomena in 3D photonic waveguide arrays that have the type-II gauge field. Figure 3a illustrates the geometry of the waveguide structure, composed of high-refractive index waveguides with circular and peanut-like cross-sections. Figure 3b displays a cross-sectional slice. It is important to note that the peanut-like waveguides are constructed by overlapping two circular waveguides. The radii and center-to-center distance can be adjusted to support eigenmodes with the same propagation constant for the fundamental mode in the individual circular waveguide (see Section 4 of the Supplementary Material). The peanut waveguides at opposite diagonal sites of a unit cell are tilted with respect to the y-axis by ± 45°, respectively. Unlike in acoustic lattices discussed in Refs. [19, 20], in the coupled photonic waveguide system, the coupling strength and sign are determined by the overlap of the decoupled base states. The coupling between the ith state Ei and the jth state |Ej is proportional to Ei|Ej [34, 35]. Figure 3c shows the fundamental modal electric field E _x (at λ = 760 nm) of the decoupled circular and peanut waveguides, in the form of monopole (s orbital) and dipole (p orbital), respectively. The results were obtained from the numerical simulations using the finite element method (FEM) in COMSOL Multiphysics^® 5.0. Clearly, p orbitals have ‘negative parts’ which contribute to the negative coupling. In the waveguide lattices shown in Figure 3b, the nearest-neighboring modal coupling can be represented by the constituent s and p orbitals, as shown in Figure 3d. Then, the positive (negative) couplings in the waveguide lattices are labeled by the solid (dashed) bonds. The coupling strength (the absolute value of the coupling coefficients) can be controlled by the site distance. The diagonal site couplings have been ignored in this work as they are relatively small (see Supplementary Material S6). By using the s and p orbitals, the π flux threading a plaquette is implemented [25]. Note that using s and d orbitals can implement the type-I gauge configuration shown in Figure 1a and e. Here in our case, each waveguide supports a pair of orbitals: s (s₁ and s₂) and p (p₁ and p₂) (see Supplementary Material Section 4). The coupling in the waveguide arrays can be in the form of s₁ − p₂ − p₁ − s₁ and s₂ − p₁ − p₂ − s₂, following the order A− B − C − D in Figure 3d.

Figure 3:

The designed photonic waveguide structures. (a) 3D schematic structure of the topological photonic Möbius waveguide array and (b) the cross-section view. (c) Fundamental mode (at λ = 760 nm) pattern Re(E _x ) of the individual circular waveguide and the peanut waveguide, resembling the s and p orbitals, respectively. (d) Illustration of the coupling behavior among the s and p orbitals in the waveguide lattices. The radius of the isolated circular waveguide is r₀ = 0.4625 μm. The radius and the center-to-center distances of the constituent circular parts forming a peanut waveguide are r₁ = 0.6 μm and d₁ = 1 μm, respectively. The zoom-in unit cell shaded in gray in (b) is shown at the bottom panel. The lattice constant is a = 6 μm. The distances between the waveguides in the x(y) directions, w₁₁, w₁₂, w₂₁, and w₂₂ (l₁₁, l₁₂, l₂₁, and l₂₂), can be varied to modulate the intra- and inter-hopping amplitude. The refractive index of the waveguide and the host are n _co = 1.59, and n _cl = 1.54, respectively. Note that solid and dashed line labeled in (d) indicates the positive and negative coupling constant due to the directional arrangement of the p orbital.

The propagating waveguide mode in the waveguide array can be strongly coupled and described by the Schrödinger-like wave equation. Within the coupled mode theory (CMT) the steady state equation reads [36, 37]:

where E(r, z) is the transverse electric field at the propagation distance z. The components of Ĥ can be obtained by the overlapping integral of the fundamental modal fields based on the CMT. Due to the consistency between the Hamiltonian Ĥ obtained by the CMT and H(k) described by Equation (2), a couple of topological phase transitions of the propagating modes are expected, which shall be intrinsic to these waveguide arrays (see Supplementary Material Section 5). Note that for the specific design in Figure 3, when there is no hopping/coupling dimerization (i.e., l₁₁ = l₁₂ = l₂₁ = l₂₂ = w₁₁ = w₁₂ = w₂₁ = w₂₂ = a/2, seen in the inset of Figure 4a), the propagation constant and the nearest-neighboring coupling between the waveguides are β₀ = 1.5514k₀ and t = 0.00027k₀, respectively, where k₀ is the wavenumber in free space at the working frequency (λ = 760 nm) (see Supplementary Material Section 6). Figure 4a shows that the bulk bands obtained by the CMT are in good agreement with the full-wave numerical data from the FEM. Clearly, the bulk bands are two-fold degenerated, and a four-fold Dirac point appears at Γ point. To break the primitive translation symmetry L _y , the dimerization in the waveguide arrays can be realized by altering the center-to-center distances of the individual waveguides, i.e., by setting l₁₁ = l₂₁ = a/2 + h and l₁₂ = l₂₂ = a/2 − h.

Figure 4:

Comparison between the band structures calculated using the full-wave numerical method (gray solid line) and the CMT (red dashed line) for: (a) the Dirac semimetal case with w₁₁ = w₁₂ = w₂₁ = w₂₂ = l₁₁ = l₁₂ = l₂₁ = l₂₂ = a/2, (b) the Möbius insulator case with w₁₁ = w₁₂ = w₂₁ = w₂₂ = a/2, l₁₁ = l₂₁ = a/2 + 0.03a, and l₁₂ = l₂₂ = a/2 − 0.03a, (c) the Weyl-like semimetal case with w₁₁ = w₁₂ = w₂₁ = w₂₂ = a/2, l₁₁ = l₂₂ = a/2 + 0.03a, and l₁₂ = l₂₁ = a/2 − 0.03a, and (d) the Weyl-like semimetal case with w₁₁ = w₂₂ = a/2 + 0.03a, and w₁₂ = w₂₁ = a/2 − 0.03a, l₁₁ = l₂₂ = a/2 + 0.03a, and l₁₂ = l₂₁ = a/2 − 0.03a. We note that a comprehensive phase diagram with respect to the geometric dimerization can be found in the Supplementary Material Table 1 (Section 5).

Firstly, we impose a staggered dimerization pattern by setting h = 0.03a (see the inset of Figure 4b). Since the primitive translation L _y is broken, a band gap is opened, as shown in Figure 4b. Here, the coupling along the y-axis obtained by the CMT are t₁ = 0.00015k₀ and t₂ = 0.0005k₀. In accordance with the tight-binding model Equation (4), the four-fold degenerate Dirac point becomes gapped, and two Möbius twisted topological edge bands in Figure 2d are expected to emerge. To observe these Möbius-twist edge states, an open boundary condition parallel to the x-axis is imposed on a finite structure truncated in the y direction. Figure 5a shows the band structures for a ‘ribbon-like’ structure with 10 periods in the y direction. It is seen that there appear four edge bands inside the band gap, with the edge states associated with L _x eigenvalues clearly indicated by the phase colorbar. As a hallmark of the Möbius topology, the phase of the eigenvalues exhibits a 4π periodicity. Since the composed waveguides intrinsically host two degenerate modes, two pairs of Möbius-twist edge states are observed at the upper boundary (waveguides in the top layer) of the ‘ribbon’ structure. The field patterns of the edge states at k _x = 0.2π are, respectively, shown in Figure 5d–g, and it is clear that these edge fields are primarily localized at the waveguides near the upper boundary, with different phases. This represents a key characteristic of the Möbius-twist edge bands. We note that actually there are also four Möbius edge bands for the lower boundary of the ‘ribbon’ structure. These edge band states residue in the lower boundary (waveguides in the lower layer) degenerate with those in the upper boundary (details are not shown).

Figure 5:

Simulated dispersions for supercell structures with open boundary condition in the y direction and the periodic boundary condition in the x direction. (a) Möbius insulator for w₁₁ = w₁₂ = w₂₁ = w₂₂ = a/2, l₁₁ = l₂₁ = a/2 + 0.03a, and l₁₂ = l₂₂ = a/2 − 0.03a, (b) Weyl-like semimetal phases for dimerization along the y direction: l₁₁ = l₂₂ = a/2 + 0.03a and l₁₂ = l₂₁ = a/2 − 0.03a, (c) Weyl-like semimetal phase for dimerization along both axis: l₁₁ = l₂₂ = a/2 + 0.03a, l₁₂ = l₂₁ = a/2 − 0.03a, w₁₁ = w₂₂ = a/2 + 0.03a, and w₁₂ = w₂₁ = a/2 − 0.03a. The eigenmode patterns at k _x = 0.2π as labeled in (a) and (b) are shown for the Mobius-twist edge states (d–g) and the Weyl-like edge states (h–k), respectively.

Secondly, to build different dimerization pattern, h = 0.03a and h = −0.03a are applied to the two neighboring columns of the unit cell simultaneously. This results in the 2D Weyl-like semimetal phase defined by Equation (8) and Figure 2e. In this case, the four-fold degenerate Dirac point is transformed into a pair of twofold nodal points along the momentum line Γ − X. The results are shown in Figure 4c. Again, the CMT results agree well with the numerical simulations. In analogy to Weyl point and Fermi arc [18, 31], [32], [33], these 2D nodal points host nontrivial topological charge, which leads to the edge band connecting the projections of the nodal points. The results in the waveguide arrays are shown in Figure 5b. It is seen that four edge bands (two degenerate ones marked by the purple curves and two isolated ones by the cyan curves) connecting the projections of the nodal points. As a matter of fact, the next-nearest-neighbor coupling between the off-diagonal sites cannot be completely neglected which gives rise to a small bulk gap between them (see Section 7 in Supplementary Material). Figure 5h–k shows the eigenmode patterns for these edge modes. Different to the Mobius-twist edge state, these Fermi-arc like edge states localized at the lower and the upper boundaries of the waveguide array are mainly confined at the circular waveguide and the peanut waveguide, respectively. Besides, a slight deviation of the edge bands is due to a slight mismatch between the fundamental mode propagation wavenumber β₀ of the isolated circular and peanut waveguides.

Finally, to further break the intrinsic symmetries in the waveguide lattices, we set w₁₁ = w₂₂ = a/2 + 0.03a, w₁₂ = w₂₁ = a/2 − 0.03a, l₁₁ = l₂₂ = a/2 + 0.03a, and l₁₂ = l₂₁ = a/2 − 0.03a as shown in the inset of Figure 4d. In this case, only the time-reversal symmetry is preserved, and two nodal points appear on the line ΓM. This corresponds to the case of θ = π/4 of the model shown in Figure 2g. Figure 5c shows the projected band structures and the edge states, which is similar to the case of Figure 5b. Both cases are able to generate 2D Weyl-like phase. Certainly, similar phenomena can be observed for another configuration with arbitrary θ (results not shown here).

2.3 Coupled mode dynamics

We have shown that the Möbius insulator and Weyl-like semimetal are characterized by distinct topological edge bands. It is possible to show the linear crossing of the Möbius-twisted edge bands at k = 0 and the flat characteristic of the Fermi-arc-like edge band can be demonstrated with full-wave simulation. In view of the good agreements between the CMT and the FEM results and for the sake of simplicity, we provide a coupled mode dynamics calculation based on Equation (10). Specifically, we construct an array comprised of 3 × 23 unit cells (see Figure 3) and launch a Gaussian beam P(x,ϕ)=exp(−(x−x0)2/Δ2)exp(ikxx+iϕ) impinging upon the top row waveguides to excite the edge state. Here x labels the waveguide position, x₀ is the beam center and Δ the beam width. A phase difference ϕ = arg(ℓ_±) − π was applied to the top layer waveguides within a unit cell to align each individual Möbius band with a fixed k _x [20]. Figure 6a–d show the results, illustrating four representative cases of beam dynamics for various combinations of k _x and ϕ. Due to the linear crossing and twist characteristics of the edge bands around k _x = 0, the propagating beams can be modulated to travel towards the right or left, or bifurcate, by adjusting the phase difference ϕ to π/2, 0, and −π/2, as depicted in Figure 6a–c. However, as shown in Figure 6d, when the launched beam momentum is set to k _x = π, the excited wave propagates almost along the z axis without obvious deflection. This clearly shows the π rotation of the twisting bands shown in Figure 1a and b. Moreover, due to the flatness of the Weyl-like semimetal edge band, the excited beam is restricted to propagation along the z axis, as demonstrated in Figure 6e and f.

Figure 6:

Wave dynamics inside the waveguide array excited by a incident Gaussian beam with various initial spatial phase difference, for (a)–(d) the Möbius twisted bands and (e)–(f) the Weyl-like semimetal flat bands: (a) k _x = 0, ϕ = π/2, (b) k _x = 0, ϕ = 0, (c) k _x = 0, ϕ = −π/2, (d) k _x = π, ϕ = 0, (e) k _x = 0, and (f) k _x = 0.1π.