Hamiltonian formulation

While metamaterials have traditionally been described in terms of macroscopic effective properties30,33,43, the importance of crystallinity is becoming increasingly apparent44. Therefore, to capture the essential physics related to complex non-local effects that arise from strong multiple-scattering45, here we study the properties of the cavity-embedded honeycomb metasurface by means of a microscopic Hamiltonian formalism. This allows us to clearly identify the distinct physical origins of the type-I and type-II Dirac points.

The full polariton Hamiltonian of this system reads H pol = H mat + H ph + H int , where the interaction Hamiltonian H int couples the matter and photonic subspaces whose free dynamics are governed by H mat and H ph , respectively. We employ the Coulomb gauge, where the instantaneous Coulomb interaction between the meta-atoms is incorporated within the matter Hamiltonian H mat , and the effects of the dynamic photonic environment—described by the transverse vector potential—are included through the principle of minimal-coupling46.

A schematic of a cavity-embedded honeycomb metasurface is depicted in Fig. 1. We model each subwavelength meta-atom by a single dynamical degree of freedom describing the electric-dipole moment associated with its (non-degenerate) fundamental eigenmode with resonant frequency ω 0 . These meta-atoms are then oriented such that their dipole moments point normal to the plane of the lattice. Furthermore, we consider subwavelength nearest-neighbor separation a such that the light cone intersects the Brillouin zone edge above ω 0 , ensuring the existence of evanescently bound, subwavelength polaritons. The strength of the Coulomb dipole–dipole interaction between neighboring meta-atoms is parametrized by Ω. Finally, the metasurface is embedded at the center of a planar photonic cavity of height L, where the cavity walls are assumed to be lossless and perfectly conducting metallic plates. Such a structure is imminently realizable across the electromagnetic spectrum from arrays of plasmonic nanoparticles to microwave helical resonators (see Fig. 1).

Fig. 1 Schematic of a cavity-embedded honeycomb metasurface. The honeycomb array of meta-atoms is composed of two inequivalent (A and B) hexagonal sublattices—defined by lattice vectors \({\mathbf{a}}_{1} = a( -\frac{\sqrt{3}}{2} , \frac{3}{2} )\) and \({\mathbf{a}}_{2} = a( \frac{\sqrt{3}}{2} , \frac{3}{2} )\)—which are connected by nearest-neighbor vectors e 1 = a(0, −1), \({\mathbf{e}}_{2} = a( \frac{\sqrt{3}}{2} , \frac{1}{2} )\), and \({\mathbf{e}}_{3} = a( - \frac{\sqrt{3}}{2} , \frac{1}{2} )\), where a is the subwavelength nearest-neighbor separation. Each subwavelength meta-atom is modeled as an electric dipole, oriented normal to the plane of the lattice. The honeycomb metasurface is then embedded inside a photonic cavity of height L, which is composed of two perfectly conducting metallic plates, enabling one to modify the photonic environment while preserving the lattice structure. This general model can be readily realized across the electromagnetic spectrum, from arrays of plasmonic nanorods to microwave helical resonators Full size image

Emergence of type-I Dirac points

The matter Hamiltonian within the nearest-neighbor approximation reads

$$H_{\mathrm{mat}}=\hbar\tilde{\omega}_0\sum_{\mathbf{q}} \left(a_{\mathbf{q}}^{\dagger} a_{\mathbf{q}}+b_{\mathbf{q}}^{\dagger} b_{\mathbf{q}}\right)+\hbar\tilde{\Omega}\sum_{\mathbf{q}}\left(f_{\mathbf{q}}b_{\mathbf{q}}^{\dagger}a_{\mathbf{q}}+\mathrm{H.c.}\right)\,,$$ (1)

where, for brevity, we have not presented the non-resonant terms (see Methods for derivation). In Eq. (1), \(\tilde{\omega}_{0}\) is the renormalized resonant frequency and \(\tilde{\Omega}\) is the renormalized Coulombic interaction strength due to the cavity-induced image dipoles (see Methods for their dependence on the cavity height). The bosonic operators \(a_{\mathbf{q}}^{\dagger}\) and \(b_{\mathbf{q}}^{\dagger}\) create quanta of the quasistatic collective-dipole modes that extend across the A and B sublattices, respectively, with wavevector q in the first Brillouin zone (see Fig. 2a). Finally, the function \(f_{\mathbf{q}} = \sum_{j=1}^{3} \exp ({\mathrm{i}}{\mathbf{q}} \cdot {\mathbf{e}}_{{j}})\) encodes the honeycomb geometry of the lattice with nearest-neighbor vectors e j (see Fig. 1).

Fig. 2 Evolution of the polariton dispersion as the cavity height is reduced. a First Brillouin zone defined by primitive reciprocal lattice vectors \({\mathbf{b}}_{1} = \frac{2\uppi}{3a} \left( - \sqrt{3} , - 1 \right)\) and \({\mathbf{b}}_{2} = \frac{2\uppi}{3a}\left(\sqrt{3} , - 1\right)\). b Quasistatic dispersion of the collective-dipole normal modes, where the upper band corresponds to a bright, symmetric dipole configuration (↑↑) and the lower band corresponds to a dark, antisymmetric dipole configuration (↑↓). The light-cone (shaded region) is bounded by the linear dispersion of the TEM mode. Due to the non-trivial winding in the light–matter interaction (see Fig. 3), the band crossings are expected to result in large (band crossings ‘1’ and ‘2’) or small (band crossings ‘3’ and ‘4’) direction-dependent anticrossings in the polariton spectrum. c–e Polariton dispersion obtained from the polariton Hamiltonian H pol (solid black lines) and the two-band Hamiltonian \(\bar{H}_{\mathrm{mat}}\) (orange dashed lines), for c subcritical (L = 5a), d critical (L = L c = 1.75a), and e supercritical (L = a) cavity heights, respectively. While type-I CDPs with an isotropic Dirac cone (see inset of c) exist even in the quasistatic dispersion (see b), new type-II SDPs with a critically tilted Dirac cone (see inset in c) emerge due to the vanishing light–matter interaction for the dark quasistatic band along the Γ−K(K′) directions (see Fig. 3). At the critical cavity height L c , three type-II SDPs merge with the type-I CDP (see Fig. 5) resulting in a quadratic band-degeneracy at K(K′) (see inset in d). After criticality, the type-II SDPs annihilate one another and the massless Dirac cone re-emerges at the type-I CDPs (see inset in e) accompanied by an inversion of chirality (see Fig. 5). Plots obtained with parameters \(\omega_{\mathrm{K}0}^{\mathrm{ph}} = 2.5\omega_{0}\) and Ω = 0.01ω 0 Full size image

We diagonalize the matter Hamiltonian (Eq. (1)) as \(H_{\mathrm{mat}}= \sum_{\tau=\pm}\sum_{\mathbf{q}}\hbar\omega_{\mathbf{q}\tau}^{\mathrm{mat}\vphantom{\dagger}}\beta_{\mathbf{q}\tau}^{\dagger}\beta_{\mathbf{q}\tau}\) where the bosonic operators \(\beta_{\mathbf{q}\tau }^{\dagger} = \psi_{\mathbf{q}}^{\dagger} | \psi _{\mathbf{q}\tau} \rangle\) create quasistatic collective-dipole normal modes with dispersion \(\omega _{\mathbf{q}\tau}^{\mathrm{mat}} = \tilde{\omega}_{0} + \tau \tilde{\Omega} |f_{\mathbf{q}}^{\vphantom{\prime}}|\). Here, τ indexes the upper (τ = +1) and lower (τ = −1) bands and \(\psi_{\mathbf{q}}^{\dagger} = (a_{\mathbf{q}}^{\dagger} , b_{\mathbf{q}}^{\dagger} )\) is a spinor creation operator. The spinors \(|\psi _{{\mathbf{q}}\tau} \rangle = (1, \tau {\mathrm{e}}^{{\mathrm{i}}\varphi_{\mathbf{q}}})^{{\mathrm{T}}}/\sqrt{2}\), where T denotes the transpose, describe an emergent pseudospin degree of freedom where the two components encode the relative amplitude and phase of the dipolar oscillations on the two inequivalent A and B sublattices, respectively, with φ q = arg(f q ). These spinors can be represented by a pseudospin vector on the Bloch sphere which reads S q τ = τ(cosφ q , sinφ q , 0).

At the high symmetry K and K′ points (see Fig. 2a), the sublattices decouple with no well-defined relative phase (i.e., f q = 0), giving rise to two inequivalent CDPs located at \(\pm {\mathbf{K}} = \pm ( \frac{4\uppi}{3\sqrt 3 a} , 0 )\) as observed in Fig. 2b. These CDPs correspond to vortices in the pseudospin vector field S q τ , which give rise to topological singularities in the Berry curvature47. Therefore, the CDPs are sources of quantized Berry \(w\uppi\), where w = ±1 is the topological charge of the Dirac point corresponding to the winding number of the vortex. As expected from the symmetry of the metasurface, the existence of the CDPs is robust against long-range Coulomb interactions as shown in Supplementary Note 1. In fact, for small cavity heights, the image dipoles quench long-range Coulomb interactions and the nearest-neighbor approximation becomes increasingly accurate as shown in Supplementary Figure 1.

To quadratic order in k = q−K (\(ka \ll 1\)), the effective matter Hamiltonian near the K point is \(H_{\mathrm{K}}^{\mathrm{eff}} = \sum_{\mathbf{k}} \psi _{\mathbf{k}}^{\dagger} {\mathcal{H}}_{\mathrm{K},\mathbf{k}}^{\mathrm{eff}}\psi _{\mathbf{k}}\), with spinor creation operator \(\psi_{{k}}^{\dagger} = \left(a_{\mathbf{k}}^{\dagger} ,b_{\mathbf{k}}^{\dagger} \right)\) and Bloch Hamiltonian

$${\cal H}_{{\mathrm{K}},{\mathbf{k}}}^{{\mathrm{eff}}} = \hbar \tilde \omega _0{\mathbb1}_2 - \hbar \tilde v{\mathbf{\sigma }} \cdot {\mathbf{k}} + \hbar \tilde t\left( {{\mathbf{\sigma }}^ \ast \cdot {\mathbf{k}}} \right)^{ \circ 2}$$ (2)

Here, \({\mathbb{1}}_2\) is the 2 × 2 identity matrix, σ = (σ x , σ y ) and σ* = (σ x , −σ y ) are vectors of Pauli matrices, and °2 represents the Hadamard (element-wise) square. Note that the image dipoles do not qualitatively affect the physics, but simply lead to a renormalization of the group velocity \(\tilde{v} = 3\tilde{\Omega} a/2\) and trigonal warping parameter \(\tilde{t} = 3\tilde{\Omega} a^2 / 8\). Apart from a global energy shift, Eq. (2) is equivalent to a 2D massless Dirac Hamiltonian to leading order in k, with an isotropic Dirac cone spectrum \(\omega_{\mathbf{k}\tau}^{\mathrm{mat}} = \tilde{\omega}_{0} + \tau \tilde{v} \left| {\mathbf{k}} \right|\) that is characterized by closed isofrequency contours. Therefore, as expected from the honeycomb symmetry, the CDP belongs to the type-I class of 2D Dirac points, and the corresponding spinors \(| \psi _{\mathbf{k}\tau} \rangle = (1, - \tau {\mathrm{e}}^{\mathrm{i}\theta _{\mathbf{k}}} )^{\mathrm{T}}/\sqrt 2\), where θ k = arctan(k y /k x ), represent massless Dirac collective-dipoles with a topological Berry phase of π. The effective Hamiltonian near the K′ point is given by \({\mathcal{H}}_{\mathrm{K}{\prime},{\mathbf{k}}}^{\mathrm{eff}} = ( {\mathcal{H}}_{\mathrm{K},-\mathbf{k}}^{\mathrm{eff}} )^{\ast}\), where the corresponding massless Dirac collective-dipoles have a topological Berry phase of −π, as required by time-reversal symmetry.

Hybridization with the photonic environment

Given the subwavelength nearest-neighbor separation, it is tempting to assert that the near-field Coulomb interactions in H mat capture the essential physics. In fact, we will show that this quasistatic description misses the profound influence of the surrounding photonic environment, which has a remarkably non-trivial effect on the Berry curvature and, therefore, on the corresponding nature of the polaritons.

Crucially, the metallic cavity supports a fundamental transverse electromagnetic (TEM) mode whose polarization (parallel to the dipole moments) and linear dispersion (see Fig. 2b) are independent of the cavity height. For brevity, in what follows we do not present the contributions from the other cavity modes since the essential physics emerges from the interaction with the fundamental TEM mode (see Methods for the full expressions). In fact, the higher order cavity modes become increasingly negligible for smaller cavities as they are progressively detuned from the dipole resonances.

The effects of the photonic environment are encoded in the free photonic Hamiltonian

$$H_{\mathrm{ph}}= \hbar \sum_{{\mathbf{q}}n} \omega _{{\mathbf{q}}n}^{\mathrm{ph}}c_{{\mathbf{q}}n}^{\dagger} c_{{\mathbf{q}}n}$$ (3)

and in the light–matter interaction Hamiltonian \(H_{\mathrm{int}}^{\vphantom{()}} = H_{\mathrm{int}}^{(1)} + H_{\mathrm{int}}^{(2)}\), with

$$\begin{array}{*{20}{l}} {H_{{\mathrm{int}}}^{(1)}} \hfill & = \hfill & {\hbar \mathop {\sum}\limits_{{\mathbf{q}}n} {\mathrm{i}}\xi _{{\mathbf{q}}n}\phi _n^ \ast \left( {a_{\mathbf{q}}^\dagger c_{{\mathbf{q}}n} + a_{\mathbf{q}}^\dagger c_{ - {\mathbf{q}}n}^\dagger } \right)} \hfill \\ {} \hfill & {} \hfill & { + \hbar \mathop {\sum}\limits_{{\mathbf{q}}n} {\mathrm{i}}\xi _{{\mathbf{q}}n}\phi _n\left( {b_{\mathbf{q}}^\dagger c_{{\mathbf{q}}n} + b_{\mathbf{q}}^\dagger c_{ - {\mathbf{q}}n}^\dagger } \right) + {\mathrm{H}}{\mathrm{.c}}{\mathrm{.}}} \hfill \end{array}$$ (4)

and

$$\begin{array}{*{20}{l}} {H_{{\mathrm{int}}}^{(2)}} \hfill & = & { \hbar \mathop {\sum}\limits_{{\mathbf{q}}nn\prime } \frac{{2\xi _{{\mathbf{q}}n}\xi _{{\mathbf{q}}n\prime }}}{{\omega _0}}{\rm Re}\left( {\phi _n\phi _{n\prime }^ \ast } \right)} \hfill \\ {} & {} & {\left( {c_{{\mathbf{q}}n}^\dagger c_{{\mathbf{q}}n\prime } + c_{{\mathbf{q}}n}^\dagger c_{ - {\mathbf{q}}n\prime }^\dagger } \right) + {\mathrm{H}}{\mathrm{.c}}.} \hfill \end{array}$$ (5)

where \(\xi _{{\mathbf{q}}n} \propto L^{ -1/2}\) parametrizes the strength of the light–matter interaction (see Methods for analytical expression). The bosonic operator \(c_{{\mathbf{q}}n}^{\dagger}\) creates a TEM photon with wavevector q in the first Brillouin zone and dispersion \(\omega_{{\mathbf{q}}n}^{\mathrm{ph}} = c\left| {\mathbf{q}} - {\mathbf{G}}_{n}\right|\), where n indexes the set of reciprocal lattice vectors G n . The complex phase factors \(\phi_{n} = \exp \left( {\mathrm{i}}a{\mathbf{G}}_{n} \cdot \hat{\mathbf{y}} \right)\) are associated with Umklapp processes that arise due to the discrete, in-plane translational symmetry of the metasurface, and must be retained as they are critical for maintaining the point-group symmetry of the polariton Hamiltonian.

We diagonalize H pol using a generalized Hopfield–Bogoliubov transformation48 (see Methods for details), and in Fig. 2c–e, we present the resulting polariton dispersion for different cavity heights. Also, in Supplementary Figure 2, we present the full polariton dispersion which includes long-range Coulomb interactions. For small cavity heights, the full polariton dispersion is almost indistinguishable from that obtained in the nearest-neighbor approximation, and therefore one can conclude that long-range Coulomb interactions do not qualitatively affect the physics presented here. It is important to stress that our general model captures the essential physics that will emerge in a variety of different experimental setups. To show this, in Supplementary Figure 3 and Supplementary Figure 4 we present the polariton dispersions obtained from finite element simulations of a honeycomb array of plasmonic nanorods and microwave helical resonators, respectively. Indeed, these entirely different physical realizations show the same evolution of the polariton spectrum as presented in Fig. 2c–e and Supplementary Figure 2.

Emergence of type-II Dirac points

Given the elementary nature of the individual resonant elements, one may be tempted to assume that nothing peculiar could emerge from the ordinary dipole–dipole interactions between the meta-atoms which are mediated by the electromagnetic field. However, by expressing the interaction Hamiltonian (Eq. (4)) in terms of the β q τ and \(\beta_{\mathbf{q}\tau}^{\dagger}\) operators that diagonalize the matter Hamiltonian,

$$H_{\mathrm{int}}^{(1)} = \hbar \sum_{\tau=\pm} \sum_{{\mathbf{q}}n} \mathrm{i}\Lambda_{\mathbf{q}\it n\tau} \left( \beta _{\mathbf{q}\tau }^{\dagger} c_{{\mathbf{q}}\it n} + \beta_{\mathbf{q}\tau }^{\dagger} c_{ -{\mathbf{q}}\it n}^{\dagger} \right) + \mathrm{H.c.}\,,$$ (6)

we find that complex non-local interactions, which arise from strong multiple-scattering in the bipartite structure, result in a non-trivial winding of the light–matter coupling as a function of wavevector direction

$$\Lambda_{{\mathbf{q}}n\tau } \propto \xi_{{\mathbf{q}}n}\left( \phi_{n}^{\ast\vphantom{\dagger}} {\mathrm{e}}^{\mathrm{i}\varphi_{\mathbf{q}}} + \tau\phi_{n} \right)\,.$$ (7)

Naively, one may expect all of the band crossings in Fig. 2b to be avoided as a result of the hybridization between the collective-dipole and photonic modes, as it is a characteristic feature of polaritonic systems48,49. Indeed, this is the case for the crossings with the upper quasistatic band where \(\Lambda_{\mathbf{q}\mathrm 0+} \propto \left( {\mathrm{e}}^{\mathrm{i}\varphi_{\mathbf{q}}} + 1 \right)\) (see red line in Fig. 3a) due to the constructive interference between the sublattices of this bright (↑↑) configuration (see Fig. 3b, c). This results in a large anticrossing for all wavevector directions, as observed in Fig. 2c. In stark contrast, for the lower quasistatic band the coupling constant is significantly reduced \(\Lambda_{\mathbf{q}\mathrm 0-} \propto \left( {\mathrm{e}}^{\mathrm{i}\varphi_{\mathbf{q}}} - 1 \right)\) (see blue line in Fig. 3a) due to the destructive interference between the sublattices of this dark (↑↓) configuration (see Fig. 3e). Consequently, this results in a small anticrossing for a general wavevector direction.

Fig. 3 Non-trivial winding in the light–matter interaction. a Dependence of the magnitude of the light–matter coupling constant \(| \Lambda_{\mathbf{q}\mathrm 0\tau} | \propto |{\mathrm{e}}^{\mathrm{i}\varphi _{\mathbf{q}}} + \tau |\) on the direction of q from Γ−M to Γ−K(K′) (see inset), for the upper (red line) and lower (blue line) quasistatic bands. Plots obtained with |q| = |K|/2. b–e Schematics of the bright (↑↑) and dark (↑↓) configurations of the two sublattices interacting with the photonic mode which is indicated by the field profile. Panels b and d represent the crossings labeled ‘1’ and ‘3’ in Fig. 2 along the Γ−K(K′) directions, respectively, while c and e represent the crossings labeled ‘2’ and ‘4’ along the Γ−M directions, respectively. Crucially, the light–matter interaction strength for the dark mode vanishes (|Λ q0 − | = 0) along the Γ−K(K′) directions due to the complete destructive interference between the two sublattices (see d), leading to the emergence of six inequivalent type-II Dirac points in the polariton spectrum Full size image

Crucially, however, the light–matter interaction for the lower quasistatic band completely vanishes (Λ q0− = 0) along the high-symmetry Γ−K(K′) directions, where φ q = 0, due to the complete destructive interference between the two sublattices (see Fig. 3d). As a result, along these high-symmetry directions the crossings are protected, leading to six inequivalent Dirac points emerging in the polariton spectrum—we call these satellite Dirac points (SDPs) to distinguish them from the CDPs. As we will see below, these SDPs belong to the type-II class of 2D Dirac points where the dispersion takes the form of a critically tilted Dirac cone (see inset of Fig. 2c), characterized by open, hyperbolic isofrequency contours.

Effective Hamiltonian in the matter subspace

To explore the nature of the polaritons in the vicinity of the different Dirac points, we first neglect non-resonant terms in the matter Hamiltonian and perform a unitary Schrieffer–Wolff transformation50 on H pol to integrate out the photonic degrees of freedom (see Methods for details). Finally, we extract the two-band Hamiltonian in the matter sublattice space

$$\begin{array}{*{20}{l}}\bar{H}_{\mathrm{mat}}=\,\,\,&H_{\mathrm{mat}}-2\hbar\sum_{{\mathbf{q}}n}\frac{\xi_{{\mathbf{q}}n}^{2\vphantom{\dagger}} \omega_{{\mathbf{q}}n}^{\mathrm{ph}} }{\left(\omega_{{\mathbf{q}}n}^{\text{ph}}\right)^2-\tilde{\omega}_{0}^{2}} \\ & \left(a_{\mathbf{q}}^{\dagger} a_{\mathbf{q}}+b_{\mathbf{q}}^{\dagger} b_{\mathbf{q}}+\phi_{n}^{2\vphantom{\dagger}}b_{\mathbf{q}}^{\dagger} a_{\mathbf{q}}+\phi_{n}^{*2\vphantom{\dagger}}a_{\mathbf{q}}^{\dagger}b_{\mathbf{q}} \right)\,.\end{array}$$ (8)

Diagonalizing the two-band Hamiltonian (Eq. (8)) leads to an effective dispersion (see Methods) which provides an excellent description of the polaritons as indicated by the orange dashed lines in Fig. 2c–e. Finally, we expand the two-band Hamiltonian (Eq. (8)) up to quadratic order in k and obtain the effective Hamiltonian near the K point \(\bar{H}_{\mathrm{K}}^{\mathrm{eff}} = \sum_{\mathbf{k}} \psi _{\mathbf{k}}^{\dagger} \bar{\mathcal{H}}_{\mathrm{K},\mathbf{k}}^{\mathrm{eff}}\psi _{\mathbf{k}}\) (see Supplementary Note 2 for derivation) with Bloch Hamiltonian

$$\bar{\mathcal{H}}_{\mathrm{K},\mathbf{k}}^{\mathrm{eff}} = \hbar \bar{\omega}_{0} {\mathbb1}_{2} - \hbar \bar{v}{\boldsymbol{\sigma }} \cdot {\mathbf{k}} + \hbar \bar{t} \left( {\boldsymbol{\sigma }}^{\ast} \cdot {\mathbf{k}} \right)^{\circ 2} - \hbar D \left| {\mathbf{k}} \right|^{2} {\mathbb1}_{2}\,.$$ (9)

Similarly, the effective Hamiltonian near the K′ point is given by \(\bar{\mathcal{H}}_{\mathrm{K}{\prime} ,\mathbf{k}}^{\mathrm{eff}} = (\bar{\mathcal{H}}_{\mathrm{K}, - {\mathbf{k}}}^{\mathrm{eff}} )^{\ast}\). In Eq. (9), the resonant frequency \(\bar{\omega}_{0}\), group velocity \(\bar{v}\), and trigonal warping parameter \(\bar{t}\), now encode non-trivial contributions from the hybridization with the photonic environment. There is also an additional wavevector-dependent diagonal term parametrized by D, which breaks the symmetry between the upper and lower polariton bands. The dependence of these parameters on the cavity height is shown in Fig. 4 (see Methods for analytical expressions). To leading order in k, one can observe that the effective Hamiltonian (Eq. (9)) near the CDP is equivalent to a 2D massless Dirac Hamiltonian. Therefore, the polariton CDPs remain in the type-I class and are robust against the coupling with the photonic environment—this is not surprising given that their physical origin is intrinsically linked to the lattice structure alone, which is preserved here.

Fig. 4 Tunable parameters in the effective Hamiltonian. Dependence of the parameters in the effective polariton Hamiltonian (Eq. (9)) on the inverse cavity height. The blue dashed line shows the variation of the group velocity \(\bar{v}\) which changes sign at the critical cavity height L c , leading to the inversion of chirality. The orange dot-dashed line shows the variation of the trigonal warping parameter \(\bar{t}\) which becomes dominant close to criticality. These parameters have been normalized to v = 3Ωa/2 and t = 3Ωa2/8 which are the group velocity and trigonal warping parameters, respectively, in the absence of image dipoles and light–matter interactions. The orange dot-dashed line in the inset shows the variation of the CDP frequency \(\bar{\omega}_{0}\), while the blue dashed line in the inset shows the variation of the wavevector-dependent diagonal term D. Plots obtained with \(\omega_{\mathbf{K}\mathrm 0}^{\mathrm{ph}} = 2.5\omega_{0}\) and Ω = 0.01ω 0 Full size image

To elucidate the nature of the SDPs, we expand the effective Hamiltonian (Eq. (9)) near one of the SDPs located at \({\mathbf{K}}_{\mathrm{S}} = \left(\bar{v}/\bar{t},0 \right)\) and obtain

$$\bar{\mathcal{H}}_{{\mathrm{K}}_{\mathrm{S}},{\mathbf{k}}{\prime}}^{\mathrm{eff}} = \hbar \left( \bar{\omega}_{0}- \frac{D\bar{v}^2}{\bar{t}^2} - \frac{2D\bar{v}}{\bar{t}}k_x\kern-4pt {\prime} \right) {\mathbb1}_{2} + \hbar {\boldsymbol{\sigma }}^{\ast} \cdot \bar{\mathbf{v}} \cdot {\mathbf{k}}{\prime}\,,$$ (10)

where k′ measures the deviation from K S and \(\overline {\mathbf{v}} = \bar v\left( {\begin{array}{*{20}{c}} 1 & 0 \\ 0 & 3 \end{array}} \right)\) is the velocity tensor. Apart from a global energy shift, the effective Hamiltonian (Eq. (10)) near the SDP takes the form of a generalized 2D massless Dirac Hamiltonian \({\cal H}_{\mathbf{k}} = \mathop {\sum}

olimits_{i = x,y} \hbar u_ik_i{\mathbb1}_2 + \mathop {\sum}

olimits_{i = x,y} \hbar v_ik_i\sigma _i\). If the parameters u i and v i satisfy the condition \(u_x^2 / v_x^2 + u_y^2/ v_y^2 < 1\), then the Dirac cone becomes tilted and anisotropic51 but still belongs to the type-I class with closed isofrequency contours. However, the condition \(u_x^2/v_x^2 + u_y^2/v_y^2 > 1\) defines a distinct type-II class of 2D Dirac points, giving rise to a critically tilted Dirac cone with open, hyperbolic isofrequency contours. Hence, the type-I and type-II classes are related via a Lifshitz transition in the topology of the isofrequency contours. Indeed, since we have u y = 0 and \(u_x^2/v_x^2 = 4D^2/\bar{t}^2 > 1\), the SDPs belong to the type-II class of 2D Dirac points. Furthermore, since the Hamiltonian (Eq. (10)) is expressed in terms of σ*, the pseudospin winds in the opposite direction around the SDPs as compared to the CDP, and therefore the SDPs located along the Γ−K directions are sources of −π Berry flux. As required by time-reversal symmetry, the SDPs located along the Γ−K′ directions are sources of π Berry flux (opposite to the CDP located at the K′ point).

Manipulation of type-I and type-II Dirac points

We have thus demonstrated that the honeycomb metasurface simultaneously exhibits two distinct species of massless Dirac polaritons, namely type-I and type-II. In contrast to the type-I CDPs, the existence of the type-II SDPs is intrinsically linked to the hybridization between the light and matter degrees of freedom, and thus one can manipulate their location within the Brillouin zone by simply modifying the light–matter interaction via the cavity height. As a result, the polariton spectrum evolves into qualitatively distinct phases as highlighted in Fig. 2c–e. To elucidate the differences between these phases, we study the spinor eigenstates (see Methods) of the two-band Hamiltonian (Eq. (8)). In Fig. 5a–c we plot the pseudospin vector field near the K point for different cavity heights and schematically depict the location of the Dirac points, along with their associated Berry flux. Finally, in Fig. 5d–f, we illustrate the corresponding effective polariton spectrum to leading order in k. Note that similar analysis can be performed near the K′ point.

Fig. 5 Merging of type-I and type-II Dirac points with chirality inversion. a–c Pseudospin vector field and isofrequency contours (see Methods for the specific isofrequency values) for the upper polariton band near the K point for a subcritical (L = 2.3a), b critical (L = L c = 1.78a), and c supercritical (L = 1.4a) cavity heights, respectively, as obtained from \(\bar{H}_{\mathrm{mat}}\). The Dirac points (corresponding to vortices) are depicted by orange circles along with their associated Berry flux. Before criticality, three type-II SDPs (−π Berry flux) are driven towards the type-I CDP (π Berry flux) along the Γ−K directions as the cavity height is reduced (see inset in a). At the critical cavity height L c , they merge together forming a quadratic band-degeneracy with combined Berry flux of −2π (see b). After criticality, the type-II SDPs re-emerge and are driven past the type-I CDP along the K−M directions (see inset in c). After a small decrease in cavity height, these SDPs annihilate other SDPs that migrate along the opposite direction and have opposite Berry flux, leaving only the type-I CDP remaining in the spectrum with π Berry flux (see c). d-f, Effective polariton spectrum near the K point to leading order in k. The colors of the bands correspond to the chirality of the Dirac polaritons as defined in the main text, where the orange and blue bands indicate a chirality of +1 and −1, respectively. The spinor eigenstates, represented by pseudospin vectors (gold arrows), describe d massless Dirac polaritons with a linear dispersion and Berry phase of π, e massive chiral polaritons with a parabolic dispersion and Berry phase of −2π, and f massless Dirac polaritons with a linear dispersion and Berry phase of π, but with inverted chirality. All pseudospin and contour plots are obtained with \(\omega_{\mathbf{K}\mathrm 0}^{\mathrm{ph}} = 2.5\omega_{0}\) and Ω = 0.01ω 0 Full size image

In the subcritical phase (L > L c ), three type-II SDPs are located along the Γ−K directions, each with −π Berry flux surrounding a type-I CDP with π Berry flux (see Fig. 5a). To leading order in k, the polariton spectrum disperses linearly about the type-I CDPs (see Fig. 2c) forming an isotropic Dirac cone with a group velocity \(\bar{v}\) that is tunable via the cavity height (see Fig. 4). Here, the effective Hamiltonian (Eq. (9)) is equivalent to a 2D massless Dirac Hamiltonian with spinor eigenstates \(| \psi_{\mathbf{k}\tau} \rangle = (1, -\tau {\mathrm{e}}^{\mathrm{i}\theta _{\mathbf{k}} } )^{\mathrm{T}} /\sqrt{2}\). These represent massless Dirac polaritons with chirality \(\langle \psi_{\mathbf{k}\tau } | {\boldsymbol{\sigma}} \cdot \hat{\mathbf{k}} | \psi _{\mathbf{k}\tau }\rangle = - \tau\), resulting in a pseudospin that winds once around the CDP and a topological Berry phase of π (see Fig. 5d).

At the critical cavity height (L = L c ), the group velocity of the massless Dirac polaritons vanishes \(\bar{v}(L_{\mathrm{c}}) = 0\) (see Fig. 4) as the type-II SDPs merge with the type-I CDP, forming a quadratic band-degeneracy (see Fig. 2d) with combined −2π Berry flux (see Fig. 5b). The leading order term in the effective Hamiltonian (Eq. (9)) is now quadratic in k with corresponding spinor eigenstates \(| \psi_{\mathbf{k}\tau} \rangle = (1, - \tau {\mathrm{e}}^{ - \mathrm{i}2\theta _{\mathbf{k}}} )^{\mathrm{T}}/\sqrt 2\). Therefore, during this critical merging transition, the massless Dirac polaritons morph into massive chiral polaritons with qualitatively distinct physical properties. These include a parabolic spectrum and chirality \(\langle \psi _{\mathbf{k}\tau} | ( {\boldsymbol{\sigma}}^{\ast} \cdot \hat{\mathbf{k}} )^{\circ 2} |\psi _{\mathbf{k}\tau} \rangle = - \tau\), resulting in a pseudospin that winds twice as fast compared to the subcritical phase and a topological Berry phase of −2π (see Fig. 5e).

Since the point-group symmetry is preserved, the type-II SDPs do not annihilate the type-I CDP, but they re-emerge and continue to migrate along the K−M directions as the cavity height is reduced past criticality (L < L c ) (see inset of Fig. 5c). After a small decrease in cavity height, these SDPs annihilate with other SDPs that migrate along the opposite direction and have opposite Berry flux. This topological transition leaves only the type-I CDP remaining in the polariton spectrum with π Berry flux (see Figs. 2e and 5c).

In this supercritical phase, we recover the linear dispersion near the type-I CDP to leading order in k (see Fig. 2e), and the effective Hamiltonian (Eq. (9)) is equivalent to a 2D massless Dirac Hamiltonian with corresponding spinor eigenstates \(| \psi _{\mathbf{k}\tau} \rangle = (1,\tau {\mathrm{e}}^{\mathrm{i}\theta _{\mathbf{k}}} )^{\mathrm{T}}/\sqrt 2\). Remarkably, massless Dirac polaritons thus re-emerge past criticality with an environment-induced inversion of chirality \(\langle \psi _{\mathbf{k}\tau} | {\boldsymbol{\sigma }} \cdot \hat{\mathbf{k}} | \psi _{\mathbf{k}\tau} \rangle = \tau\) (see Fig. 5f). Physically, this corresponds to a π-rotation in the relative phase between the dipole oscillations on the two inequivalent sublattices, which is also accompanied by a π-rotation in the isofrequency domains (compare Fig. 5a and Fig. 5c).

We emphasize that it is the chirality of massless Dirac fermions that is responsible for most of the remarkable properties of monolayer graphene, including the Klein tunneling phenomenon13. Consequently, this unique environment-induced inversion of chirality could give rise to unconventional phenomena such as anomalous Klein transport. For example, near the K point, the right-propagating polaritons correspond to an antisymmetric dipole configuration \(| \psi _{\mathbf{k}\tau} \rangle = (1, - 1 )^{\mathrm{T}}/\sqrt 2\) in the subcritical phase and to a symmetric configuration \(| \psi _{\mathbf{k}\tau}\rangle = (1,1)^{\mathrm{T}}/\sqrt 2\) in the supercritical phase. Thus, due to the orthogonality between these two spinor eigenstates, the inversion of chirality removes the channel responsible for the perfect transmission in the conventional Klein tunneling effect13 (see Fig. 5d, f). Such a scenario could be realized in a simple setup characterized by two neighboring regions with different cavity heights.

As a side remark, we note that the polariton spectrum near criticality bears some resemblance with the low-energy spectrum of bilayer graphene with its central Dirac point and three satellite Dirac points, which all belong to the type-I class52,53. However, given the type-II nature of the polariton SDPs, the topology of the polariton isofrequency contours are markedly different from that of the bilayer spectrum. This is further highlighted at criticality where the polariton bands have the same curvature, which is in stark contrast to the electronic bands in bilayer graphene.

We also note that recent works explored the possibility to manipulate the (3+1) type-I Dirac points in bilayer graphene through the application of lattice deformations54,55,56,57, leading to the merging and annihilation of pairs of Dirac points. In addition, a multimerging transition of all (3+1) type-I Dirac points has been proposed theoretically within tight-binding models involving the artificial tuning of third-nearest-neighbor hopping amplitudes in a graphene-like honeycomb structure58,59,60,61. However, these proposals have no physical realization so far. In stark contrast, the imminently realizable metasurfaces in our work enable the exploration of rich Dirac phases with ease by simply modifying the photonic environment via an enclosing cavity.

As a final remark, we briefly comment on how one might probe the Dirac physics presented in this work. Given that the Dirac points exist in a polaritonic excitation spectrum, one must drive the system with photons at the required frequency in order to probe them. In fact, both of the type-I and type-II Dirac points lie outside of the light-line and therefore one must overcome the momentum mismatch with photons. The specific experimental technique that one would employ will depend on the nature of the metasurface and the corresponding frequency regime. For example, techniques for plasmonic systems have traditionally involved coupling via evanescent waves with prisms, gratings, and local scatterers62, or more recent techniques such as non-linear wave-mixing63. In contrast, realizations in the microwave regime can be probed using point-like antenna sources and detectors33. In fact, microwave metamaterials are proving to be a versatile platform for exploring Dirac/Weyl physics, as one can directly probe the field distributions using near-field scanning techniques33, and thus one could directly probe the environment-induced chirality inversion predicted here.