Photonic topological phases in dispersive metamaterials

We analyze the photonic topological phases in dispersive metamaterials which satisfy the degenerate condition at a reference frequency. The electromagnetic duality allows for the hybrid modes to be decoupled and described by the spin-orbit Hamiltonians with pseudospin 1, which result in nonzero spin Chern numbers that characterize the topological phases. In particular, the combined Hamiltonian of the hybrid modes complies with a fermionic-like pseudo time-reversal symmetry that ensures the Kramers degeneracy, leading to the topological protection of helical edge states. The transverse spin generated by the evanescent surface waves is perpendicular to the wave vector, which exhibits the spin-momentum locking as in the surface states for three-dimensional topological insulators. The topological properties of the helical edge states are further illustrated with the robust transport of a pair of counterpropagating surface waves with opposite polarization handedness at an irregular boundary of the metamaterial.

Kramers degeneracy may even exist in 2D dielectric photonic crystals 20,42,43 without the hybridization of TM and TE modes, where the pseudospin states are represented by the combination of p and d orbital-like basis functions. In photonic systems, the Kramers degeneracy is guaranteed by the pseudo TR symmetry 20,21 , which can be constructed from the symmetry in crystal structure 20,44 or constitutive relation 21 .
In the present study, we analyze the photonic topological phases in dispersive metamaterials which satisfy the degenerate condition (ε μ = − ) at a reference frequency ω 0 , around which a frequency gap between the bulk modes may exist. The electromagnetic duality allows for the hybrid modes (E ± η 0 H) to be decoupled at the reference frequency, which are determined by two subsystems with degenerate eigenvalues. By introducing the pseudospin states as the eigenfield basis, the hybrid modes are described by the spin-orbit Hamiltonians with pseudospin 1, which result in nonzero spin Chern numbers that characterize the topological phases. In particular, the combined Hamiltonian of the two subsystems is TR invariant under a fermionic-like pseudo TR operator T p with = − T I p 2 , which ensures the Kramers degeneracy of the hybrid modes, leading to the topological protection of helical edge states. The Kramers doublet consists of a pseudospin state from the plus hybrid mode (E + η 0 H) and another state from the minus hybrid mode (E − η 0 H), which are TR partners and orthogonal to each other.
For illumination, the surface modes at the interface between vacuum and the metamaterial are analytically formulated based on Maxwell's boundary conditions, which are represented by a dispersion surface in the frequency gap and reduced to an ellipse or a circle at the reference frequency. The evanescent surface wave generates a transverse spin perpendicular to the wave vector, which exhibits the spin-moment locking as in the surface states for 3D topological insulators. The topological properties of the surface modes are further illustrated with the electromagnetic radiation excited by an appropriately phased point dipole at the interface. The surface waves with opposite polarization handedness counterpropagate toward different directions, which are able to bend around sharp corners without backscattering.

Results
Bulk modes. Consider a dispersive medium characterized by the frequency-dependent permittivity tensor ε ε ω ( ) 0 and permeability tensor μ μ ω ( ) 0 . Treating the combined electric field E = (E x , E y , E z ) and magnetic field H = (H x , H y , H z ) as a six-component vector, Maxwell's equations for the time-harmonic fields (with the convention e −iωt ) are written in matrix form as where I is the 3 × 3 identity matrix, H′ = η 0 H, and η μ ε = / 0 0 0 . Assume that the medium is uniaxially anisotropic with the material parameters: ε ε ε ε = diag( , , ) t t z and μ μ μ μ = diag( , , ) t t z . The existence of a nontrivial solution of E and H′ requires that the determinant of the 6 × 6 matrix in Eq. (1) be zero, which gives the characteristic equation of the bulk modes as 2 and k 0 = ω/c. The above bi-quadratic equation shows that the bulk modes consist of two parts with dual symmetry between ε n and μ n (n = t, z).
For analyzing the topological phases in the present medium, we assume that the medium satisfies the 'degenerate' condition: ε μ = − at a reference frequency ω 0 . In the neighborhood of ω 0 , ε n and μ n can be approximated, respectively, as ε ε ω ω ε ε δω ω , a frequency gap between the two bulk modes exists around ω 0 . The band gap size can be estimated by the solutions of ω at k t = k z = 0 [cf. Eq. (2)], between which the bulk modes do not exist. For ε < 0 n0 (μ > 0 n0 ), the upper and lower band edges are approximated by ω ω and μ  z are assumed to be positive definite 45 .
Spin-orbit Hamiltonians. The electromagnetic duality of Maxwell's equations dictates that the matrix in Eq. (1) is symmetric with the same diagonal elements when the degenerate condition is satisfied. This enables us to rewrite the wave equations for the hybrid modes, defined by E ± H′, at the reference frequency as . Note that the hybrid modes are completely decoupled and determined by two subsystems (3 × 3 matrix) with the same eigenvalues, their matrix determinants being equal: , the wave equation for the plus hybrid mode F = E + H′ can be rearranged as (see Methods A) by introducing the pseudospin state ψ ψ where x y z , and S n (n = x, y, z) are the spin matrices for spin 1.
Note here that Eq. (4) is formulated as an eigensystem with δω as the eigenvalue. The Hamiltonian +  in Eq. (5) represents the spin-orbit interaction S⋅κ with pseudospin 1 45 , which is mathematically equivalent to the Hamiltonian of a magnetic dipole moment in a magnetic field. The imaginary wave vector κ in  + reflects the fact that the wave is evanescent in the frequency gap. For the minus hybrid mode F = E − H′ in Eq. (3), it is straightforward to show that the corresponding spin-orbit Hamiltonian is given by Note that the hybrid modes can be slightly different if the corresponding Hamiltonian is rewritten in a different format 37,46,47 . Topological invariants. The topological properties of the spin-orbit Hamiltonians  ± can be characterized by the topological invariants based on their eigenfields. For this purpose, we calculate the Berry flux over a closed surface S: 2 , corresponding to the bulk mode at the reference frequency ω 0 in the imaginary wave vector space. The eigensystem for the plus hybrid mode in Eq. (5) ψ λ ψ = .
σ σ σ + (7)  is solved to give the eigenvalues λ σ and eigenvectors ψ σ (σ = ±1, 0), based on which the Chern numbers are calculated to give C σ = 2σ (see Methods B). The nonzero topological invariants C σ (σ = ±1) reveal the topological nature of the pseudospin states, where σ describes the helicity (or handedness) of the states. For this subsystem, the total Chern number = ∑ σ σ C C and the spin Chern number σ = ∑ σ σ C C spin 33 are given by spin which are consistent with the quantum spin Hall effect of light 48 . On the other hand, the eigensystem for the minus hybrid mode in Eq. (6) is given by where the eigenvalues λ σ are the same as those of  + . The helicity of the eigenvectors, however, has been flipped from σ to −σ. The Chern numbers are therefore change signs as C σ = −2σ. For this subsystem, the total and spin Chern numbers are given by spin The vanishing total Chern number C in the subsystems reflects the TR symmetry of Maxwell's equations and the absence of QH states in free photons, while the spin Chern numbers C spin = ±4 indicate that there exist two pairs of QSH edge states which are doubly-degenerate with respect to the helicity σ. The existence of surface modes in Maxwell's equations, however, requires the presence of an interface (between two different media) that breaks the duality symmetry of electromagnetic fields as in an unbounded region, and therefore only one pair of edge modes survives at the interface 48 . This feature will be confirmed by the characteristic equation of surface modes at the interface between vacuum and the metamaterial.
According to Eqs (7) and (9), the eigensystem of the combined hybrid modes is written as which states that the eigenstates ψ σ and ψ −σ of their respective subsystems are degenerate with the eigenvalue λ σ . The combined Hamiltonian is therefore considered two copies of the spin-orbit Hamiltonian with opposite helicity. This feature will play a crucial role in constructing the Kramers degeneracy in the present problem, leading to the topological protection of helical edge states.

Pseudo time-reversal symmetry. The Hamiltonian for Maxwell's equations [cf. Eq. (1)] in a lossless
is the bosonic TR operator for photons and K is the complex conjugation 18 . The is the fermionic TR operator for electrons 18 . Therefore, the Kramers degeneracy does not hold in a general photonic system, unless other symmetry such as polarization degeneracy 15,21 or spatial symmetry 20 has been imposed in the system. In another aspect, the combined Hamiltonian of the hybrid modes for Maxwell's equations with the degenerate condition: ε μ = − (at the reference frequency ω 0 ) is TR invariant under T p , that is, where T p is the fermionic-like pseudo TR operator having the same form of T f . The pseudo TR operator T p is inspired by noticing that E + H′ ↔ E − H′ during the TR operation. The pseudo TR operator is thus defined as are the Pauli matrices. The pseudo TR symmetry of the combined Hamiltonian  c ensures the Kramers degeneracy and guarantees the appearance of a Kramers doublet.
As revealed in Eq. (11), the photonic Kramers doublet consists of an eigenstate ψ σ of  + and another eigenstate ψ −σ of  − , with the same eigenvalue λ σ (σ = ±1). The states ψ σ and ψ −σ become TR partners under T p , that , which implies that it is impossible to introduce any backscattering between the two states, unless the TR symmetry has been broken. The two states therefore counterpropagate toward opposite directions without backscattering, a typical feature of the helical edge states that appear in the QSH system. As indicated in the field basis ψ [cf. Eq. (4)], the photonic Kramers doublet is a pair of two pseudospin states for the hybrid modes, analogous to the spin-up and spin-down states in electronic systems. The nonzero Chern numbers C σ associated with the eigenstate ψ σ further assert that the helical edge states are topologically protected, their existence being guaranteed by the difference of topological invariants on two sides of the interface. As the topological invariants remain constant under arbitrary continuous deformations of the system, the topological properties of the isotropic medium will be preserved when a certian anisotropy is included in the medium. The exact calculation of topological invariants for the anisotropic medium may resort to the numerical integration of Berry curvatures 49 . Surface modes. Let the xz plane be the interface between a dielectric with the relative material parameters ε d , μ d and a uniaxially anisotropic medium characterized by ε t , ε z , μ t , and μ z . The characteristic equation of surface modes is formulated based on Maxwell's boundary conditions: the continuity of tangential electric and magnetic field components at the interface, which is given by (see Methods C) ( 2) 0 2 2 2 are the normal (to interface) wave vector components in the dielectric, are the normal components in the anisotropic medium. For the surface waves to be valid on the dielectric side (y > 0), k y (1) and k y (2) should be purely imaginary with a positive value, which requires that ε μ 2 . On the anisotropic medium side (y < 0), k y (3) and k y (4) should be purely imaginary with a negative value, which requires that In the presence of square roots in k y n ( ) (n = 1, 2, 3, 4), the surface modes are represented by a part of the bi-quadratic surface in the frequency-wave vector space. In the isotropic case, where ε t = ε z ≡ ε and μ t = μ z ≡ μ, Eq. (14) is simplified to εε εμ με ε ε μμ με εμ depending on whether ε ω < ( ) 0 0 or μ ω < ( ) 0 0 , respectively. In this situation, the characteristic equation represents a quadratic surface of revolution about the frequency axis.
At the reference frequency ω 0 , where ε t = −μ t and ε z = −μ z , Eq. (14) is reduced to which is a part of the bi-quadratic curve. If ε d = μ d , Eq. (16) is further simplified to a standard quadratic curve as . Here, we assume that μ t and μ z are of the same sign. In the isotropic case, Eq. (15) at the reference frequency ω 0 is simplified to ε ε ε μ ε ε μ μ ε μ depending on whether ε < 0 or μ < 0, respectively. In this situation, the characteristic equation represents a circle when ε ε > d 2 2 or μ ε > d 2 2 , which is valid even when ε d ≠ μ d . In case all the materials considered are nonmagnetic, that is, μ = 1 (ε = −1) and μ d = 1 at ω = ω 0 , Eq. (16) is simplified to , which represents a circle when ε < < 0 1 d . For a fixed k z (k x ), Eq. (16) allows for two solutions of k x (k z ) with reflection symmetry about the k z (k x ) axis, which means that there exist a pair of surface modes at the interface counterpropagating toward the positive and negative k x (k z ) directions. Figure 1 shows the dispersion of bulk modes in the frequency-wave vector space for the dispersive medium based on Eq. (2). The Drude-type dispersive model, which is usually employed in the analysis of metamaterials, is assumed for the permittivity and permeability components: . If ε z = 0 or ε t = 0 at the reference frequency ω 0 , the upper and the lower modes touch at ω = ω 0 [ Fig. 1(a)]. In particular, the upper mode has a positive index (along the optical axis): ε μ ≡ > n 0 z z z , while the lower mode has a negative index: < n 0 z . If both ε t ≠ 0 and ε z ≠ 0 at ω = ω 0 , a frequency gap is opened between the two bulk modes, with the gap size dependent on the material parameters [ Fig. 1(b)]. As the anisotropy of the medium is reduced, the bulk modes tend to be dispersionless, that is, independent of the frequency. In the isotropic case, where ε t∞ = ε z∞ ≡ ε ∞ and μ t∞ = μ z∞ ≡ μ ∞ , the bulk modes are basically represented by two flat surfaces at ω ω ε = ∞ / p and ω ω μ = ∞ / p , leaving in between a frequency gap when ε ∞ ≠ μ ∞ . Figure 2(a) shows the dispersion of surface mode (in red color) at the interface between vacuum and the dispersive metamaterial based on Eq. (14), which exists in the frequency gap between two bulk modes (in gray color). The surface mode is represented by a funnel-shaped surface with a lower frequency at the center that connects to the lower bulk mode. As k x or k z increases, the surface mode raises its frequency and approaches asymptotically to

Discussion
, the lowest frequency of the upper bulk mode is given by . At the reference frequency ω 0 , where the degenerate condition ε = −μ is satisfied, the surface mode is described by an ellipse (denoted by black curve) with ε d = μ d and μ μ μ μ > > t tz d 2 2 [cf. Eq. (17)]. In the isotropic case where ε t∞ = ε z∞ = ε ∞ , the surface mode approaches asymptotically to a flat surface for a sufficiently large k x or k z , as shown in Fig. 2(b). The asymptotic frequency becomes ω ω ε Eq. (19)], which has the same form of the surface plasma frequency at the interface between a dielectric and the Drude-type metal. In this situation, there is always a gap between the surface mode and the upper bulk mode. At the reference frequency ω 0 , the surface mode is represented by a circle [cf. Eq. (18)].
Since the surface waves are evanescent in the direction normal to the interface, their normal wave vector components are purely imaginary, that is, k y = ±iκ y and κ y is real. The transversality condition (k ⋅ E = 0) dictates that  the normal (to interface) electric field component E y has a π/2 phase difference relative to the tangential component E x or E z . As the wave propagates at the interface (y = 0), the electric fields rotate in the perpendicular (xy or yz) plane, leading to elliptically polarized waves. The elliptical polarization of the surface waves is demonstrated in Fig. 3(a), where the instantaneous electric fields (blue arrows) and magnetic fields (green arrows) at the interface result in helical trajectories (traced by the tips of field vectors) along the propagation direction (the helical trajectory for the magnetic fields is less obvious as the normal magnetic field componet H y is relatively smaller than the electric field component E y ). In particular, the rotating electric field generates a transverse spin S ⊥ perpendicular to the wave vector k as 54 2 which is considered as the spin-momentum locking in the evanescent surface waves. The transverse spin of the surface mode exhibits a vortex spin texture [denoted by white arrows in Fig. 2(b)] that occurs in the surface states for 3D topological insulators. For the evanescent surface wave with the transverse spin, the handedness can be evaluated by first calculating the electric field components in a new coordinate system (x′, y′, z′), obtained by rotating the original system (x, y, z)  about the y axis such that the x′ axis is oriented to the time-averaged Poynting vector on the xz plane. Denoting θ the angle from the x axis to the x′ axis, the electric field components in the new coordinate system are given by E x′ = E x cosθ + E z sinθ, E y′ = E y , and E z′ = −E x sinθ + E z cosθ. The polarization handedness is then determined by the phase difference δ between E y′ and E z′ : is the argument of a complex number z. The wave is right-handed elliptically polarized (REP) or left-handed elliptically polarized (LEP) if δ = π/2 or −π/2, respectively, that is, the phase of E z′ is delayed or advanced by 90° relative to that of E y′ (under the time-harmonic convention e −iωt ). In Fig. 3(b), the polarization handedness of the surface waves at the reference frequency ω 0 is shown to have odd symmetry with respect to either of the two in-plane wave vector components (k x and k z ). For a constant k z (indicated by black dashed line), there exist a LEP wave propagating toward the positive k x direction and a REP wave toward the negative k x direction. The two surface waves with opposite handedness counterpropagate at a given edge (y = 0 in the xy plane), showing the feature of spin-momentum locking that occurs in the helical edge states. This feature holds when the surface mode is described by either an ellipse (for the anisotropic medium) or a circle (for the isotropic medium).
Finally, the topological properties of the surface modes are illustrated with the electromagnetic wave propagation at the interface between vacuum and the metamaterial 55 , as shown in Fig. 4. Here, a dipole source with k z /k 0 = 1.2 is placed at the interface (marked by asterisk symbol) to excite the surface wave, where the field is evanescent both in vacuum (outside the dispersion circle: 2 ) and the metamaterial (inside the frequency gap). By appropriately adjusting the phase of the point dipole, the REP (LEP) surface wave at the reference frequency ω 0 propagates unidirectionally toward the left (right), which is consistent with the direction of surface wave energy flow [cf. Figure 3(b)] and exhibits the typical feature of spin-polarized helical edge states. In particular, the surface waves are able to bend around sharp corners without backscattering, showing the robust transport of edge states against disorder. As the frequency deviates from ω 0 , the degenerate condition (ε = −μ) will be violated to a certain extent. The removal of degeneracy, however, is tolerated to some degree for the system to remain in the photonic topological phase 15,23 .
In conclusion, we have analyzed the photonic topological phases in dispersive metamaterials with the degenerate condition at a reference frequency. The topological phases are characterized by the nonzero spin Chern numbers for the hybrid modes described by the spin-orbit Hamiltonians with pseudospin 1. In particular, the hybrid modes comply with a fermionic-like pseudo TR symmetry that ensures the Kramers degeneracy, leading to the topological protection of helical edge states. The transverse spin generated by the evanescent surface wave is perpendicular to the wave vector, which exhibits the spin-momentum locking as in the surface states for 3D topological insulators. The topological features of helical edge states are further demonstrated by the robust transport of surface waves at an irregular boundary between vacuum and the dispersive metamaterial.

Methods
Spin-orbit Hamiltonians. The wave equation for the plus hybrid mode F = E + H′ in Eq. (3) can be rewritten as  Figure 4. Electromagnetic wave simulation for the surface modes at the interface between vacuum and the dispersive metamaterial with the same paramsters in Fig. 3(b). Asteristic symbols denote the dipole sources for exciting the (a) REP and (b) LEP surface waves at the reference frequency ω 0 .
x y x y is the basis of the pseudospin states that include a π/2 phase difference in the transverse field components (with respect to the optical axis of the anisotropic medium) 45 . In the neighborhood of ω 0 , ε n is approximated as ωε ω ε ε δω ≈ +  n n n 0 0 (n = t, z) and Eq. (21) can be rearranged as where ε ε with α ε =  c/ , κ κ κ κ = + +ˆx y z x y z , = + +ˆŜ x S y S z S Topological invariants. In terms of spherical coordinates, the Hamiltonian for the plus hybrid mode where κ x = a sin θ cos φ, κ y = a sin θ sin φ, and κ z = a cos θ with ε = | | a k 0 . Here, θ and φ are the polar angle and azimuthal angle, respectively, on the closed surface S: κ κ κ ε + + = k x y z 2 2 2 2 0 2 , corresponding to the bulk mode at the reference frequency ω 0 in the imaginary wave vector space. The eigensystem of the plus hybrid mode