I 1. Diagonalization of the Matter Hamiltonian

# Unconventional Dirac Polaritons in Cavity-Embedded Honeycomb Metasurfaces

The symmetries that dictate the existence of relativistic Dirac quasiparticles in condensed-matter systems (2) have been exploited in the realization of a plethora of artificial Dirac materials (5); (4); (6); (8); (7); (3); (9); (10); (11). In these artificial systems, the ability to design and manipulate the lattice structure has enabled the exploration of Dirac physics in new regimes (12); (14); (15); (16); (13). However, little attention has been paid to the effect of the surrounding environment on the nature of the Dirac quasiparticles. Here we theoretically investigate honeycomb arrays of meta-atoms embedded inside a planar photonic cavity. Massless Dirac polaritons emerge near the conventional Dirac points located at the corners of the Brillouin zone, in analogy with graphene (2). However, this analogy breaks down as the interaction with the photonic environment generates additional satellite Dirac points with Berry flux. Reducing the cavity height induces the merging of the satellite Dirac points with the conventional ones, forming a quadratic band-degeneracy with combined Berry flux. As a result, the massless Dirac polaritons with a linear spectrum morph into massive ones with a parabolic spectrum. Remarkably, this merging is not followed by Dirac point annihilation (17), but instead, massless Dirac polaritons re-emerge with an unprecedented inversion of chirality which has no analog in real or artificial graphene systems. This novel tunability could open up a new realm of unexplored Dirac-related physics, such as unconventional tunnelling and pseudo-magnetic related phenomena, in readily realizable experimental set-ups.

The discovery of monolayer graphene (18) has inspired an extensive quest for artificial graphene systems that mimic its underlying honeycomb symmetry. As a result, relativistic Dirac quasiparticles have emerged in an expanding catalog of distinct physical systems, ranging from ultracold atoms in optical lattices (12) to evanescently-coupled photonic waveguide arrays (14). The Dirac quasiparticles are endowed with an additional pseudospin degree of freedom due to the presence of two inequivalent sublattices forming the honeycomb structure (see Fig. 1). At the and points in the first Brillouin zone (see Fig. 2a), the two sublattices decouple giving rise to topological vortices in the pseudospin vector field (19). For systems with inversion and time-reversal symmetry, these vortices coincide with linear band-degeneracies which we call conventional Dirac points (CDPs). Since the Berry curvature is singular at the CDPs, the latter are sources of quantized Berry flux , where is the pseudospin winding number (i.e. the topological charge of the CDPs) (19). Consequently, the massless Dirac quasiparticles acquire a topological Berry phase when adiabatically transported around a loop in momentum space that encloses a CDP, giving rise to remarkable properties such as the anomalous quantum Hall effect (20); (21) and the iconic Klein tunnelling effect (22).

To qualitatively alter the fundamental properties of Dirac quasiparticles, one conventionally modifies the underlying lattice structure and symmetries that dictate their existence. For example, breaking the inversion or time-reversal symmetry lifts the degeneracy at the CDPs (23); (24), giving rise to massive Dirac quasiparticles with an energy-dependent Berry phase (19). Alternatively, introducing lattice anisotropy via uniaxial strain can induce the merging of CDPs with opposite Berry flux (25). This merging is followed by the annihilation of the CDPs and the corresponding Berry curvature in a topological phase transition that gaps out the quasiparticle spectrum (17). Furthermore, inhomogeneous strain can create large pseudomagnetic fields which give rise to quantized Landau levels (26). These regimes have been the focus of many artificial graphene studies which exploit the meticulous control over the underlying lattice structure (12); (14); (15); (16); (13).

Here we investigate a conceptually different scenario by exploring the fate of the Dirac quasiparticles under the influence of a tunable photonic environment whilst preserving the underlying symmetry. In particular, we theoretically investigate the polariton modes supported by a cavity-embedded metasurface that is composed of a honeycomb array of identical meta-atoms (see Fig. 1). These hybrid states, resulting from the coupling between light and matter degrees of freedom (27), are endowed with signatures of the honeycomb structure giving rise to CDPs in the polariton spectrum at the and points. However, the polariton modes are inextricably linked to the surrounding photonic environment which has a non-trivial effect on the Berry curvature. Specifically, additional linear band-degeneracies with Berry flux are generated—which we call satellite Dirac points (SDPs)—whose location within the Brillouin zone can be manipulated by adjusting the cavity height. Increasing the strength of the light-matter interaction by reducing the cavity height drives the SDPs towards the CDPs which remain pinned at the and points. At a critical cavity height, the SDPs merge with the CDPs forming a quadratic band-degeneracy with combined Berry flux. Crucially, the SDPs do not annihilate the CDPs since the point-group symmetry is preserved. Instead, by decreasing the cavity height further, massless Dirac polaritons with a linear spectrum re-emerge near the CDPs accompanied by an unprecedented inversion of chirality.

The essential physics is captured by a minimal model of the metasurface which is analytically tractable. We model each meta-atom with a single dynamical degree of freedom describing the electric-dipole moment associated with its (non-degenerate) fundamental eigenmode, with resonant frequency . Furthermore, we consider sub-wavelength nearest-neighbor separation , and orient the meta-atoms such that their corresponding dipole moments point normal to the plane of the lattice. The metasurface is embedded at the centre of a planar photonic cavity of height (see Fig. 1), where the cavity walls are assumed to be lossless and perfectly-conducting metallic plates. This general model can be readily realized in a variety of experimental set-ups including, e.g., arrays of plasmonic nanorods or microwave helical resonators (see Supplementary Section 5 for numerical simulations of these systems).

Within the Coulomb gauge the polariton Hamiltonian reads , where the matter Hamiltonian (see Methods) within the rotating-wave approximation (RWA) reads

 HRWAmat=ℏω0∑q{[1−ΩS(L)](a†qaq+b†qbq)+Ω[1−I(L)](fqb†qaq+f∗†qa†qbq)}. (1)

Equation (1) is a bosonic analog of the electronic tight-binding Hamiltonian in graphene (2). In this analogy, the kinetic term related to hopping of electrons between carbon atoms is replaced by Coulomb dipole-dipole interactions between the meta-atoms, where parametrizes the strength of the nearest-neighbor interaction. We neglect Coulomb interactions beyond nearest-neighbors as they do not qualitatively affect the physics near the and points and are negligible for small cavity heights (see Supplementary Section 3 where we retain all interactions). Coulomb interactions with the cavity-induced image dipoles are encoded in the parameters

 Extra open brace or missing close brace (2)

which renormalize the resonant frequency and interaction strength , respectively. The bosonic operators and create quanta of the quasistatic collective-dipole modes that extend across the and sublattices, respectively, with wavevector in the first Brillouin zone (see Fig. 2a). The function encodes the honeycomb geometry of the lattice with nearest-neighbor vectors (see Fig. 1). Diagonalizing (see Supplementary Section 1 for details) leads to the quasistatic collective-dipole dispersion , where indexes the upper () and lower () bands. As observed in Fig. 2b, this dispersion resembles the electronic band structure of graphene (2) with two inequivalent CDPs at the () and () points located at , where indexes the two valleys.

The free photonic environment inside the cavity is described by the Hamiltonian

 Hph=∑qmnℏωqmnc†qmncqmn (3)

while the light-matter interaction Hamiltonian (see Methods and Supplementary Section 2 for details) reads

 Missing \left or extra \right (4)

Here, the bosonic operator creates a transverse magnetic (TM) cavity photon with dispersion , where is a non-negative integer indexing the different TM cavity modes, and indexes the TM photons associated with the reciprocal lattice vector . Only TM photons with even couple to the dipoles due to the parity selection rule (see Methods), and we neglect transverse electric photons whose polarization is purely in-plane (orthogonal to dipole moments). The strength of the light-matter interaction in equation (4) is parametrized by

 ξqmn=F(ωqmn)ωq0nωqmn(8π3√3NmΩaLω0ωqmn)12, (5)

where . The dipole approximation employed in equation (4) is valid only for wavelengths larger than the size of the meta-atom. Therefore, to take into account the finite size of the meta-atoms we have introduced a phenomenological function that provides a smooth cut-off for the interaction with short-wavelength photonic modes (see Methods for expression). The complex phase factors in equation (4) are associated with Umklapp processes which must be retained as they are critical for maintaining the point-group symmetry of the polariton Hamiltonian. 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. This allows one to continuously increase the light-matter interaction strength with the TEM photons by decreasing , while the higher-order TM photons () become increasingly detuned with the dipole resonances, i.e., for small cavity heights. We diagonalize using a generalized Hopfield-Bogoliubov transformation (27) (see Methods), and in Fig. 2c-e we present the resulting polariton dispersion for different cavity heights.

Naively, one may expect all of the band crossings between the photon and quasistatic collective-dipole dispersions to be avoided. This is the case for the crossings with the upper quasistatic band where the constructive interference between the sublattices of this bright () configuration results in a large characteristic anticrossing for all wavevector directions (see Fig. 2f and Fig. 2h). In contrast, the lower quasistatic band corresponds to a dark () configuration which results in a small anticrossing for a general wavevector direction (see Fig. 2i). Crucially, the light-matter interaction vanishes along the high-symmetry directions (see Supplementary Section 2) due to the perfect destructive interference between the two sublattices (see Fig. 2g). As a result of this protected band crossing for wavevector directions perpendicular to the nearest-neighbor directions , six inequivalent SDPs emerge in the polariton spectrum within the first Brillouin zone. Near the SDPs, the dispersion takes the form of an anisotropic and critically-tilted Dirac cone (see inset of Fig. 2c).

Fig. 2c-e highlights the evolution of the polariton spectrum into qualitatively distinct phases as the cavity height is reduced. To elucidate the topological properties of the band-degeneracies and gain insight into the nature of the polaritons near the and points, we first perform a unitary Schrieffer-Wolff transformation (28) (see Methods) and extract the effective two-band Hamiltonian in the matter subspace

 ¯Hmat=HRWAmat−2ℏω20∑qmnξ2†qmnωqmnω2†qmn−ω20[1−ΩS(L)]2(a†qaq+b†qbq+ϕ2†nb†qaq+ϕ∗2†na†qbq). (6)

Here the photonic degrees of freedom have been integrated out to obtain an effective dipole-dipole interaction to quadratic order in . These quadratic renormalizations quantitatively capture the evolution of the polariton spectrum as shown by the orange-dashed lines in Fig. 2c-e (see Methods). Next, we expand the two-band Hamiltonian (6) near the and points and, up to quadratic order in (with ), obtain the effective Hamiltonian (see Methods and Supplementary Section 4 for details) with spinor creation operator and Bloch Hamiltonian

 ¯Heffkζ=ℏω0[1−Δ0(L)]12−ℏv|k|[1−Δ1(L)]σ⋅^n1ζ+ℏ2|k|22m[1−Δ2(L)]σ⋅^n2ζ−ℏ2|k|22m∗Δ3(L)12. (7)

In equation (7) is the vector of Pauli matrices operating in the sublattice space and is the identity matrix. The projection unit vectors for the linear and quadratic terms in are defined as and , respectively, where , and the parameters , , and are fixed by the lattice structure. As shown in Fig. 3, the tunable parameters () are all real, positive, and increase as the cavity height reduces (see Methods for expressions). In Fig. 4a-c we plot the pseudospin vector field near the point for different cavity heights (see Methods) and schematically depict the location of the band-degeneracies along with their associated Berry flux. In Fig. 4d-f we plot the corresponding effective polariton spectrum to leading order in . Despite the preserved honeycomb symmetry of the metasurface, the generalized Dirac Hamiltonian (7), and thus the very nature of its chiral quasiparticles, are highly tunable with the strength of the interaction with the photonic environment.

For sub-critical cavity heights (), three SDPs are located along the directions, each with Berry flux surrounding the CDP with Berry flux (see Fig. 4a). The polariton spectrum disperses linearly near the CDPs (see Fig. 2c) forming Dirac cones with a group velocity that is tunable via the cavity height (see Fig. 3). This is in stark contrast to graphene where the Fermi velocity is fixed by the lattice structure. To leading order in , the effective Hamiltonian (7) is equivalent to a 2D massless Dirac Hamiltonian with spinor eigenstates , where indexes the upper () and lower () polariton bands. These represent massless Dirac polaritons with chirality , resulting in a pseudospin that winds once around the CDPs and a topological Berry phase of (see Fig. 4d).

The coupling to the environment induces a redshift of the CDP energy (see Fig. 3) allowing one to investigate the polaritonic analogue of Klein tunnelling (22) at interfaces separating regions with different cavity heights. One could exploit this for polaritonic Veselago lensing in the same spirit as negative refraction in graphene (29). However, here the effective negative index is induced by a simple variation in cavity height as opposed to the precise tuning of the Fermi energy that is required in graphene (30).

At the critical cavity height () where (see Fig. 3), the group velocity of massless Dirac polaritons vanishes as the SDPs merge with the CDP forming a quadratic band-degeneracy (see Fig. 2d) with combined Berry flux (see Fig. 4b). The leading order term in the effective Hamiltonian (7) is now quadratic in with corresponding spinor eigenstates . Therefore, during this merging transition the massless Dirac polaritons morph into massive Dirac polaritons with chirality , resulting in a pseudospin that winds twice as fast compared to the sub-critical case and a topological Berry phase of (see Fig. 4e).

Since the point-group symmetry is preserved, the SDPs do not annihilate the CDPs but re-emerge and continue to migrate along the directions as the cavity height is reduced past criticality (see inset of Fig. 4c). After a small decrease in cavity height, SDPs with opposite Berry flux merge together and annihilate their corresponding Berry curvature in a topological transition that leaves only the Berry flux from the CDPs (see Fig. 2e and Fig. 4c). For super-critical coupling (), we recover the linear dispersion near the CDPs to leading order in (see Fig. 2e) and the effective Hamiltonian (7) is equivalent to a 2D massless Dirac Hamiltonian with corresponding spinor eigenstates . Remarkably, massless Dirac polaritons thus re-emerge past criticality but with an inversion of chirality , i.e., a flip in pseudospin direction (see Fig. 4f). One could observe this directly by scanning the near-fields (31) and measuring the change in phase linking the two sublattices, or indirectly by measuring the inversion of trigonal warping ( rotation of the isofrequency domains) that accompanies the chirality inversion (cf. Fig. 4a and Fig. 4c).

Despite having identical lattice symmetries, the honeycomb metasurface studied here does not provide a trivial analogy with graphene. The hybridization with the dynamic photonic environment has profound consequences on the Berry curvature and the fundamental nature of the emergent Dirac polaritons. Exploiting the photonic environment—rather than the lattice structure—as a means to modify the fundamental properties of Dirac polaritons, could open a new realm of unexplored Dirac-related phenomena which have no analogue in real or artificial graphene systems. For example, novel transport phenomena could occur through interfaces separating regions with different cavity heights. In particular, the predicted inversion of chirality would have drastic consequences for the Klein tunnelling effect (22) since it removes the conventional channel responsible for the perfect transmission (see Fig. 4d and Fig. 4f), and could give rise to strongly valley-dependent transport properties. Furthermore, combining the effects of the cavity with inhomogeneous strain deformations could lead to unique pseudo-magnetic related effects, such as a Landau level spectrum and novel quantum Hall states that are tunable via the environment.

## Acknowledgments

C.-R.M. and T.J.S. acknowledge financial support from the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom. In addition C.-R.M acknowledges the EPSRC Centre for Doctoral Training in Metamaterials (Grant No. EP/L015331/1) and QinetiQ. G.W. acknowledges financial support from Agence Nationale de la Recherche (Project ANR-14-CE26-0005 Q-MetaMat) and the Centre National de la Recherche Scientifique through the Projet International de Coopération Scientifique program (Contract No. 6384 APAG). W.L.B. acknowledges the financial support from EPSRC (Grant No. EP/K041150/1). E.M. acknowledges financial support from the Leverhulme Trust (Research Project Grant RPG-2015-101), and the Royal Society (International Exchange Grant No. IE140367, Newton Mobility Grants 2016/R1 UK-Brazil, and Theo Murphy Award TM160190).

## Author contributions

C.-R.M. performed the theoretical calculations and wrote the manuscript with E.M.; T.J.S. contributed to the theoretical calculations; E.M. and G.W. conceived the project; and E.M. and W.L.B. supervised the project. All authors commented on the manuscript.

Supplementary information accompanies this work. Correspondence and requests for materials should be addressed to C.-R.M or E.M.

## Competing financial interests

The authors declare no competing financial interests.

## methods

Derivation of polariton Hamiltonian. The cavity-embedded metasurface is composed of a honeycomb array of identical meta-atoms located at and on the the and sublattices, respectively. Here, is an in-plane lattice translation vector with primitive vectors and (see Fig. 1), and integers and . Each meta-atom supports a non-degenerate fundamental eigenmode with resonant frequency , and we neglect all higher-order eigenmodes. We employ the Coulomb gauge where the non-dynamical scalar potential is eliminated in favor of the instantaneous Coulomb interaction between the meta-atoms. The nearest-neighbor separation is considered to be large enough such that the Coulomb interaction is dominated by the dipole-dipole term. Therefore, each meta-atom is modelled by a single dynamical degree of freedom (with dimensions of length) where the dipole moment associated with its fundamental eigenmode is with effective charge . The Coulomb potential energy between two dipole moments and located at and , respectively, is given by

 VCoul=p⋅p′−3(p⋅^n)(p′⋅^n)4πε0|R−R′|3, (8)

where and is the vacuum permittivity. The presence of the perfectly-conducting metallic plates, placed at and , modifies the boundary conditions on the scalar potential. Using the method of images to ensure the vanishing of the scalar potential at the cavity walls (33), we introduce an infinite series of image dipoles located outside the cavity at positions , where indexes the two sublattices and is a non-zero integer. Noting that the Coulomb potential energy between a real and image dipole is of that given by equation (8) (34), the matter Hamiltonian within the nearest-neighbor approximation reads

 Misplaced & (9)

Here, the dynamical coordinate and conjugate momentum corresponding to the meta-atom located at satisfy the canonical commutation relations and , and is an effective mass. Next we introduce the bosonic operators

 aRA=√Mω02ℏhA(RA)+i√12ℏMω0ΠA(RA) (10)

and

 bRB=√Mω02ℏhB(RB)+i√12ℏMω0ΠB(RB) (11)

that annihilate quanta of the fundamental eigenmode on the meta-atom located at and on the and sublattices, respectively, and satisfy the commutation relations

 [aRA,a†R′A]=δRAR′A,[bRB,b†R′B]=δRBR′B,[aRA,b†R′B]=0. (12)

In terms of these operators the matter Hamiltonian (9) reads

 Missing or unrecognized delimiter for \Biggl (13)

where parametrizes the strength of the nearest-neighbor Coulomb interaction, and we assume . We apply Born-von Kármán boundary conditions over a lattice with unit cells and introduce the Fourier transform of the bosonic operators

 aRA=1√N∑qaqeiq⋅RA,bRB=1√N∑qbqeiq⋅RB (14)

which transforms the matter Hamiltonian (13) into a local and block-diagonal form

 Missing or unrecognized delimiter for \Biggl (15)

Finally, neglecting non-resonant terms in equation (15) we obtain the matter Hamiltonian (1) within the RWA.

In the Coulomb gauge the light-matter interaction is described by the minimal-coupling Hamiltonian (35) which, within the dipole approximation, reads

 Hint=QM∑s=A,B∑RsΠs(Rs)Az(Rs)HΠ⋅A+Q22M∑s=A,B∑RsA2z(Rs)HA2, (16)

where we have used . Due to the discrete in-plane translational symmetry, the interactions between the photons and quasistatic collective-dipole modes only conserve in-plane momentum modulo a reciprocal lattice vector. Therefore, it is convenient to expand the -component of the vector potential (see Supplementary Section 2 for details) as

 Az(r,z)=∑qmn√ℏε0NmNALωqmn|q−Gn||q−Gn+mπL^z|cos(mπLz)[cqmnei(q−Gn)⋅r+c†qmne−i(q−Gn)⋅r], (17)

where is the area of a unit cell. Here, is a reciprocal lattice vector with primitive vectors and , and is a single integer indexing the set of ordered pairs of integers and . From equation (17) one observes that the dipoles do not couple to TM modes with odd due to the odd parity about the centre of the cavity where the meta-atoms are located. Substituting the vector potential (17) into equation (16) we obtain the light-matter interaction Hamiltonian (4). We choose the phenomenological cut-off function in the light-matter interaction parameter (5) to be of the Fermi-Dirac distribution form

 F(ωqmn)=11+e2(ωqmn−3ω0)/ω0, (18)

which is smooth enough to avoid spurious artifacts appearing in the polariton spectrum. The main physical reason motivating the introduction of a cut-off is the break down of the dipole-approximation employed in equation (4) due to the finite size of the meta-atoms.

Hopfield-Bogoliubov diagonlization. The polariton Hamiltonian , where is given by equation (15), by equation (3), and by equation (4), can be recast into matrix form as where . Here is the spinor creation operator in the matter sublattice space and is the vector of creation operators for TM cavity photons, where is a single integer indexing the set of ordered triplets of integers , , and . The Hermitian matrix (where is the number of photon operators) can be written in block form

 Hpolq=(H+qH−q−Wq(H−q−Wq)†(H+−q)∗), (19)

where

 Wq=ℏDiag(ω0,ω0,ωq1,ωq2,…,ωqp,…) (20)

is the diagonal matrix of resonant frequencies of the free oscillators. The block matrices can be sub-divided into block matrices

 H±q=⎛⎜⎝H% matqHintq±(Hintq)†Hphq⎞⎟⎠, (21)

where the upper-diagonal block is the matrix in the matter subspace

 Unknown environment '% (22)

and the lower-diagonal block is the matrix in the photonic subspace with components

 Extra open brace or missing close brace (23)

Finally, the off-diagonal block in equation (21) is the interaction matrix, where the column reads

 Unknown environment '% (24)

The polariton Hamiltonian is diagonalized by a generalized Hopfield-Bogoliubov transformation (27) , where and . To ensure the invariance of the bosonic commutation relations for the transformed operators, must be a para-unitary matrix (36) that satisfies , where and is the Pauli matrix. The transformed bosonic operators and that diagonalize the polariton Hamiltonian as

 Hpol=∑qνℏωqνγ†qνγqν, (25)

create and annihilate polaritons in the band, respectively. The polariton dispersion (black solid lines in Fig. 2c-e) and the corresponding linearly-independent eigenvectors (first two columns of ) are determined from the positive eigenvalue solutions to the non-Hermitian eigenvalue equation .

Schrieffer-Wolff transformation. To obtain an effective two-band Hamiltonian in the matter sublattice space, we first perform the RWA in the matter Hamiltonian (15) but not in the light-matter interaction Hamiltonian (4) since the photons are not resonant with the dipoles near the corners of the Brillouin zone (see Fig. 2b). Next, we perform a unitary transformation

 ¯Hpol=eSHpole−S=Hpol+[S,H%pol]+12[S[S,Hpol]]+… (26)

and impose the Schrieffer-Wolff condition (28)

 [S,Hmat+Hph]=−HΠ⋅A (27)

which eliminates the light-matter interaction to first order in . From equation (27) the particular form of the anti-Hermitian operator reads

 S=−∑qmniω0ξqmnωqmn−¯ω0(ϕ∗†na†q+ϕnb†q)cqmn+∑qmniω0ξqmnωqmn+¯ω0(ϕ∗†na†q+ϕnb†q)c†−qmn−H% .c. (28)

where we have defined and used the approximation that is valid near the and points. Retaining leading-order terms in , the transformed polariton Hamiltonian (26) reads

 ¯Hpol≃Hmat+12[S,HΠ⋅A]¯Hmat+Hph+HA2¯Hph, (29)

where the matter and photonic subspaces are decoupled to quadratic order in . Calculating the commutator in equation (29) and extracting the Hamiltonian within the matter subspace we obtain the effective two-band Hamiltonian (6). We can recast the Hamiltonian (6) into matrix form as with Bloch Hamiltonian

 ¯Hmatq=ℏω0(WqΩF∗qΩFqWq). (30)

Here and with

 δqmn(L)=16π3√3Nm(aL)⎛⎝ω20ω2†qmn−¯ω20⎞⎠(ωq0nωqmn)2F2(ωqmn). (31)

Diagonalizing leads to the effective two-band polariton dispersion (orange dashed lines in Fig. 2c-e). The corresponding spinor eigenstates , with , can be represented by a pseudospin vector on the Bloch sphere

 Unknown environment '% (32)

from which we obtain the pseudospin vector field plots in Fig 4a-c.

Expansion of the effective two-band Hamiltonian. Near the and points the function , given by equation (31), expands as

 δkmn(L)≃δ(0)Kmn−a2δ(1)Kmn[ζ(K−Gn)xkx+(K−Gn)yky]+12[−a2δ(1)Kmn+a4δ(2)Kmn(K−Gn)2x]k2x+12[−a2δ(1)Kmn+a4δ(2)Kmn(K−Gn)2y]k2y+ζa4δ(2)Kmn(K−Gn)x(K−Gn)ykxky (33)

to quadratic order in , where the real parameters () depend only on the photon frequencies at the and points. Collecting the contributions from the degenerate photons (see Supplementary Section 4 for details) we obtain the tunable Dirac Hamiltonian (7), where the tunable parameters () are given by

 Δ0(L)=Ω(S(L)+∑mnδ(0)Kmn(L)), (34)
 Δ1(L)=I(L)+4π27∑mnAnδ(1)Kmn(L), (35)
 Δ2(L)=I(L)+8π281∑mnBnδ(2)Kmn(L), (36)

and

 Δ3(L)=Ω∑mn(4π227Cnδ(2)Kmn(L)−12δ(1)Kmn(L)), (37)

with

 An=√32(2−3n1)cos[4π3(n1+n2)]+12(6n2−3n1)sin[4π3(n1+n2)], (38)
 Missing or unrecognized delimiter for \right (39)

and

 Cn=1+3n1(n1−1)+3n2(n2−n1). (40)

For brevity, we retain only the dominant contribution from the TEM mode (), where the coefficients in equation (33) are given by

 δ(0)K0n(L)=2(aL)(a|K|)⎛⎝ω20ω2†K0n−¯ω20⎞⎠F2(ωK0n), (41)
 Missing or unrecognized delimiter for \left (42)

and

 δ(2)K0n(L)=16(aL)(a|K|)−3(ωK00ω0)4⎛⎝ω20ω2†K0n−¯ω20⎞⎠3⎧⎪ ⎪⎨⎪ ⎪⎩F2(ωK0n)−∂[F2(ωK0n)]∂ωK0n⎡⎢ ⎢⎣(5ω2†K0n−¯ω20)(ω2†K0n−¯ω20)8ω3†K0n⎤⎥ ⎥⎦+∂2[F2(ωK0n)](∂ωK0n)2⎡⎢ ⎢ ⎢⎣(ω2†K0n−¯ω20)28ω2†K0n⎤⎥ ⎥ ⎥⎦⎫⎪ ⎪ ⎪⎬⎪ ⎪ ⎪⎭. (43)

## References

Supplementary Information: Unconventional Dirac Polaritons in Cavity-Embedded Honeycomb Metasurfaces

This supplementary material is organized as follows. In Sec. 1 we diagonalize the matter Hamiltonian (15) (see Methods) via a Bogoliubov transformation and obtain an effective massless Dirac Hamiltonian describing quasistatic collective-dipole modes near the and points. In Sec. 2 we discuss the electromagnetic modes supported by the planar metallic cavity. In Sec. 3 we include all Coulomb interactions in the matter Hamiltonian and discuss their effects on the polariton dispersion. In Sec. 4 we expand the two-band Hamiltonian (6) (see main text) and obtain the effective Dirac Hamiltonian (7) (see main text) describing polaritons near the and points. Finally, to illustrate the applicability of this general model, in Sec. 5 we present numerical simulations of the polariton dispersion for honeycomb arrays of plasmonic nanorods and microwave helical resonators.

## Appendix I 1. Diagonalization of the Matter Hamiltonian

The matter Hamiltonian (15) within the nearest-neighbor approximation (see Methods) can be recast into matrix form as , where and . The Hermitian matrix can be written in block form as

 Hmatq=(HmatqHmatq−ℏω012Hmatq−ℏω012Hmatq), (S1)

where

 Unknown environment '% (S2)

We diagonalize via a Bogoliubov transformation , where and , with being a paraunitary matrix (see Methods). The Bogoliubov operators and