Photonic topological semimetals in bianisotropic metamaterials

We analyze the photonic topological phases in bianisotropic metamaterials characterized by a chirality tensor with zero trace. We found that the strength of chirality component determines the topological character of the metamaterial. The underlying medium can be considered as a topological semimetal with the nontrivial band gap in the momentum space. The topological properties are described by the spin-orbit Hamiltonians with spin 1 and characterized by the nonzero topological invariants. In particular, photonic quantum Hall states exist when the longitudinal chirality component exceeds the permittivity, whereas photonic quantum spin Hall states are present when the chiral nihility occurs. Considering the dispersion in the frequency domain, the bianisotropic metamaterial is regarded as a photonic Weyl system that supports the Weyl points and Fermi arcs. The topological features are further illustrated with the robust transport of edge states at an irregular boundary of the metamaterial.

The surface modes that connect two distinct topological phases are analytically formulated at the interface between vacuum and the metamaterial based on Maxwell's boundary conditions. The surface modes are represented by dispersion surfaces in the frequency-wave vector space, which are reduced to surface-state arcs at a reference frequency. The topological features of the surface states are further demonstrated by the electromagnetic radiation excited by a point dipole placed at the interface. The surface waves propagate unidirectionally at an irregular boundary and are immune to backscattering.

Results
Bulk modes. Consider a bianisotropic metamaterial characterized by the constitutive relations: where ε, μ, ξ and ζ are the frequency-dependent permittivity, permeability, and magnetoelectric coupling tensors, respectively. 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 reciprocal and the material tensors are uniaxial: ε ε ε ε = diag( , , ) t t z , μ μ μ μ = diag( , , ) t t z , and ξ ζ γ γ γ γ = − = = i i i i diag ( , , ) t t z , where ε n , μ n , and γ n (n = t, z) are real quantities. The existence of a nontrivial solution of E and H requires that the determinant of the 6 × 6 matrix in Eq. (3) be zero, which gives the characteristic equation of the bulk modes as where k t 2 = k x 2 + k y 2 and k 0 = ω/c. This is a bi-quadratic equation that incorporates the coupling between transverse electric and transverse magnetic modes.
We now assume that the metamaterial satisfies the duality condition: ε μ = 16,20,21 . Equation (4) can then be factorized as which is a product of two quadratic equations. In this situation, the bulk modes are a combination of two decoupled modes with opposite signs of γ t and γ z . In the neighborhood of a reference frequency ω ref , ε n (n = t, z) can be approximated as ε ε ω ω ε ε δω ω where ε  n is positive definite 34 . We assume that the chirality parameters γ t and γ z vary smoothly around ω ref and can be treated as constants as a first-order approximation. We further propose that the chirality tensor γ has a zero trace so that γ z = −2γ t . In case ε t = ε z , the underlying metamaterial is equivalent to the pseudochiral medium 38,39 characterized by the magnetoelectric coupling tensors with zero diagonal and symmetric off-diagonal entries (see Methods A), which can be synthesized by the lattice of intersecting split-ring resonators 40 or mutually perpendicular Ω-shaped metal wires 41 .
Spin-orbit Hamiltonians. The electromagnetic duality of Maxwell's equations dictates that the matrix in Eq. (3) holds a symmetric pattern when the duality condition ε μ = is satisfied. This enables us to rewrite the wave equations for the hybrid modes, defined by F ± = E ± iH′, as Note that F + and F − are completely decoupled and determined by two subsystems (3 × 3 matrix) with a similar structure. In the isotropic case, where , and γ t = γ z ≡ γ, the wave equation can be rearranged as (see Methods B) www.nature.com/scientificreports www.nature.com/scientificreports/ by introducing the pseudospin state ψ where x y z , and S n (n = x, y, z) are the spin matrices for spin 1. Note that Eq. (7) is formulated as an eigensystem with δω as the eigenvalue. The Hamiltonian ± H in Eq. (8) represents the spin-orbit interaction k · S with spin 1, which is mathematically equivalent to the Hamiltonian of a magnetic dipole moment in the magnetic field 34,35 . topological invariants. The topological properties of the spin-orbit Hamiltonians ± H can be characterized by the topological invariants based on the eigenfields. For this purpose, we calculate the Berry flux over a closed surface S: 2 , corresponding to the bulk mode at the reference frequency ω ref in the wave vector space. The eigensystem for the Hamiltonian ± H in Eq. (8): is solved to give the eigenvalues λ ± σ and eigenvectors ψ ± σ (σ = ±1, 0), based on which the Chern numbers are calculated to give (see Methods C) The nonzero C σ (σ = ±1) characterize the topological properties of the system, where σ refers to the helicity (or handedness) of the pseudospin states. In particular, the helical edge states are topologically protected, which means that their existence is guaranteed by the difference in band topology on two sides of the interface. As long as the band gap is not closed, the topological invariants remain unchanged under arbitrary continuous deformations of the system. The topological properties of the isotropic medium will be retained when a certain anisotropy is included. For the anisotropic medium, the exact calculation of topological invariants can be obtained by the numerical integration of Berry curvatures 27 . In this system, the total Chern number = ∑ = The eigensystem of the combined hybrid modes is written as Another condition of degenerate eigenvalues occurs at ε = 0 (and δω = 0), which is known as the chiral nihility 43,44 . The eigenvalues of the Hamiltonian ± H are given by λ + σ = −λ − σ ≡ λ σ and the eigensystem is written as In Eqs. (11) and (12), the combined Hamiltonian consists of two copies of the spin-orbit Hamiltonian with opposite helicity, which is characteristic of the Bernevig-Hughes-Zhang (BHZ) model for the QSH system 3 . In particular, the combined Hamiltonian is TR invariant under T p (see Methods D): c and T p is the fermionic-like pseudo TR operator with T p 2 = −1 20 . Here, β = −1 for γ = 0 (ε ≠ 0) and β = 1 for ε = 0 (γ ≠ 0) [cf. Equations. (8), (11), and (12)]. The pseudo TR symmetry of the combined Hamiltonian H c is crucial in determining the QSH phases in photonic systems of spin 1, which allows the coexistence of counterpropagating spin-polarized edge states as in electronic systems 12 . For a nonzero γ (and ε ≠ 0), + H and + H are no longer degenerate and the pseudo TR symmetry is not preserved.
Surface modes. Let the xz plane be an interface between vacuum and the bianisotropic metamaterial 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 E) 2 are the normal (to interface) wave vector components in vacuum, and 2 are the normal wave vector components in the bianisotropic medium, with α ± = ε ± γ t and β ± = ε ± γ z . For the surface waves to exist on the vacuum side (y > 0), k y (1,2) should be purely imaginary with a positive value, so that the waves decay exponentially away from the interface. On the bianisotropic medium side (y < 0), k y (3) and k y (4) should be purely imaginary with a negative value for a similar reason.

Discussion
Bulk modes. Let the permittivity and transverse chirality component of the bianisotropic metamaterial be positive (ε t = ε z ≡ ε > 0, γ t > 0, and γ z < 0) at the reference frequency ω ref , without loss of generality. The bulk modes of the metamaterial can be classified into three phases: (I) If |γ z | < ε and γ t < ε, the bulk modes are described by two intersecting ellipsoids [cf. Equation (5)], with the major and minor axes given by ε γ ε γ ± ± k ( )( ) z t 0 and (ε ± γ t )k 0 , respectively. The equifrequency surfaces of the bulk modes in the wave vector space are shown in Fig. 1(a) for ε = 1.3 and γ z = −2γ t = −1.
In this phase, the present medium behaves like an ordinary anisotropic dielectric material. (II) If |γ z | < ε and γ t < ε, the bulk modes are described by an ellipsoid and a two-sheeted hyperboloid [cf. Equation (5)]. Two band gaps (one for k z > 0 and the other for k z < 0) exist between three well separated bulk modes. The gap size Δk z = 2γ t k 0 is determined by the band edges: k z /k 0 = ±(ε ± γ t ), between which the bulk modes do not exist. The equifrequency surfaces of the bulk modes are shown in Fig. 1(b) for ε = 1.3 and γ z = −2γ t = −1.5. In this phase, the present medium is equivalent to the chiral hyperbolic metamaterial 24 . The chirality components in the former (γ t γ z < 0) play a similar role of the permittivity components in the latter (ε t ε z < 0), giving rise to the same bulk dispersion. (III) If |γ z | < ε and γ t < ε, the bulk modes are described by two asymmetric two-sheeted hyperboloids [cf.
Equation (5)] with a band gap Δk z = 2|ε − γ t |k 0 in between. The equifrequency surfaces of the bulk modes are shown in Fig. 1(c) for ε = 0.2 and γ z = −2γ t = −3. In this phase, the present medium behaves like an anisotropic hyperbolic material. (IV) If ε = 0 (and μ = 0), which is referred to as the chiral nihility 43,44 , the two hyperboloids are degenerate with an identical band gap Δk z = 2γ t k 0 . The equifrequency surfaces of the bulk modes are shown in Fig. 1(d) for ε = 0 and γ z = −2γ t = −3. In this phase, the present medium is equivalent to the double hyperbolic metamaterial 26 . The chirality components at nihility in the former (ε = 0 and γ t γ z < 0) play a similar role of the permittivity and permeability components in the latter (ε t ε z < 0, μ t μ z < 0, and γ t = γ z = 0), giving rise to the same bulk dispersion.
Edges modes. Figure 2 shows the surface modes at the interface (xz plane) between vacuum and the bianisotropic metamaterial based on Eq. (15) for different phases stated in Sec. 3.1. For a small chirality parameter such that the system is in the phase (I), the bulk modes consist of two intersecting ellipsoids. As the common band gap between the bulk modes on both sides of the interface do not exist, there are no surface modes at the boundary of the metamaterial [cf. Fig. 2(a)]. For a larger chirality parameter such that the system is in the phase (II), one of the two ellipsoids is transformed into a two-sheeted hyperboloid. The change in topology of the bulk dispersion opens a nontrivial band gap, leading to the topological phase transition in the momentum space 45 . The topological invariants of the bulk modes are determined by integrating the Berry curvatures on the dispersion surface 27 , which give C = −2sgn(γ t ) for the ellipsoid and C = sgn(γ t ) for each sheet of the hyperboloid 40 .
There exists a single surface mode at the boundary of the metamaterial (adjacent to vacuum), which connects the hyperboloid dispersion of the metamaterial to the spherical dispersion (dashed circle) of vacuum [cf. Fig. 2(b)]. The surface mode lies inside the common band gap of the metamaterial and vacuum for either k z > 0 or k z < 0, which are antisymmetric with respect to the k x axis. If the chirality parameter changes sign, the surface mode flips to the other side as a mirror reflection with respect to the k z axis. In particular, the surface modes propagate unidirectionally at the boundary (either +k x or −k x axis), which is characteristic of the chiral edge states in the QH system, their existence being consistent with the bulk-edge correspondence 40 . In this situation, the bianisotropic medium can be regarded as a photonic analogue of the QH system. Note that the two subsystems (2019) 9:18312 | https://doi.org/10.1038/s41598-019-54523-1 www.nature.com/scientificreports www.nature.com/scientificreports/ in the combined Hamiltonian [cf. Equation (11)] are not degenerate in the presence of chirality, and the pseudo TR symmetry is broken in this phase. As there are no states available for backscattering, the chiral edge states are insensitive to disorder 5 .
For an even larger chirality parameter such that the system is in the phase (III), both ellipsoids are transformed into two-sheeted hyperboloids. There are a pair of surface modes at the boundary of the metamaterial (adjacent to vacuum) inside the band gap for either k z > 0 or k z < 0 [cf. Fig. 2(c)], which are in general asymmetric about the k z axis. A particular situation occurs when the system is in the phase (IV), where the two hyperboloids are degenerate and the pair of surface modes is symmetric [cf. Fig. 2(d)]. The characteristic equation of surface modes [cf. Equation (15)] is simplified to Note that the degenerate bulk dispersions correspond to two copies of the subsystems for the hybrid modes [cf. Equation (6)], which are represented by the spin-orbit Hamiltonians [cf. Equation (8)]. The topological invariants of the Hamiltonians are given by C ± = ±2 [cf. Equation (10)]. For a transition of the bulk dispersion from an ellipsoid to a two-sheeted hyperboloid, the closed surface in the wave vector space is cut open along the equator into two sheets, on which the Berry curvatures are flipped and the Chern numbers are equally split into half for each sheet: C ± = ±1 26 . In particular, the pair of surface modes with opposite spins counterpropagates at the same boundary, which is characteristic of the helical edge states in the QSH system, their existence being consistent with the bulk-edge correspondence. In this situation, the bianisotropic medium is regarded as a photonic www.nature.com/scientificreports www.nature.com/scientificreports/ analogue of the QSH system. The combined Hamiltonian respects the pseudo TR symmetry, leading to the topological protection of helical edge states in the photonic system.
Weyl system. We now consider the explicit dispersion of the bianisotropic medium in the frequency domain.
The Lorentz-type dispersive model, which is usually employed in the study of metamaterials, is adopted for the permittivity and permeability components: ε ε ω ω 26 , where ω 0 is the resonance frequency of the resonators and ω p is the effective plasma frequency of the medium. The chirality components are given by γ n = Ω γn ωω p /(ω 2 − ω 0 2 ), where Ω γn 2 = Ω μn 46,47 . Such a dispersion guarantees that the energy density of the present medium is positive definite (see Methods F). In the present study, the bianisotropic response can be modelled by metallic helices oriented along three perpendicular directions 39 . In particular, the handedness of the helix in z direction is flipped so that the corresponding chirality component changes sign. The size of z-directed helix is adjusted in order to satisfy the traceless condition of the chirality tensor (2γ t + γ z = 0). Figure 3(a) shows the bulk and surface modes of the bianisotropic metamaterial in the frequency-wave vector space. In the frequency range: ω > ω 1 , the bianisotropic medium is in the phase (I) as a dielectric system, where ω 1 is the frequency across which the transition between the phase (I) (two ellipsoids) and the phase (II) (an ellipsoid and a two-sheeted hyperboloid) occurs, that is, ε(ω 1 ) = μ n (ω 1 ) = |γ z (ω 1 )|. It is noted that one ellipsoid in (I) and the hyperboloid in (II) touch at a pair of Weyl points symmetrically displaced on the k z axis: (k t , k z ) = (0, ±(γ t − γ z )k 0 ) [cf. green dot in Fig. 3(a)], which resembles the crossing of valence and conduction bands in the Weyl semimetal 33 . The frequency dispersion represents a tilted Weyl cone 29 , corresponding to the transition between type I Weyl points with spherical or ellipsoid dispersion surfaces and type II Weyl points with hyperbolic dispersion surfaces 48 . The surface modes (between the bianisotropic medium and vacuum) at the Weyl point frequency form the Fermi-arc-like edge states that connect the Weyl points [cf. red line in Fig. 3(a)]. The characteristic equation of surface modes [cf. Equation (15)] at this frequency is simplified to b a d c www.nature.com/scientificreports www.nature.com/scientificreports/ which is a combination of two linear equations. In the frequency range: ω 2 < ω < ω 1 , the bianisotropic medium is in the phase (II) as a QH system, where ω 2 is the frequency across which the transition between the phase (II) (an ellipsoid and a two-sheeted hyperboloid) and the phase (III) (two asymmetric hyperboloids) occurs, that is, ε(ω 2 ) = μ n (ω 2 ) = γ t (ω 2 ). In this range, the surface modes are approximate flat surfaces that slightly change their orientations with the frequency [cf. orange surface in Fig. 3(a)]. In the frequency range: ω < ω 2 , the bianisotropic medium is in the phase (III) as a bi-hyperbolic system. The surface modes are a pair of approximate flat surfaces in a somewhat smaller frequency range, where a common band gap of the metamaterial and vacuum exists. In this range, the surface modes are intersecting approximate flat surfaces [orange and green surfaces in Fig. 3(b)]. At a particular frequency ω = ω 3 , the bianisotropic medium is in the phase (IV) as a QSH system, where ω 3 is the frequency at which the chiral nihility occurs, that is, ε(ω 3 ) = μ n (ω 3 ) = 0. The two bulk modes are degenerate [cf. thick gray curves in Fig. 3(b)], which form the nodal lines 49 at the chiral nihility frequency.
Finally, the topological features of edge states are illustrated with the surface wave propagation at the interface between vacuum and the metamaterial 24,27 , as shown in Fig. 4. A dipole source is placed at the interface (marked by dot symbol) to excite the surface mode in the band gap of the metamaterial but outside the light cone of vacuum, so that the fields are evanescent on both sides (see Methods G). In Fig. 4(a,b), the surface modes are excited at k z /k 0 = 1.2 [cf. green dot in Fig. 2(b)], which propagate unidirectionally along an irregular boundary and are able to bend around sharp corners without backscattering. The propagation direction is reversed when the sign of chirality parameter is changed, showing the chiral nature of edge states in the phase (II), where the metamaterial behaves like a QH system. In Fig. 4(c,d), the surface modes are excited at the same k z , which correspond to a pair of k x 's [cf. green and blue dots in Fig. 2(d)]. In particular, the surface waves counterpropagate at the boundary for different handednesses of circular or elliptical polarization, which are immune to backscattering from disorder. The robust transport of the surface modes exhibits the helical nature of spin-polarized edge states in the phase (IV), where the metamaterial is regarded a QSH system.
In conclusion, we have investigated the photonic topological phases in bianisotropic metamaterials characterized by a chirality tensor with zero trace. The underlying medium is regarded as a photonic Weyl semimetal that supports the Weyl points and Fermi arcs. In particular, the photonic QH and QSH states exist in the momentum gap of the bulk modes, which are analytically formulated based on Maxwell's equations. The topological properties are described by the spin-orbit Hamiltonians and characterized by the nozero topological invariants. The topological features are further illustrated with the robust transport of edge states at an irregular boundary of the metamaterial.

Methods
Bianisotropic and pseudochiral media. Consider a pseudochiral medium characterized by the magnetoelectric coupling tensor ξ ζ ′ = − ′ with zero diagonal and symmetric off-diagonal entries as www.nature.com/scientificreports www.nature.com/scientificreports/ The above matrices can be diagonalized through a similarity transformation: is the transformation matrix formed by the column eigenvectors of ξ′. Note that ε ε ε ε ′ = diag( , , ) and μ μ μ μ ′ = diag( , , ) remain unchanged under the transformation: ε ε ε ε ε = ′ = − P P diag( , , ) 1 and μ μ μ μ μ = ′ = − P P diag( , , ) 1 . The above transformation shows that the pseudochiral medium characterized by the zero-diagonal ξ′ ( x y z is equivalent to the bianisotropic medium characterized by the traceless ξ ( ξ In the present study, the magnetoelectric coupling matrix contains nonzero diagonal elements with opposite signs [cf. Equation (19)], indicating that the waves propagating along the optical axis (z) and its transverse directions (x and y) exhibit opposite optical activities in the underlying medium 40 .
Spin-orbit Hamiltonians. The wave equation for the hybrid modes F ± = E ± iH′ in Eq. (6) can be rewritten www.nature.com/scientificreports www.nature.com/scientificreports/ x y x y is the basis of the pseudospin states that include a π/2 phase difference between the transverse field components (with respect to the optical axis of the medium) 34 . In the neighborhood of a reference frequency ω ref , ε n (n = t, z) can be approximated as ε ε ω ω ε ε δω ω where ε  n is positive definite 34 . Equation (21) is rearranged as where ε ε topological invariants. In terms of the spherical coordinates, the Hamiltonian ± H [cf. Equation (28)] is rewritten as where k x = asinθcosφ, k y = asinθsinφ, and k z = acosθ with a = |ε ± γ|k 0 . Here, θ and φ are the polar and azimuthal angles, respectively, on the closed surface S: 2 , corresponding to the bulk mode at the reference frequency ω ref in the wave vector space. The eigensystem for the Hamiltonian ± H : is solved to give the eigenvalues λ ± σ = |d ± |σ (σ = ±1, 0) and the normalized eigenvectors as  Note here that the eigenvalue λ ± σ is related to δω in Eq. (7) as λ ± σ = d ± + δω. Based on Eqs. (32) and (33), the Berry connections The Berry curvatures F σ = ∇ × A σ are then given by is the bosonic TR operator for photons, and K is the complex conjugation 12 . The Hamiltonian H m , however, is not TR invariant under T f , that is, is the fermionic TR operator for electrons 12 . The combined Hamiltonian for the hybrid modes with the duality condition: ε μ = , nevertheless, is TR invariant under T p , that is,   Simulation. The simulation domain is on the xy plane and k z is the out-of-plane wave vector component, which is kept fixed in the simulation so that the eigenwaves possess the same k z 24 . The surface mode is excited at a point on the boundary of the metamaterial, which can be implemented experimentally by a dipole antenna 14,18 . For the dipole to serve as the source of circular or elliptical polarization, the dipole has two in-plane components and the phase difference in between is set to be π/2 or −π/2 to mimic the right-handed or left-handed wave 50 . In experiment, the dipole source will excite the electromagnetic radiation for all wave numbers, and the measurement results are projected along the k z direction. The surface modes with a particular value of k z can be found in the contour plot of measurement 51 .