Manipulating type-I and type-II Dirac polaritons in cavity-embedded honeycomb metasurfaces

Pseudorelativistic Dirac quasiparticles have emerged in a plethora of artificial graphene systems that mimic the underlying honeycomb symmetry of graphene. However, it is notoriously difficult to manipulate their properties without modifying the lattice structure. Here we theoretically investigate polaritons supported by honeycomb metasurfaces and, despite the trivial nature of the resonant elements, we unveil rich Dirac physics stemming from a non-trivial winding in the light–matter interaction. The metasurfaces simultaneously exhibit two distinct species of massless Dirac polaritons, namely type-I and type-II. By modifying only the photonic environment via an enclosing cavity, one can manipulate the location of the type-II Dirac points, leading to qualitatively different polariton phases. This enables one to alter the fundamental properties of the emergent Dirac polaritons while preserving the lattice structure—a unique scenario which has no analog in real or artificial graphene systems. Exploiting the photonic environment will thus give rise to unexplored Dirac physics at the subwavelength scale.

T he groundbreaking discovery of monolayer graphene 1 has inspired an extensive quest to emulate massless Dirac quasiparticles in a myriad of distinct artificial graphene systems [2][3][4][5][6][7][8][9][10][11] , ranging from ultracold atoms in optical lattices 3 to evanescently coupled photonic waveguide arrays 4 . Owing to their honeycomb symmetry, linear band-degeneracies manifest in the quasiparticle spectrum which we call conventional Dirac points (CDPs). These belong to the ubiquitous type-I class of twodimensional (2D) Dirac points that are characterized by Dirac cones with closed isofrequency contours. As a result, the corresponding quasiparticles are described by the rather exotic 2D massless Dirac Hamiltonian 12 , and thus offer fundamental insight into pseudorelativistic phenomena such as the iconic Klein paradox 13 . The latter is responsible for the suppression of backscattering and for the antilocalization of massless Dirac quasiparticles, which are highly desirable properties for efficient quasiparticle propagation in novel devices.
Since the existence of type-I CDPs is intrinsically linked to the honeycomb structure, the fundamental properties of the massless Dirac quasiparticles are notoriously robust and difficult to manipulate. However, by exploiting meticulous control over the lattice structure, artificial graphene systems have enabled the exploration of Dirac quasiparticles in new regimes that are difficult, if not impossible to achieve in graphene itself [14][15][16][17][18][19] . Among others, the archetypal example which has attracted considerable interest is the paradigm of strain-engineering, where it has been shown that lattice anisotropy can induce the merging and annihilation of type-I CDPs 3, [14][15][16][20][21][22][23] , and that aperiodicity can generate large pseudomagnetic fields 17,24 .
Moreover, the recent discovery of type-II Dirac/Weyl semimetals [25][26][27][28][29] sparked a burgeoning exploration into the prospects of a rarer type-II class of three-dimensional Dirac/Weyl points. As the latter are characterized by critically tilted Dirac/Weyl cones with open, hyperbolic isofrequency contours, the corresponding Lorentz-violating Dirac/Weyl quasiparticles exhibit markedly different properties from their type-I counterparts [25][26][27][28][29] . Soon after their realization, electromagnetic analogs emerged [30][31][32][33][34] , and this exploration has recently been extended to 2D systems where a distinct type-II class of 2D Dirac points were theoretically predicted 35,36 . However, since their existence is predicated on strong anisotropy in judiciously engineered photonic structures, one cannot manipulate their properties without modifying the lattice structure.
This hunt for exotic quasiparticles has recently entered the realm of polaritonics [37][38][39][40][41][42] . The true potential of polaritons lies in their hybrid nature, where their light and matter constituents can be manipulated independently, thereby providing additional tunable degrees of freedom. Among other examples, recent works have shown the tantalizing prospect of engineering novel topological polaritons by introducing a winding coupling between ordinary photons and excitons 39,41 .
In this work, we exploit the hybrid nature of polaritons in a different setting, namely metasurfaces, and we unveil unique Dirac physics by shifting the focus from the lattice structure and its deformations to the effect of manipulating the surrounding photonic environment. In particular, we theoretically study the polaritons supported by imminently realizable, crystalline metasurfaces consisting of a honeycomb array of resonant, dipolar meta-atoms. Despite the elementary nature of these metasurfaces, we unveil the simultaneous existence of both type-I and type-II massless Dirac polaritons which have distinct physical origins. Crucially, the existence of the latter is not a result of anisotropy but is intrinsically linked to the hybrid nature of the polaritons, emerging from a non-trivial winding in the light-matter interaction. Furthermore, we show that by embedding the honeycomb metasurface inside a planar photonic cavity and simply changing the cavity height, one can induce multiple phase transitions including the multimerging of type-I and type-II Dirac points and the annihilation of type-II Dirac points. This striking tunability results in qualitatively different polariton phases, despite the preserved lattice structure. In particular, we unveil a morphing between a linear and a parabolic spectrum accompanied by a change in the topological Berry phase, and an environment-induced inversion of chirality, all of which have no analog in graphene or artificial graphene systems studied thus far. Therefore, this unique paradigm will give rise to unexplored Dirac-related phenomena at the subwavelength scale, such as anomalous Klein tunneling, negative refraction, and pseudomagnetic Landau levels, which can all be tuned via the photonic environment alone.

Results
Hamiltonian formulation. While metamaterials have traditionally been described in terms of macroscopic effective properties 30,33,43 , the importance of crystallinity is becoming increasingly apparent 44 . Therefore, to capture the essential physics related to complex non-local effects that arise from strong multiple-scattering 45 , 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-coupling 46 .
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).
Emergence of type-I Dirac points. The matter Hamiltonian within the nearest-neighbor approximation reads where, for brevity, we have not presented the non-resonant terms (see Methods for derivation). In Eq. (1),ω 0 is the renormalized resonant frequency andΩ 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 y q and b y q 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 q ¼ P 3 j¼1 expðiq Á e j Þ encodes the honeycomb geometry of the lattice with nearest-neighbor vectors e j (see Fig. 1).
We diagonalize the matter Hamiltonian (Eq. (1)) as H mat ¼ P τ¼ ± P q hω mat qτ β y qτ β qτ where the bosonic operators β y qτ ¼ ψ y q jψ qτ i create quasistatic collective-dipole normal modes with dispersion ω mat qτ ¼ω 0 þ τΩjf q j. Here, τ indexes the upper (τ = +1) and lower (τ = −1) bands and ψ y q ¼ ða y q ; b y q Þ is a spinor creation operator. The spinors jψ qτ i ¼ ð1; τe iφ q Þ T = ffiffi ffi 2 p , 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 ± K ¼ ± ð 4π 3 ffiffi 3 p 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 curvature 47 . Therefore, the CDPs are sources of quantized Berry wπ, 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 longrange Coulomb interactions as shown in Supplementary Note 1. In fact, for small cavity heights, the image dipoles quench longrange Coulomb interactions and the nearest-neighbor approximation becomes increasingly accurate as shown in Supplementary Figure 1.
To quadratic order in k = q−K (ka ( 1), the effective matter Hamiltonian near the K point is H eff K ¼ P k ψ y k H eff K;k ψ k , with spinor creation operator ψ y k ¼ a y k ; b y k and Bloch Hamiltonian Here, 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ṽ ¼ 3Ωa=2 and trigonal warping parametert ¼ 3Ω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 ω mat kτ ¼ω 0 þ τṽ k j j 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 jψ kτ i ¼ ð1; Àτe iθ k Þ T = ffiffi ffi 2 p , 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 H eff K′;k ¼ ðH eff K;Àk Þ Ã , 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 and in the light-matter interaction Hamiltonian where ξ qn / L À1=2 parametrizes the strength of the light-matter interaction (see Methods for analytical expression). The bosonic operator c y qn creates a TEM photon with wavevector q in the first Brillouin zone and dispersion ω ph qn ¼ c q À G n j j , where n indexes the   2 Þ, where a is the subwavelength nearest-neighbor separation. Each subwavelength metaatom 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 set of reciprocal lattice vectors G n . The complex phase factors ϕ n ¼ exp iaG n Áŷ ð Þ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 pointgroup symmetry of the polariton Hamiltonian.
We diagonalize H pol using a generalized Hopfield-Bogoliubov transformation 48 (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 nearestneighbor 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  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 β y qτ operators that diagonalize the matter Hamiltonian, 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 Naively, one may expect all of the band crossings in Fig. 2b to be avoided as a result of the hybridization between the collectivedipole and photonic modes, as it is a characteristic feature of polaritonic systems 48,49 . Indeed, this is the case for the crossings with the upper quasistatic band where Λ q0þ / e iφ q þ 1 À Á (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 Λ q0À / e iφ q À 1 À Á (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.
Crucially, however, the light-matter interaction for the lower quasistatic band completely vanishes (Λ q0− = 0) along the  Fig. 2 Evolution of the polariton dispersion as the cavity height is reduced. a First Brillouin zone defined by primitive reciprocal lattice vectors 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 H 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 ω 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 transformation 50 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 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 H Similarly, the effective Hamiltonian near the K′ point is given by H eff K′;k ¼ ðH eff K;Àk Þ Ã . In Eq. (9), the resonant frequency ω 0 , group velocity v, and trigonal warping parameter t, now encode nontrivial 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.
To elucidate the nature of the SDPs, we expand the effective Hamiltonian (Eq. (9)) near one of the SDPs located at where k′ measures the deviation from K S and v ¼ v 1 0 0 3 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 H k ¼ P i¼x;y hu i k i 1 2 þ P i¼x;y hv i k i σ i . If the parameters u i and v i satisfy the condition u 2 x =v 2 x þ u 2 y =v 2 y <1, then the Dirac cone becomes tilted and anisotropic 51 but still belongs to the type-I class with closed isofrequency contours. However, the condition u 2 x =v 2 x þ u 2 y =v 2 y >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 2 x =v 2 x ¼ 4D 2 =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.
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 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 jψ kτ i ¼ ð1; Àτe iθ k Þ T = ffiffi ffi 2 p . These represent massless Dirac polaritons with chirality hψ kτ jσ Ákjψ kτ i ¼ Àτ, 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 vðL 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 jψ kτ i ¼ ð1; Àτe Ài2θ k Þ T = ffiffi ffi 2 p . 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 hψ kτ jðσ Ã ÁkÞ 2 jψ kτ i ¼ Àτ, 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 jψ kτ i ¼ ð1; τe iθ k Þ T = ffiffi ffi 2 p . Remarkably, massless Dirac polaritons thus re-emerge past criticality with an environment-induced inversion of chirality hψ kτ jσ Ákjψ kτ i ¼ τ (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 phenomenon 13 . 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 jψ kτ i ¼ ð1; À1Þ T = ffiffi ffi 2 p in the subcritical phase and to a symmetric configuration jψ kτ i ¼ ð1; 1Þ T = ffiffi ffi 2 p 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 effect 13 (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 class 52,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 deformations [54][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 structure [58][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 scatterers 62 , or more recent techniques such as non-linear wave-mixing 63 . In contrast, realizations in the microwave regime can be probed using point-like antenna sources and detectors 33 . 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 techniques 33 , and thus one could directly probe the environment-induced chirality inversion predicted here.

Discussion
To conclude, we have revealed rich and unique Dirac physics that emerges even in the most elementary honeycomb metasurfaces. In particular, we have unveiled the simultaneous existence of both type-I and type-II massless Dirac polaritons, where the latter emerge from a non-trivial winding in the light-matter interaction. We would like to emphasize that it is this unique physical origin of the type-II SDPs, together with the truly 2D nature of the metasurface, that enables one to qualitatively modify the fundamental properties of these emergent Dirac polaritons by manipulating the surrounding photonic environment alone. This stands in stark contrast to conventional artificial graphene systems where the fundamental properties are dictated by the lattice structure. Therefore, exploiting the rich tunability of the polariton spectrum with the environment offers a new paradigm that opens a variety of opportunities to explore unique Dirac-related physics at the subwavelength scale.
For example, one can simultaneously probe the dynamics of type-I and type-II massless Dirac quasiparticles, where the latter are predicted to exhibit intriguing anomalous refraction behavior 34,35 . Furthermore, the environment-induced redshift of the CDP frequency ω 0 (see Fig. 4) will allow the investigation of polaritonic Klein tunneling through interfaces separating regions with different cavity heights. Consequently, negative refraction can be induced by simple variations in the cavity height, which could be exploited in novel schemes for guiding and manipulating light at the subwavelength scale, including polaritonic Veselago lensing 64,65 . Moreover, the tunable group velocity will enable the exploration of velocity barriers for the unprecedented guiding and localization of massless Dirac quasiparticles 66,67 , which is extremely difficult to achieve in real graphene. One could also combine the effects of the environment with inhomogeneous strain deformations, giving rise to unique pseudomagnetic-related effects, including the intriguing ability to induce a pseudo-Landau level spectrum for polaritons that can be qualitatively tuned via the cavity height. Finally, the ability to controllably invert the chirality of the massless Dirac polaritons opens new perspectives for anomalous pseudorelativistic transport through interfaces separating regions in distinct polaritonic phases.

Methods
Derivation of the polaritonic Hamiltonian. The cavity-embedded metasurface is composed of a honeycomb array of identical meta-atoms located at Here, R = l 1 a 1 + l 2 a 2 is an in-plane lattice translation vector with primitive vectors a 1 and a 2 (see Fig. 1) and integers l 1 and l 2 . Each meta-atom is modeled by a single dynamical degree of freedom h (with dimensions of length), where the electricdipole moment associated with its fundamental eigenmode is p ¼ ÀQhẑ, with effective charge Q. The Coulomb potential energy between two dipole moments p and p′ located at generic positions r and r′, respectively, is given by wheren ¼ r À r′ ð Þ = r À r′ j jand ε 0 is the vacuum permittivity. The presence of the perfectly conducting metallic plates, placed at z = 0 and z = L, modifies the boundary conditions on the scalar potential and, therefore, the Coulomb interaction between the meta-atoms. Using the method of images to ensure the vanishing of the scalar potential at the cavity walls 68 , we introduce an infinite series of image dipoles located outside the cavity at positions R s þ lLẑ, where s = A, B labels the two sublattices and l is a non-zero integer. Noting that the Coulomb potential energy between a real and image dipole is 1/2 of that given by Eq. (11) 69 , the matter Hamiltonian within the nearest-neighbor approximation reads where the primed summations exclude the l = 0 term. Here, Π R s is the conjugate momentum to the dynamical coordinate h R s corresponding to the meta-atom located at R s , and M is an effective mass. Next, we introduce the bosonic operators that annihilate quanta of the fundamental eigenmode on the meta-atom located at R A and R B , respectively, and satisfy the commutation relations ½a R ; a y R′ ¼ δ RR′ , ½b R ; b y R′ ¼ δ RR′ , and ½a R ; b y R′ ¼ 0. In terms of these operators, the matter Hamiltonian (Eq. (12)) reads where Ω = Q 2 /8πε 0 Mω 0 a 3 parametrizes the strength of the nearest-neighbor Coulomb interaction, and the parameters encode renormalizations due to the cavity-induced image dipoles. We apply Bornvon Kármán boundary conditions over a lattice with N ) 1 unit cells and introduce the Fourier transform of the bosonic operators a R A ¼ N À1=2 P q a q e iqÁR A and b R B ¼ N À1=2 P q b q e iqÁR B , which transforms the matter Hamiltonian (Eq. (15)) into the local and block-diagonal form Ωð1 À IÞ f q b y q a q þ a y Àq þ H:c: In the main text, we do not present the non-resonant terms (e.g., b y q a y Àq ), leading to Eq. (1) whereω 0 ¼ ω 0 À ΩS andΩ ¼ Ω 1 À I ð Þ . In the Coulomb gauge, the light-matter interaction is described by the minimalcoupling Hamiltonian 46 which, within the dipole approximation, reads where we have used Π R ¼ Π Rẑ . The vector potential can be decomposed into transverse electric (TE) and transverse magnetic (TM) modes of the cavity. However, the photons corresponding to the TE modes have an in-plane polarization, and therefore only TM modes contribute to the z-component of the vector potential c qmn e iðqÀG n ÞÁr þ c y qmn e ÀiðqÀG n ÞÁr is the area of a unit cell and N m = 1 + δ m0 . The bosonic operator c y qmn creates a TM photon with wavevector q in the first Brillouin zone and dispersion ω ph qmn ¼ c q À G n þẑmπ=L j j . Here, G n = n 1 b 1 + n 2 b 2 is a reciprocal lattice vector with primitive vectors b 1 and b 2 , where n indexes the set of ordered pairs of integers (n 1 , n 2 ), and m is a non-negative integer indexing the different TM cavity modes. Only TM photons with even m couple to the dipoles due to the parity selection rule at the center of the cavity.
Substituting the vector potential (Eq. (19)) into Eq. (18) we obtain the light-matter interaction Hamiltonian given by The strength of the light-matter interaction is parametrized by where, to take into account the finite size of the meta-atoms, we have introduced a phenomenological function F ðω ph qmn Þ that provides a smooth cut-off for the interaction with short-wavelength photonic modes where the dipole approximation breaks down. We choose the phenomenological cut-off function to be of the Fermi-Dirac distribution form which is smooth enough to avoid spurious artifacts appearing in the polariton spectrum. Finally, the free photonic Hamiltonian of the cavity reads In Eqs. (3), (4), (5) in the main text, we only present the contribution from the TEM mode (m = 0), dropping the corresponding index. In Supplementary Note 1, we discuss the effect of the higher order (m ≠ 0) TM cavity modes for larger cavities.
Hopfield-Bogoliubov diagonalization. The polariton Hamiltonian H pol = H mat + H ph + H int , where H mat is given by Eq. (17), H ph by Eq. (24), and H int by Eqs. (20) and (21), can be recast into matrix form as H pol ¼ 1 Here, ψ y q ¼ ða y q ; b y q Þ is the spinor creation operator in the matter sublattice space and C y q ¼ ðc y q1 ; c y q2 ; ; c y qp ; ; c y qN Þ is the vector of TM photon creation operators, where p indexes the set of ordered triplets of integers (n 1 , n 2 , m), and N is the total number of photonic operators considered. The Hermitian [2(N + 2)] × [2(N + 2)] matrix H pol q can be written in block form as where is the (N + 2) × (N + 2) diagonal matrix of resonant frequencies of the free oscillators. The (N + 2) × (N + 2) block matrices H ± q can be sub-divided into block matrices where the upper-diagonal block is the 2 × 2 matrix in the matter subspace, and the lower-diagonal block H ph q is the N × N matrix in the photonic subspace with components Finally, the off-diagonal block H int q in Eq. (27) is the 2 × N interaction matrix, where the pth column reads The polariton Hamiltonian H pol is diagonalized by a generalized Hopfield-Bogoliubov transformation 48 Ψ q = T q X q , where X y q ¼ ðχ y q ; χ T Àq Þ and χ y q ¼ ðγ y q1 ; γ y q2 ; ; γ y qν ; ; γ y qNþ2 Þ. To ensure the invariance of the bosonic commutation relations for the transformed operators, T q must be a [2(N + 2)] × [2 (N + 2)] paraunitary matrix 70 that satisfies T q η z T y q ¼ T y q η z T q ¼ η z , where η z ¼ σ z 1 2 and σ z is the Pauli matrix. The transformed bosonic operators γ y qν ¼ Ψ y q η z jΨ qν i and γ qν = 〈Ψ qν |η z Ψ q that diagonalize the polariton Hamiltonian as create and annihilate polaritons in the vth band, respectively. The polariton dispersion ω pol qν (black solid lines in Fig. 2c-e) and the corresponding linearly independent eigenvectors |Ψ qν 〉 (first two columns of T q ) are determined from the positive eigenvalue solutions to the non-Hermitian eigenvalue equation η z H pol q jΨ qν i ¼ hω pol qν jΨ qν i.

Schrieffer-Wolff transformation.
To obtain an effective two-band Hamiltonian in the matter sublattice space, we neglect non-resonant terms in the matter Hamiltonian (since Ω=ω 0 ( 1 for practical realizations of the metasurface), but not in the light-matter interaction Hamiltonian since the photons are not resonant with the collective-dipoles near the corners of the Brillouin zone (see Fig. 2b). Next, we perform a unitary transformation and impose the Schrieffer-Wolff condition 50 where we have used the approximation jω ph qmn ± ω 0 j ) Ωjf q j that is valid near the K and K′ points. Retaining leading-order terms in ξ qmn , the transformed polariton where the matter and photonic subspaces are decoupled to quadratic order in ξ qmn . Calculating the commutator in Eq. (35) and extracting the Hamiltonian within the matter sublattice space, we obtain the two-band Hamiltonian H mat ¼ H mat À 2 h P qmn ξ 2 qmn ω ph qmn ðω ph qmn Þ 2 Àω 2 0 a y q a q þb y q b q þ ϕ 2 n b y q a q þ ϕ Ã2 n a y q b q : ð36Þ In Eq. (8) in the main text, we only present the contribution from the TEM mode (m = 0), dropping the corresponding index. We can recast the Hamiltonian (Eq. (36)) into matrix form as H mat ¼ P q ψ y q H mat q ψ q , with Bloch Hamiltonian Here W q ¼ω 0 À Ω P mn Δ qmn and F q ¼Ωf q À Ω P mn Δ qmn ϕ 2 n with Δ qmn ¼ Diagonalizing H mat leads to the two-band dispersion ω mat qτ ¼ W q þ τjF q j, which is indicated by the orange-dashed lines in Fig. 2c-e. The corresponding spinor eigenstates jψ qτ i ¼ ð1; τe iφ q Þ T = ffiffi ffi 2 p , where φ q ¼ argðF q Þ, can be represented by the pseudospin vector S qτ ¼ τðcos φ q ; sin φ q ; 0Þ from which we obtain the pseudospin vector field plots in Fig. 5a-c.
Expansion of the effective two-band Hamiltonian. Near the K point, the function Δ qmn , given by Eq. (38), expands as Kmn ðK À G n Þ x k x þ ðK À G n Þ y k y h i Kmn þ a 4 Δ ð2Þ Kmn þ a 4 Δ Kmn K À G n ð Þ 2 Kmn K À G n ð Þ x K À G n ð Þ y k x k y ð39Þ to quadratic order in k, where the real parameters Δ ðυÞ Kmn (υ = 0, 1, 2) depend only on the photon frequencies ω ph Kmn at the K point. Collecting the contributions from the degenerate photons (see Supplementary Note 2 for details), we obtain the effective Hamiltonian (Eq. (9)), where parameters are given by and C n ¼ 1 þ 3n 1 n 1 À 1 ð Þþ3n 2 n 2 À n 1 ð Þ: ð46Þ For brevity, we retain only the dominant (m = 0) TEM contribution for the plots in Fig. 4