ARTICLE DOI: 10.1038/s41467-018-03982-7 OPEN Manipulating type-I and type-II Dirac polaritons in cavity-embedded honeycomb metasurfaces 1 1 2 1 1 Charlie-Ray Mann , Thomas J. Sturges , Guillaume Weick , William L. Barnes & Eros Mariani Pseudorelativistic Dirac quasiparticles have emerged in a plethora of artiﬁcial graphene systems that mimic the underlying honeycomb symmetry of graphene. However, it is notoriously difﬁcult 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 artiﬁcial graphene systems. Exploiting the photonic environment will thus give rise to unexplored Dirac physics at the subwavelength scale. 1 2 EPSRC Centre for Doctoral Training in Metamaterials (XM2), Department of Physics and Astronomy, University of Exeter, Exeter EX4 4QL, UK. Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504, F-67000 Strasbourg, France. Correspondence and requests for materials should be addressed to C.-R.M. (email: cm433@exeter.ac.uk) or to E.M. (email: E.Mariani@exeter.ac.uk) NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications 1 | | | 1234567890():,; ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 he groundbreaking discovery of monolayer graphene has photonic cavity and simply changing the cavity height, one can inspired an extensive quest to emulate massless Dirac induce multiple phase transitions including the multimerging of Tquasiparticles in a myriad of distinct artiﬁcial graphene type-I and type-II Dirac points and the annihilation of type-II Dirac 2–11 3 systems , ranging from ultracold atoms in optical lattices to points. This striking tunability results in qualitatively different evanescently coupled photonic waveguide arrays . Owing to their polariton phases, despite the preserved lattice structure. In parti- honeycomb symmetry, linear band-degeneracies manifest in the cular, we unveil a morphing between a linear and a parabolic quasiparticle spectrum which we call conventional Dirac points spectrum accompanied by a change in the topological Berry phase, (CDPs). These belong to the ubiquitous type-I class of two- and an environment-induced inversion of chirality, all of which dimensional (2D) Dirac points that are characterized by Dirac have no analog in graphene or artiﬁcial graphene systems studied cones with closed isofrequency contours. As a result, the corre- thus far. Therefore, this unique paradigm will give rise to unex- sponding quasiparticles are described by the rather exotic 2D plored Dirac-related phenomena at the subwavelength scale, such as massless Dirac Hamiltonian , and thus offer fundamental insight anomalous Klein tunneling, negative refraction, and pseudomag- into pseudorelativistic phenomena such as the iconic Klein netic Landau levels, which can all be tuned via the photonic paradox . The latter is responsible for the suppression of back- environment alone. scattering and for the antilocalization of massless Dirac quasi- particles, which are highly desirable properties for efﬁcient Results quasiparticle propagation in novel devices. Hamiltonian formulation. While metamaterials have tradition- Since the existence of type-I CDPs is intrinsically linked to the ally been described in terms of macroscopic effective 30,33,43 honeycomb structure, the fundamental properties of the massless properties , the importance of crystallinity is becoming Dirac quasiparticles are notoriously robust and difﬁcult to increasingly apparent . Therefore, to capture the essential manipulate. However, by exploiting meticulous control over the physics related to complex non-local effects that arise from lattice structure, artiﬁcial graphene systems have enabled the strong multiple-scattering , here we study the properties of the exploration of Dirac quasiparticles in new regimes that are dif- cavity-embedded honeycomb metasurface by means of a 14–19 ﬁcult, if not impossible to achieve in graphene itself . Among microscopic Hamiltonian formalism. This allows us to clearly others, the archetypal example which has attracted considerable identify the distinct physical origins of the type-I and type-II interest is the paradigm of strain-engineering, where it has been Dirac points. shown that lattice anisotropy can induce the merging and anni- The full polariton Hamiltonian of this system reads H = pol 3,14–16,20–23 hilation of type-I CDPs , and that aperiodicity can H + H + H , where the interaction Hamiltonian H mat ph int int 17,24 generate large pseudomagnetic ﬁelds . couples the matter and photonic subspaces whose free dynamics Moreover, the recent discovery of type-II Dirac/Weyl semi- are governed by H and H , respectively. We employ the mat ph 25–29 metals sparked a burgeoning exploration into the prospects Coulomb gauge, where the instantaneous Coulomb interaction of a rarer type-II class of three-dimensional Dirac/Weyl points. between the meta-atoms is incorporated within the matter As the latter are characterized by critically tilted Dirac/Weyl Hamiltonian H , and the effects of the dynamic photonic mat cones with open, hyperbolic isofrequency contours, the corre- environment—described by the transverse vector potential—are sponding Lorentz-violating Dirac/Weyl quasiparticles exhibit included through the principle of minimal-coupling . 25–29 markedly different properties from their type-I counterparts . A schematic of a cavity-embedded honeycomb metasurface is 30– Soon after their realization, electromagnetic analogs emerged depicted in Fig. 1. We model each subwavelength meta-atom by a , and this exploration has recently been extended to 2D systems single dynamical degree of freedom describing the electric-dipole where a distinct type-II class of 2D Dirac points were theoretically moment associated with its (non-degenerate) fundamental 35,36 predicted . However, since their existence is predicated on eigenmode with resonant frequency ω . These meta-atoms are strong anisotropy in judiciously engineered photonic structures, then oriented such that their dipole moments point normal to the one cannot manipulate their properties without modifying the plane of the lattice. Furthermore, we consider subwavelength lattice structure. nearest-neighbor separation a such that the light cone intersects This hunt for exotic quasiparticles has recently entered the the Brillouin zone edge above ω , ensuring the existence of 37–42 realm of polaritonics . The true potential of polaritons lies in evanescently bound, subwavelength polaritons. The strength of their hybrid nature, where their light and matter constituents can the Coulomb dipole–dipole interaction between neighboring be manipulated independently, thereby providing additional meta-atoms is parametrized by Ω. Finally, the metasurface is tunable degrees of freedom. Among other examples, recent works embedded at the center of a planar photonic cavity of height L, have shown the tantalizing prospect of engineering novel topo- where the cavity walls are assumed to be lossless and perfectly logical polaritons by introducing a winding coupling between conducting metallic plates. Such a structure is imminently 39,41 ordinary photons and excitons . realizable across the electromagnetic spectrum from arrays of In this work, we exploit the hybrid nature of polaritons in a plasmonic nanoparticles to microwave helical resonators (see different setting, namely metasurfaces, and we unveil unique Fig. 1). Dirac physics by shifting the focus from the lattice structure and its deformations to the effect of manipulating the sur- Emergence of type-I Dirac points. The matter Hamiltonian rounding photonic environment. In particular, we theoretically within the nearest-neighbor approximation reads study the polaritons supported by imminently realizable, crys- X X y y y talline metasurfaces consisting of a honeycomb array of reso- H ¼ hω ~ a a þ b b þ hΩ f b a þ H:c: ; mat 0 q q q q q q q ð1Þ nant, dipolar meta-atoms. Despite theelementarynatureof q q these metasurfaces, we unveil the simultaneous existence of both type-I and type-II massless Dirac polaritons which have where, for brevity, we have not presented the non-resonant terms distinct physical origins. Crucially, the existence of the latter is (see Methods for derivation). In Eq. (1), ω is the renormalized not a result of anisotropy but is intrinsically linked to the resonant frequency and Ω is the renormalized Coulombic inter- hybrid nature of the polaritons, emerging from a non-trivial action strength due to the cavity-induced image dipoles (see winding in the light–matter interaction. Furthermore, we show Methods for their dependence on the cavity height). The y y that by embedding the honeycomb metasurface inside a planar bosonic operators a and b create quanta of the quasistatic q q 2 NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 ARTICLE y To quadratic order in k = q−K (ka 1), the effective matter eff eff Hamiltonian near the K point is H ¼ ψ H ψ , with K K;k k k k y y y spinor creation operator ψ ¼ a ; b and Bloch Hamiltonian k k k eff H ¼ hω ~ 1 h~vσ k þ htðÞ σ k ð2Þ K;k 0 2 Here, 1 is the 2 × 2 identity matrix, σ = (σ , σ ) and σ = 2 x y °2 (σ , −σ ) are vectors of Pauli matrices, and represents the x y Hadamard (element-wise) square. Note that the image dipoles do B not qualitatively affect the physics, but simply lead to a renormalization of the group velocity ~v ¼ 3Ωa=2 and trigonal a a 2 1 2 ~ warping parameter t ¼ 3Ωa =8. Apart from a global energy shift, Meta-atom Eq. (2) is equivalent to a 2D massless Dirac Hamiltonian to leading order in k, with an isotropic Dirac cone spectrum mat ω ¼ ω ~ þ τ~vjj k that is characterized by closed isofrequency Possible experimental Microwave Plasmonic kτ 0 realizations helical resonators nanorods contours. Therefore, as expected from the honeycomb symmetry, the CDP belongs to the type-I class of 2D Dirac points, and the pﬃﬃﬃ iθ corresponding spinors jψ i¼ð1; τe Þ = 2, where θ = kτ Fig. 1 Schematic of a cavity-embedded honeycomb metasurface. The arctan(k /k ), represent massless Dirac collective-dipoles with a y x honeycomb array of meta-atoms is composed of two inequivalent (A pﬃﬃ topological Berry phase of π. The effective Hamiltonian near the 3 3 and B) hexagonal sublattices—deﬁned by lattice vectors a ¼ að ; Þ pﬃﬃ 1 2 2 eff eff K′ point is given by H ¼ðH Þ , where the corresponding 3 3 K′;k K;k and a ¼ að ; Þ—which are connected by nearest-neighbor vectors 2 2 pﬃﬃ pﬃﬃ massless Dirac collective-dipoles have a topological Berry phase 3 1 3 1 e = a(0, −1), e ¼ að ; Þ, and e ¼ að ; Þ, where a is the 2 2 2 3 2 2 of −π, as required by time-reversal symmetry. subwavelength nearest-neighbor separation. Each subwavelength meta- atom is modeled as an electric dipole, oriented normal to the plane of the Hybridization with the photonic environment. Given the sub- lattice. The honeycomb metasurface is then embedded inside a photonic wavelength nearest-neighbor separation, it is tempting to assert cavity of height L, which is composed of two perfectly conducting metallic that the near-ﬁeld Coulomb interactions in H capture the mat plates, enabling one to modify the photonic environment while preserving essential physics. In fact, we will show that this quasistatic the lattice structure. This general model can be readily realized across the description misses the profound inﬂuence of the surrounding electromagnetic spectrum, from arrays of plasmonic nanorods to photonic environment, which has a remarkably non-trivial effect microwave helical resonators on the Berry curvature and, therefore, on the corresponding nature of the polaritons. collective-dipole modes that extend across the A and B sub- Crucially, the metallic cavity supports a fundamental transverse lattices, respectively, with wavevector q in the ﬁrst Brillouin zone electromagnetic (TEM) mode whose polarization (parallel to (see Fig. 2a). Finally, the function f ¼ expðiq e Þ encodes q j¼1 j the dipole moments) and linear dispersion (see Fig. 2b) are the honeycomb geometry of the lattice with nearest-neighbor independent of the cavity height. For brevity, in what follows we vectors e (see Fig. 1). do not present the contributions from the other cavity modes We diagonalize the matter Hamiltonian (Eq. (1)) as since the essential physics emerges from the interaction with the P P mat H ¼ hω β β where the bosonic operators fundamental TEM mode (see Methods for the full expressions). mat τ¼ ± q qτ qτ qτ In fact, the higher order cavity modes become increasingly β ¼ ψ jψ i create quasistatic collective-dipole normal modes qτ q qτ negligible for smaller cavities as they are progressively detuned mat with dispersion ω ¼ ω ~ þ τΩjf j. Here, τ indexes the upper qτ 0 q from the dipole resonances. y y y (τ=+1) and lower (τ = −1) bands and ψ ¼ða ; b Þ is a spinor q q The effects of the photonic environment are encoded in the pﬃﬃﬃ iφ T free photonic Hamiltonian creation operator. The spinors jψ i¼ð1; τe Þ = 2, where qτ ph y denotes the transpose, describe an emergent pseudospin degree H ¼ h ω c c ph qn qn qn ð3Þ of freedom where the two components encode the relative qn amplitude and phase of the dipolar oscillations on the two inequivalent A and B sublattices, respectively, with φ = arg(f ). q q and in the light–matter interaction Hamiltonian ð1Þ ð2Þ These spinors can be represented by a pseudospin vector on the H ¼ H þ H , with int int int Bloch sphere which reads S = τ(cosφ , sinφ , 0). qτ q q ð1Þ y y y At the high symmetry K and K′ points (see Fig. 2a), the H ¼ h iξ ϕ a c þ a c qn qn int n q q qn qn sublattices decouple with no well-deﬁned relative phase (i.e., f = 0), ð4Þ 4π y y y pﬃﬃ giving rise to two inequivalent CDPs located at ± K ¼ ± ð ; 0Þ þh iξ ϕ b c þ b c þ H:c: qn qn 3 3a n q q qn qn as observed in Fig. 2b. These CDPs correspond to vortices in the pseudospin vector ﬁeld S , which give rise to topological qτ and singularities in the Berry curvature . Therefore, the CDPs are ð2Þ 2ξ ξ qn qn′ sources of quantized Berry wπ, where w = ±1 is the topological H ¼ h Re ϕ ϕ int ω n n′ charge of the Dirac point corresponding to the winding qnn′ ð5Þ number of the vortex. As expected from the symmetry of the y y c c þ c c þ H:c: metasurface, the existence of the CDPs is robust against long- qn qn′ qn qn′ range Coulomb interactions as shown in Supplementary Note 1. 1=2 In fact, for small cavity heights, the image dipoles quench long- where ξ / L parametrizes the strength of the light–matter qn range Coulomb interactions and the nearest-neighbor approx- interaction (see Methods for analytical expression). The bosonic imation becomes increasingly accurate as shown in Supplemen- operator c creates a TEM photon with wavevector q in the ﬁrst qn ph tary Figure 1. Brillouin zone and dispersion ω ¼ cjj q G ,where n indexes the qn n NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications 3 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 K' b b 1 2 Subcritical phase L > L Critical phase L = L Supercritical phase L < L c c c 1.05 b c d e Full polariton Quasistatic Hamiltonian dispersion Type-I CDP Two-band 1 Hamiltonian Quadratic band-degeneracy 1.00 3 4 Type-I CDP TEM mode 0.95 dispersion Type-II SDP Light-cone 0.90 Γ K(K') M Γ Γ K(K') M Γ Γ K(K') M Γ Γ K(K') M Γ Fig. 2 Evolution of the polariton dispersion as the cavity height is reduced. a First Brillouin zone deﬁned by primitive reciprocal lattice vectors b ¼ pﬃﬃﬃ pﬃﬃﬃ 2π 2π 3; 1 and b ¼ 3; 1 . b Quasistatic dispersion of the collective-dipole normal modes, where the upper band corresponds to a bright, 3a 3a symmetric dipole conﬁguration (↑↑) and the lower band corresponds to a dark, antisymmetric dipole conﬁguration (↑↓). 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 (solid black lines) and the two-band Hamiltonian H (orange dashed lines), for c pol mat subcritical (L = 5a), d critical (L= L = 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 , 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 ph parameters ω ¼ 2:5ω and Ω = 0.01ω K0 0 set of reciprocal lattice vectors G . The complex phase factors ϕ ¼ mediated by the electromagnetic ﬁeld. However, by expressing the expðÞ iaG ^y are associated with Umklapp processes that arise due interaction Hamiltonian (Eq. (4)) in terms of the β and β qτ n qτ to the discrete, in-plane translational symmetry of the metasurface, operators that diagonalize the matter Hamiltonian, X X and must be retained as they are critical for maintaining the point- ð1Þ y y y H ¼ h iΛ β c þ β c þ H:c:; int qnτ qτ qn qτ qn group symmetry of the polariton Hamiltonian. ð6Þ τ¼ ± qn We diagonalize H using a generalized Hopﬁeld–Bogoliubov pol transformation (see Methods for details), and in Fig. 2c–e, we ﬁnd that complex non-local interactions, which arise from we present the resulting polariton dispersion for different cavity strong multiple-scattering in the bipartite structure, result in a heights. Also, in Supplementary Figure 2, we present the full non-trivial winding of the light–matter coupling as a function of polariton dispersion which includes long-range Coulomb inter- wavevector direction actions. For small cavity heights, the full polariton dispersion iφ is almost indistinguishable from that obtained in the nearest- Λ / ξ ϕ e þ τϕ : ð7Þ qnτ qn n n neighbor approximation, and therefore one can conclude that long-range Coulomb interactions do not qualitatively affect Naively, one may expect all of the band crossings in Fig. 2bto the physics presented here. It is important to stress that our be avoided as a result of the hybridization between the collective- general model captures the essential physics that will emerge in a dipole and photonic modes, as it is a characteristic feature of variety of different experimental setups. To show this, in 48,49 polaritonic systems . Indeed, this is the case for the crossings iφ Supplementary Figure 3 and Supplementary Figure 4 we present with the upper quasistatic band where Λ / e þ 1 (see q0þ the polariton dispersions obtained from ﬁnite element simula- red line in Fig. 3a) due to the constructive interference between tions of a honeycomb array of plasmonic nanorods and the sublattices of this bright (↑↑) conﬁguration (see Fig. 3b, c). microwave helical resonators, respectively. Indeed, these entirely This results in a large anticrossing for all wavevector directions, different physical realizations show the same evolution of the as observed in Fig. 2c. In stark contrast, for the lower polariton spectrum as presented in Fig. 2c–e and Supplementary quasistatic band the coupling constant is signiﬁcantly reduced Figure 2. iφ Λ / e 1 (see blue line in Fig. 3a) due to the destructive q0 interference between the sublattices of this dark (↑↓) conﬁgura- Emergence of type-II Dirac points. Given the elementary nature tion (see Fig. 3e). Consequently, this results in a small of the individual resonant elements, one may be tempted to anticrossing for a general wavevector direction. assume that nothing peculiar could emerge from the ordinary Crucially, however, the light–matter interaction for the dipole–dipole interactions between the meta-atoms which are lower quasistatic band completely vanishes (Λ = 0) along the q0− 4 NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications | | | q 0 NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 ARTICLE 2.05 2 1 2.00 Inversion of chirality 0.10 Γ K(K') 4 4 0.05 v/v Γ→M Γ→K(K') Γ→M 1.00 Wavevector direction 1–D/ –1 bc d e t/t 0 0 0.92 0.2 1.0 /L 1 2 3 4 –2 0.2 0.4 0.6 0.8 1.0 /L Fig. 3 Non-trivial winding in the light–matter interaction. a Dependence of iφ the magnitude of the light–matter coupling constant jΛ j/je þ τj on q0τ Fig. 4 Tunable parameters in the effective Hamiltonian. Dependence of the the direction of q from Γ−Mto Γ−K(K′) (see inset), for the upper (red parameters in the effective polariton Hamiltonian (Eq. (9)) on the inverse line) and lower (blue line) quasistatic bands. Plots obtained with |q| = |K|/ cavity height. The blue dashed line shows the variation of the group velocity v 2. b–e Schematics of the bright (↑↑) and dark (↑↓) conﬁgurations of the two which changes sign at the critical cavity height L , leading to the inversion of sublattices interacting with the photonic mode which is indicated by the chirality. The orange dot-dashed line shows the variation of the trigonal ﬁeld proﬁle. Panels b and d represent the crossings labeled ‘1’ and ‘3’ in warping parameter t which becomes dominant close to criticality. These Fig. 2 along the Γ−K(K′) directions, respectively, while c and e represent parameters have been normalized to v= 3Ωa/2 and t= 3Ωa /8 which are the crossings labeled ‘2’ and ‘4’ along the Γ−M directions, respectively. the group velocity and trigonal warping parameters, respectively, in the Crucially, the light–matter interaction strength for the dark mode vanishes absence of image dipoles and light–matter interactions. The orange dot- (|Λ | = 0) along the Γ−K(K′) directions due to the complete destructive q0− dashed line in the inset shows the variation of the CDP frequency ω , while interference between the two sublattices (see d), leading to the emergence the blue dashed line in the inset shows the variation of the wavevector- of six inequivalent type-II Dirac points in the polariton spectrum ph dependent diagonal term D. Plots obtained with ω ¼ 2:5ω and Ω= 0.01ω K0 0 high-symmetry Γ−K(K′) directions, where φ = 0, due to the Similarly, the effective Hamiltonian near the K′ point is given by eff eff complete destructive interference between the two sublattices (see H ¼ðH Þ . In Eq. (9), the resonant frequency ω , group K′;k K;k 0 Fig. 3d). As a result, along these high-symmetry directions the velocity v, and trigonal warping parameter t, now encode non- crossings are protected, leading to six inequivalent Dirac points trivial contributions from the hybridization with the photonic emerging in the polariton spectrum—we call these satellite Dirac environment. There is also an additional wavevector-dependent points (SDPs) to distinguish them from the CDPs. As we will see diagonal term parametrized by D, which breaks the symmetry below, these SDPs belong to the type-II class of 2D Dirac points between the upper and lower polariton bands. The dependence of where the dispersion takes the form of a critically tilted Dirac these parameters on the cavity height is shown in Fig. 4 (see cone (see inset of Fig. 2c), characterized by open, hyperbolic Methods for analytical expressions). To leading order in k, one isofrequency contours. can observe that the effective Hamiltonian (Eq. (9)) near the CDP is equivalent to a 2D massless Dirac Hamiltonian. Therefore, the Effective Hamiltonian in the matter subspace. To explore the polariton CDPs remain in the type-I class and are robust against the coupling with the photonic environment—this is not sur- nature of the polaritons in the vicinity of the different Dirac points, we ﬁrst neglect non-resonant terms in the matter prising given that their physical origin is intrinsically linked to the lattice structure alone, which is preserved here. Hamiltonian and perform a unitary Schrieffer–Wolff transfor- mation on H to integrate out the photonic degrees of freedom To elucidate the nature of the SDPs, we expand the effective pol Hamiltonian (Eq. (9)) near one of the SDPs located at (see Methods for details). Finally, we extract the two-band Hamiltonian in the matter sublattice space K ¼ðÞ v=t; 0 and obtain ph Dv 2Dv X eff ξ ω qn qn H ¼ h ω k′ 1 þ hσ v k′ ; ð10Þ K ;k′ 0 x 2 H ¼ H 2h 2 mat mat 2 t ph qn ω ω qn ð8Þ where k′ measures the deviation from K and v ¼ v is y y 2 y 2 y a a þ b b þ ϕ b a þ ϕ a b : q q q q n q q n q q the velocity tensor. Apart from a global energy shift, the effective Hamiltonian (Eq. (10)) near the SDP takes the form Diagonalizing the two-band Hamiltonian (Eq. (8)) leads to an of a generalized 2D massless Dirac Hamiltonian P P effective dispersion (see Methods) which provides an excellent H ¼ hu k 1 þ hv k σ . If the parameters u and k i¼x;y i i 2 i¼x;y i i i description of the polaritons as indicated by the orange dashed 2 2 2 2 v satisfy the condition u =v þ u =v <1, then the Dirac cone x x y y lines in Fig. 2c–e. Finally, we expand the two-band Hamiltonian becomes tilted and anisotropic but still belongs to the type-I (Eq. (8)) up to quadratic order in k and obtain the effective eff eff y class with closed isofrequency contours. However, the condition Hamiltonian near the K point H ¼ ψ H ψ (see Sup- K k K;k k 2 2 2 2 u =v þ u =v >1deﬁnes a distinct type-II class of 2D Dirac x x y y plementary Note 2 for derivation) with Bloch Hamiltonian points, giving rise to a critically tilted Dirac cone with open, hyperbolic isofrequency contours. Hence, the type-I and type-II eff 2 2 ð9Þ H ¼ h ω 1 hvσ k þ htðÞ σ k hDjj k 1 : K;k 0 2 2 classes are related via a Lifshitz transition in the topology of the NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications 5 | | | |e +| / k 0 High Low ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 a b c Subcritical phase L >L Critical phase L =L Supercritical phase L < L c c c 1 1 –π Γ–K –π K–M –π +π –π +π +π - – ππ 0.5 0.5 0.5 –π –π –π 0 +π 0 –2π 0 +π –0.5 –0.5 –0.5 –π –1 –1 –1 –1 –0.5 0 0.5 1 –1 –0.5 0 0.5 1 –1 –0.5 0 0.5 1 k k k x x x d e f k k x y k k x y k k x y Massless Dirac polaritons Massless Dirac polaritons Massive chiral polaritons with inverted chirality Fig. 5 Merging of type-I and type-II Dirac points with chirality inversion. a–c Pseudospin vector ﬁeld and isofrequency contours (see Methods for the speciﬁc isofrequency values) for the upper polariton band near the K point for a subcritical (L = 2.3a), b critical (L= L = 1.78a), and c supercritical (L= 1.4a) cavity heights, respectively, as obtained from H . The Dirac points (corresponding to vortices) are depicted by orange circles along with their mat associated Berry ﬂux. Before criticality, three type-II SDPs (−π Berry ﬂux) are driven towards the type-I CDP (π Berry ﬂux) along the Γ−K directions as the cavity height is reduced (see inset in a). At the critical cavity height L , they merge together forming a quadratic band-degeneracy with combined Berry ﬂux 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 ﬂux, leaving only the type-I CDP remaining in the spectrum with π Berry ﬂux (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 deﬁned 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 ph and Berry phase of π, but with inverted chirality. All pseudospin and contour plots are obtained with ω ¼ 2:5ω and Ω = 0.01ω K0 0 0 isofrequency contours. Indeed, since we have u = 0 and of the two-band Hamiltonian (Eq. (8)). In Fig. 5a–cweplot the 2 2 2 u =v ¼ 4D =t >1, the SDPs belong to the type-II class of 2D pseudospin vector ﬁeld near the K point for different cavity heights x x Dirac points. Furthermore, since the Hamiltonian (Eq. (10)) is and schematically depict the location of the Dirac points, along with expressed in terms of σ , the pseudospin winds in the opposite their associated Berry ﬂux. Finally, in Fig. 5d–f, we illustrate the direction around the SDPs as compared to the CDP, and corresponding effective polariton spectrum to leading order in k. therefore the SDPs located along the Γ−K directions are sources Note that similaranalysiscan be performednearthe K′ point. of −π Berry ﬂux. As required by time-reversal symmetry, the In the subcritical phase (L > L ), three type-II SDPs are located SDPs located along the Γ−K′ directions are sources of π Berry along the Γ−K directions, each with −π Berry ﬂux surrounding ﬂux (opposite to the CDP located at the K′ point). a type-I CDP with π Berry ﬂux (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 Manipulation of type-I and type-II Dirac points. Wehavethus velocity v that is tunable via the cavity height (see Fig. 4). demonstrated that the honeycomb metasurface simultaneously Here, the effective Hamiltonian (Eq. (9)) is equivalent to a exhibits two distinct species of massless Dirac polaritons, namely 2D massless Dirac Hamiltonian with spinor eigenstates pﬃﬃﬃ iθ type-I and type-II. In contrast to the type-I CDPs, the existence of jψ i¼ð1; τe Þ = 2. These represent massless Dirac polar- kτ the type-II SDPs is intrinsically linked to the hybridization between itons with chirality hψ jσ kjψ i¼τ, resulting in a pseudos- kτ kτ the light and matter degrees of freedom, and thus one can pin that winds once around the CDP and a topological Berry manipulate their location within the Brillouin zone by simply phase of π (see Fig. 5d). modifying the light–matter interaction via the cavity height. As a At the critical cavity height (L = L ), the group velocity of the result, the polariton spectrum evolves into qualitatively distinct massless Dirac polaritons vanishes vðL Þ¼ 0 (see Fig. 4) as the phases as highlighted in Fig. 2c–e. To elucidate the differences type-II SDPs merge with the type-I CDP, forming a quadratic between these phases, we study the spinor eigenstates (see Methods) band-degeneracy (see Fig. 2d) with combined −2π Berry ﬂux (see 6 NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications | | | y NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 ARTICLE Fig. 5b). The leading order term in the effective Hamiltonian As a ﬁnal remark, we brieﬂy comment on how one might (Eq. (9)) is now quadratic in k with corresponding spinor probe the Dirac physics presented in this work. Given that the pﬃﬃﬃ i2θ eigenstates jψ i¼ð1; τe Þ = 2. Therefore, during this Dirac points exist in a polaritonic excitation spectrum, one kτ critical merging transition, the massless Dirac polaritons morph must drive the system with photons at the required frequency in into massive chiral polaritons with qualitatively distinct physical order to probe them. In fact, both of the type-I and type-II Dirac properties. These include a parabolic spectrum and chirality points lie outside of the light-line and therefore one must hψ jðσ kÞ jψ i¼τ, resulting in a pseudospin that winds overcome the momentum mismatch with photons. The speciﬁc kτ kτ twice as fast compared to the subcritical phase and a topological experimental technique that one would employ will depend on Berry phase of −2π (see Fig. 5e). the nature of the metasurface and the corresponding frequency Since the point-group symmetry is preserved, the type-II SDPs regime. For example, techniques for plasmonic systems have do not annihilate the type-I CDP, but they re-emerge and traditionally involved coupling via evanescent waves with prisms, continue to migrate along the K−M directions as the cavity gratings, and local scatterers , or more recent techniques such as height is reduced past criticality (L < L ) (see inset of Fig. 5c). non-linear wave-mixing . In contrast, realizations in the After a small decrease in cavity height, these SDPs annihilate with microwave regime can be probed using point-like antenna other SDPs that migrate along the opposite direction and have sources and detectors . In fact, microwave metamaterials are opposite Berry ﬂux. This topological transition leaves only the proving to be a versatile platform for exploring Dirac/Weyl type-I CDP remaining in the polariton spectrum with π Berry ﬂux physics, as one can directly probe the ﬁeld distributions using (see Figs. 2e and 5c). near-ﬁeld scanning techniques , and thus one could directly In this supercritical phase, we recover the linear dispersion probe the environment-induced chirality inversion predicted near the type-I CDP to leading order in k (see Fig. 2e), and here. the effective Hamiltonian (Eq. (9)) is equivalent to a 2D massless Dirac Hamiltonian with corresponding spinor eigen- pﬃﬃﬃ Discussion iθ states jψ i¼ð1; τe Þ = 2. Remarkably, massless Dirac polar- kτ To conclude, we have revealed rich and unique Dirac physics that itons thus re-emerge past criticality with an environment-induced emerges even in the most elementary honeycomb metasurfaces. inversion of chirality hψ jσ kjψ i¼ τ (see Fig. 5f). Physically, kτ kτ In particular, we have unveiled the simultaneous existence of both this corresponds to a π-rotation in the relative phase between the type-I and type-II massless Dirac polaritons, where the latter dipole oscillations on the two inequivalent sublattices, which is emerge from a non-trivial winding in the light–matter interac- also accompanied by a π-rotation in the isofrequency domains tion. We would like to emphasize that it is this unique physical (compare Fig. 5a and Fig. 5c). origin of the type-II SDPs, together with the truly 2D nature of We emphasize that it is the chirality of massless Dirac fermions the metasurface, that enables one to qualitatively modify the that is responsible for most of the remarkable properties of fundamental properties of these emergent Dirac polaritons by monolayer graphene, including the Klein tunneling phenomenon . manipulating the surrounding photonic environment alone. This Consequently, this unique environment-induced inversion of stands in stark contrast to conventional artiﬁcial graphene sys- chirality could give rise to unconventional phenomena such as tems where the fundamental properties are dictated by the lattice anomalous Klein transport. For example, near the K point, the structure. Therefore, exploiting the rich tunability of the polariton right-propagating polaritons correspond to an antisymmetric dipole spectrum with the environment offers a new paradigm that opens pﬃﬃﬃ conﬁguration jψ i¼ð1; 1Þ = 2 in the subcritical phase and to a variety of opportunities to explore unique Dirac-related physics kτ pﬃﬃﬃ a symmetric conﬁguration jψ i¼ð1; 1Þ = 2 in the supercritical at the subwavelength scale. kτ phase. Thus, due to the orthogonality between these two For example, one can simultaneously probe the dynamics spinor eigenstates, the inversion of chirality removes the channel of type-I and type-II massless Dirac quasiparticles, where responsible for the perfect transmission in the conventional Klein the latter are predicted to exhibit intriguing anomalous refraction 34,35 tunneling effect (see Fig. 5d, f). Such a scenario could be realized behavior . Furthermore, the environment-induced redshift in a simple setup characterized by two neighboring regions with of the CDP frequency ω (see Fig. 4) will allow the investigation different cavity heights. of polaritonic Klein tunneling through interfaces separating As a side remark, we note that the polariton spectrum near regions with different cavity heights. Consequently, negative criticality bears some resemblance with the low-energy spectrum refraction can be induced by simple variations in the cavity of bilayer graphene with its central Dirac point and three satellite height, which could be exploited in novel schemes for guiding and 52,53 Dirac points, which all belong to the type-I class . However, manipulating light at the subwavelength scale, including polari- 64,65 given the type-II nature of the polariton SDPs, the topology of the tonic Veselago lensing . Moreover, the tunable group velocity polariton isofrequency contours are markedly different from that will enable the exploration of velocity barriers for the unprece- of the bilayer spectrum. This is further highlighted at criticality dented guiding and localization of massless Dirac 66,67 where the polariton bands have the same curvature, which is in quasiparticles , which is extremely difﬁcult to achieve in real stark contrast to the electronic bands in bilayer graphene. graphene. One could also combine the effects of the environment We also note that recent works explored the possibility to with inhomogeneous strain deformations, giving rise to unique manipulate the (3+1) type-I Dirac points in bilayer graphene pseudomagnetic-related effects, including the intriguing ability to 54–57 through the application of lattice deformations , leading to induce a pseudo-Landau level spectrum for polaritons that can be the merging and annihilation of pairs of Dirac points. In addition, qualitatively tuned via the cavity height. Finally, the ability to a multimerging transition of all (3+1) type-I Dirac points controllably invert the chirality of the massless Dirac polaritons has been proposed theoretically within tight-binding models opens new perspectives for anomalous pseudorelativistic trans- involving the artiﬁcial tuning of third-nearest-neighbor hopping port through interfaces separating regions in distinct polaritonic 58–61 amplitudes in a graphene-like honeycomb structure . How- phases. ever, these proposals have no physical realization so far. In stark contrast, the imminently realizable metasurfaces in our work Methods enable the exploration of rich Dirac phases with ease by simply Derivation of the polaritonic Hamiltonian. The cavity-embedded metasurface modifying the photonic environment via an enclosing cavity. is composed of a honeycomb array of identical meta-atoms located at L L R ¼ R þ a^y þ ^z and R ¼ R a^y þ ^z on the A and B sublattices, respectively. A B 2 2 NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications 7 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 the local and block-diagonal form Here, R = l a + l a is an in-plane lattice translation vector with primitive vectors 1 1 2 2 a and a (see Fig. 1) and integers l and l . Each meta-atom is modeled by a single 1 2 1 2 y y H ¼ h ðÞ ω ΩS a a þ b b mat 0 q q q q dynamical degree of freedom h (with dimensions of length), where the electric- dipole moment associated with its fundamental eigenmode is p ¼Qh^z, with hi y y effective charge Q. The Coulomb potential energy between two dipole moments þh Ωð1 IÞ f b a þ a þ H:c: q q q q ð17Þ p and p′ located at generic positions r and r′, respectively, is given by 1 y y y y h ΩS a a þ b b þ H:c: : 2 q q q q p p′ 3ðp n ^Þðp′ n ^Þ V ¼ ; ð11Þ Coul 4πε jr r′j y y In the main text, we do not present the non-resonant terms (e.g., b a ), leading to q q Eq. (1) where ω ~ ¼ ω ΩS and Ω ¼ ΩðÞ 1 I . 0 0 In the Coulomb gauge, the light–matter interaction is described by the minimal- where n ^ ¼ðÞ r r′ =jj r r′ and ε is the vacuum permittivity. coupling Hamiltonian which, within the dipole approximation, reads The presence of the perfectly conducting metallic plates, placed at z = 0 and z = L, modiﬁes the boundary conditions on the scalar potential and, therefore, the X X X X Q Q Coulomb interaction between the meta-atoms. Using the method of images to H ¼ Π A ðR Þ þ A ðR Þ ; int R z s z s M 2M ensure the vanishing of the scalar potential at the cavity walls , we introduce an s¼A;B R s¼A;B R s s ð18Þ |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} inﬁnite series of image dipoles located outside the cavity at positions R þ lL^z, s ð1Þ ð2Þ H H int int 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 where we have used Π ¼ Π ^z. The vector potential can be decomposed into Eq. (11) , the matter Hamiltonian within the nearest-neighbor approximation R R transverse electric (TE) and transverse magnetic (TM) modes of the cavity. reads 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 P P s M 2 2 H ¼ þ ω h vector potential mat 2M 2 0 R s¼A;B R sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ h jj q G mπ 2 PP A ðr; zÞ¼ cos z þ h h 3 z R R þe ph mπ 4πε a 0 B B j ^ L q G þ z ε N NALω n R j¼1 qmn qmn L ð19Þ B 0 m hi þ1 P P P iðqG Þr y iðqG Þr 2 3 n n Q a 2 c e þ c e ; qmn qmn ′2 h 8πε a lL R s ð12Þ s¼A;B R l¼1 pﬃﬃﬃ 3 þ1 2 lL 2 PP P 2jj 1 where A¼ 3 3a =2 is the area of a unit cell and N = 1 + δ . The bosonic Q a m m0 ′ h h 3 5 y R R þe 8πε a 0 2 B B j operator c creates a TM photon with wavevector q in the ﬁrst Brillouin zone lL qmn R j¼1 l¼1 B 1þjj ph and dispersion ω ¼ cjj q G þ zmπ=L . Here, G = n b + n b is a reciprocal n n 1 1 2 2 qmn 3 þ1 lL 2 PP P 2jj 1 Q a lattice vector with primitive vectors b and b , where n indexes the set of ordered 1 2 ′ h h 3 5 R R e 8πε a A A j 0 2 2 lL pairs of integers (n , n ), and m is a non-negative integer indexing the different TM R j¼1 l¼1 1þ 1 2 A jj cavity modes. Only TM photons with even m couple to the dipoles due to the parity selection rule at the center of the cavity. where the primed summations exclude the l= 0 term. Here, Π is the conjugate Substituting the vector potential (Eq. (19)) into Eq. (18) we obtain the light–matter interaction Hamiltonian given by momentum to the dynamical coordinate h corresponding to the meta-atom located at R ,and M is an effective mass. Next, we introduce the bosonic operators s P ð1Þ y y y H ¼ h iξ ϕ a c þ a c int qmn n q qmn q qmn qmn sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð20Þ Mω 1 P y y y a ¼ h þ i Π ð13Þ þh iξ ϕ b c þ b c þ H:c: R R R qmn n q qmn q qmn A A A 2h 2hMω 0 qmn and and ð2Þ 2ξ ξ qmn qm′n′ H ¼ h Re ϕ ϕ int n n′ sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ 0 rﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ qmm′nn′ ð21Þ Mω 1 b ¼ h þ i Π ð14Þ y y y R R R B B B c c þ c c þ H:c: 2h 2hMω qmn qm′n′ qmn qm′n′ The strength of the light–matter interaction is parametrized by that annihilate quanta of the fundamental eigenmode on the meta-atom located at !1 ph R and R , respectively, and satisfy the commutation relations ½a ; a ¼ δ , A B ω R R′ RR′ 8π a Ω q0n ph y y ð22Þ ξ ¼ ω Fðω Þ pﬃﬃﬃ ; ½b ; b ¼ δ ,and ½a ; b ¼ 0. In terms of these operators, the matter qmn 0 qmn R RR′ R ph ph R′ R′ 3 3N ω ω qmn m qmn Hamiltonian (Eq. (12)) reads P P where, to take into account the ﬁnite size of the meta-atoms, we have introduced a y y H ¼ hω a a þ hω b b ph mat 0 R R 0 R R A A B B phenomenological function Fðω Þ that provides a smooth cut-off for the qmn R R A B interaction with short-wavelength photonic modes where the dipole approximation h i PP y y breaks down. We choose the phenomenological cut-off function to be of the Fermi- þhΩð1 IÞ b a þa þ H:c: R R þe R þe B B j B j R j¼1 Dirac distribution form ð15Þ hi y y h 1 ph ΩS a a þ a þ H:c: 2 R R R A A A Fðω Þ¼ ; ð23Þ ph qmn 2ðω 3ω Þ=ω A qmn 0 0 1 þ e hi y y ΩS b b þ b þ H:c: 2 R R R B B B R which is smooth enough to avoid spurious artifacts appearing in the polariton spectrum. Finally, the free photonic Hamiltonian of the cavity reads ph y 2 3 H ¼ h ω c c : where Ω= Q /8πε Mω a parametrizes the strength of the nearest-neighbor ph qmn 0 0 qmn qmn ð24Þ qmn Coulomb interaction, and the parameters 1 1 lL X X a 2 1 In Eqs. (3), (4), (5) in the main text, we only present the contribution from the S¼ 4 ; I¼ 2 hi5 ð16Þ lL 2 lL TEM mode (m = 0), dropping the corresponding index. In Supplementary Note 1, l¼1 l¼1 1 þ we discuss the effect of the higher order (m ≠ 0) TM cavity modes for larger cavities. encode renormalizations due to the cavity-induced image dipoles. We apply Born- von Kármán boundary conditions over a lattice withN 1 unit cells and introduce Hopﬁeld–Bogoliubov diagonalization. The polariton Hamiltonian H = H + pol mat 1=2 iqR the Fourier transform of the bosonic operators a ¼N a e and H + H , where H is given by Eq. (17), H by Eq. (24), and H by Eqs. (20) ph int mat ph int R q q P A P 1=2 y pol iqR 1 b ¼N b e , which transforms the matter Hamiltonian (Eq. (15)) into and (21), can be recast into matrix form as H ¼ Ψ H Ψ where R q q pol q q q q 8 NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 ARTICLE y y y T T y y y Ψ ¼ðψ ; C ; ψ ; C Þ. Here, ψ ¼ða ; b Þ is the spinor creation operator in Hamiltonian (Eq. (32)) reads q q q q q q q q y y y y the matter sublattice space and C ¼ðc ; c ; ¼ ; c ; ¼ ; c Þ is the vector of TM qp q q1 q2 qN hi ð1Þ ð2Þ photon creation operators, where p indexes the set of ordered triplets of integers H ’ H þ S; H þ H þ H ; pol mat int ph int |ﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄ} ð35Þ (n , n , m), and N is the total number of photonic operators considered. The 1 2 |ﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ{zﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄﬄ} pol Hermitian [2(N + 2)] × [2(N + 2)] matrix H can be written in block form as ph q mat 0 1 H H W q q pol where the matter and photonic subspaces are decoupled to quadratic order in ξ . @ A qmn H ¼ y ; ð25Þ Calculating the commutator in Eq. (35) and extracting the Hamiltonian within the H W H q q q matter sublattice space, we obtain the two-band Hamiltonian 2 ph where P ξ ω qmn qmn y y 2 y 2 y H ¼ H 2h a a þb b þ ϕ b a þ ϕ a b : ph 2 ð36Þ mat mat 2 q q q q n q q n q q ðω Þ ω ~ qmn ph ph ph qmn 0 ph W ¼ hDiag ω ; ω ; ω ; ω ; ¼ ; ω ; ¼ ; ω ð26Þ q 0 0 q1 q2 qp qN is the (N + 2) × (N + 2) diagonal matrix of resonant frequencies of the free oscil- In Eq. (8) in the main text, we only present the contribution from the TEM lators. The (N + 2) × (N + 2) block matrices H can be sub-divided into block mode (m = 0), dropping the corresponding index. We can recast the Hamiltonian mat matrices (Eq. (36)) into matrix form as H ¼ ψ H ψ , with Bloch Hamiltonian mat q q q q 0 1 mat int H H ! q q W F @ A mat q q H ¼ y ; ð27Þ q H ¼ h : ð37Þ int ph ± H H F W q q q q P P where the upper-diagonal block ~ Here W ¼ ω ~ Ω Δ and F ¼ Ωf Ω Δ ϕ with q 0 mn qmn q q mn qmn n ω ~ Ωf 2 0 q ph mat H ¼ h ð28Þ 16π a ω q0n 0 2 ph ~ Δ ¼ pﬃﬃﬃ F ðω Þ: ð38Þ Ωf ω ~ qmn qmn ph 2 ph q 0 L 2 3 3N ðω Þ ω ~ ω m qmn qmn ph is the 2 × 2 matrix in the matter subspace, and the lower-diagonal block H is the mat Diagonalizing H leads to the two-band dispersion ω ¼ W þ τjF j, which is N × N matrix in the photonic subspace with components mat qτ q q indicated by the orange-dashed lines in Fig. 2c–e. The corresponding spinor pﬃﬃﬃ no ξ ξ iφ qp qp′ ph ph eigenstates jψ i¼ ð1; τe Þ = 2, where φ ¼ argðF Þ, can be represented by the qτ q q ð29Þ H ¼ hω δ þ 4h Re ϕ ϕ : pp′ q qp p p′ pseudospin vector S ¼ τðcos φ ; sin φ ; 0Þ from which we obtain the pseudospin pp′ ω qτ q q vector ﬁeld plots in Fig. 5a–c. int Finally, the off-diagonal block H in Eq. (27) is the 2 × N interaction matrix, where the pth column reads Expansion of the effective two-band Hamiltonian. Near the K point, the func- tion Δ , given by Eq. (38), expands as qmn iξ ϕ qp p int H ¼ h : ð30Þ hi ð0Þ ð1Þ p iξ ϕ 2 qp Δ ’ Δ a Δ ðK G Þ k þðK G Þ k kmn Kmn Kmn n x x n y y hi ð1Þ ð2Þ 1 2 4 2 þ a Δ þ a Δ ðÞ K G k Kmn Kmn n x x ð39Þ hi The polariton Hamiltonian H is diagonalized by a generalized pol ð1Þ ð2Þ 2 1 2 4 2 48 y y T þ a Δ þ a Δ ðÞ K G k Hopﬁeld–Bogoliubov transformation Ψ = T X , where X ¼ðχ ; χ Þ and n q q q 2 Kmn Kmn y y q q q y y y y y χ ¼ðγ ; γ ; ¼ ; γ ; ¼ ; γ Þ. To ensure the invariance of the bosonic q q1 q2 qν qNþ2 ð2Þ þa Δ ðÞ K G ðÞ K G k k commutation relations for the transformed operators, T must be a [2(N + 2)] × [2 Kmn n x n y x y 70 y y (N + 2)] paraunitary matrix that satisﬁes T η T ¼ T η T ¼ η , where η ¼ q z q q z q z z y ðυÞ σ 1 and σ is the Pauli matrix. The transformed bosonic operators γ ¼ z 2 z to quadratic order in k, where the real parameters Δ (υ = 0, 1, 2) depend only qν Kmn y ph Ψ η jΨ i and γ = ⟨Ψ |η Ψ that diagonalize the polariton Hamiltonian as qν qν z q on the photon frequencies ω at the K point. Collecting the contributions from q z qν Kmn the degenerate photons (see Supplementary Note 2 for details), we obtain the pol y H ¼ h ω γ γ ; pol qν qν qν ð31Þ effective Hamiltonian (Eq. (9)), where parameters are given by qν ω Ω Ω ð0Þ ¼ 1 S Δ ; ð40Þ create and annihilate polaritons in the vth band, respectively. The polariton Kmn ω ω ω 0 0 0 pol mn dispersion ω (black solid lines in Fig. 2c–e) and the corresponding linearly qν independent eigenvectors |Ψ ⟩ (ﬁrst two columns of T ) are determined from the qν q positive eigenvalue solutions to the non-Hermitian eigenvalue equation v 4π ð1Þ pol pol η H jΨ i¼ hω jΨ i. ¼ 1 I A Δ ; ð41Þ z q qν qν qν n Kmn v 27 mn Schrieffer–Wolff transformation. To obtain an effective two-band Hamiltonian in the matter sublattice space, we neglect non-resonant terms in the matter 2 t 8π ð2Þ ¼ 1 I B Δ ; ð42Þ Hamiltonian (since Ω=ω 1 for practical realizations of the metasurface), but n Kmn t 81 mn 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 1 X 2 S S D Ω 4π 1 ð2Þ ð1Þ H ¼ e H e ¼ H þ½S; H þ ½S½S; H þ ¼ ð32Þ pol pol pol pol pol ¼ C Δ Δ ; ð43Þ n Kmn Kmn 2 2 ω a ω 27 2 0 0 mn and impose the Schrieffer–Wolff condition hi with ð1Þ S; H þ H ¼H ð33Þ pﬃﬃﬃ mat ph int 3 4π A ¼ ðÞ 2 3n cos ðÞ n þ n n 1 1 2 2 3 which eliminates the light–matter interaction to ﬁrst order in ξ . From Eq. (33), qmn ð44Þ 1 4π the particular form of the anti-Hermitian operator S reads þðÞ 6n 3n sin ðÞ n þ n ; 2 1 1 2 P 2 3 iξ qmn y y S¼ ðϕ a þ ϕ b Þc ph qmn n q n q ω ω ~ qmn qmn ð34Þ iξ qmn y y y 4π þ ðϕ a þ ϕ b Þc H:c:; ph 2 2 n q n q qmn ω þω ~ qmn 0 B ¼ 3n 6n 6n þ 6n n þ 2 cos ðÞ n þ n qmn n 1 1 2 1 2 1 2 ð45Þ pﬃﬃﬃ 4π ph where we have used the approximation jω ± ω j Ωjf j that is valid near the K þ 3 2n 4n þ 6n n 3n sin ðÞ n þ n ; qmn 0 q 1 2 1 2 1 1 2 and K′ points. Retaining leading-order terms in ξ , the transformed polariton qmn NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications 9 | | | ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 and 16. Rechtsman, M. C. et al. Topological creation and destruction of edge states in photonic graphene. Phys. Rev. Lett. 111, 103901 (2013). C ¼ 1 þ 3nðÞ n 1þ 3nðÞ n n : ð46Þ n 1 1 2 2 1 17. Rechtsman, M. C. et al. Strain-induced pseudomagnetic ﬁeld and photonic Landau levels in dielectric structures. Nat. Photonics 7, 153–158 (2013). 18. Rechtsman, M. C. et al. Photonic Floquet topological insulators. Nature 496, For brevity, we retain only the dominant (m = 0) TEM contribution for the 196–200 (2013). ðυÞ plots in Fig. 4, where the coefﬁcients Δ in Eq. (39) are given by K0n 19. Ni, X., Purtseladze, D., Smirnova, D. A., Slobozhanyuk, A. & Alù, A. Spin and "# 2 valley polarized one-way Klein tunneling in photonic topological insulators. 4π a ω ð0Þ ph 0 2 Δ ¼ 2 pﬃﬃﬃ F ðω Þ; ð47Þ Preprint at https://arxiv.org/abs/1707.05802 (2017). K0n K0n ph 2 3 3 ðω Þ ω ~ K0n 0 20. Hasegawa, Y., Konno, R., Nakano, H. & Kohmoto, M. Zero modes of tight- binding electrons on the honeycomb lattice. Phys. Rev. B 74, 033413 (2006). 2 2 21. Wunsch, B., Guinea, F. & Sols, F. Dirac-point engineering and topological phase 1 ph ð1Þ ω ω 4π a K00 0 pﬃﬃ Δ ¼ 4 ph 2 transitions in honeycomb optical lattices. New J. Phys. 10, 103027 (2008). K0n 2 3 3 L ω 0 ðω Þ ω ~ K0n 0 ð48Þ 22. Pereira, V. M., Castro Neto, A. H. & Peres, N. M. R. Tight-binding approach ph 2 ph 2 ∂ F ðω Þ ~ ph ½ ðω Þ ω 2 K0n K0n 0 to uniaxial strain in graphene. Phys. Rev. B 80, 045401 (2009). F ðω Þ ; ph ph K0n ∂ω 2ω K0n K0n 23. Montambaux, G., Piéchon, F., Fuchs, J. N. & Goerbig, M. O. Merging of Dirac points in a two-dimensional crystal. Phys. Rev. B 80, 153412 (2009). and 24. Guinea,F., Katsnelson, M.I.&Geim, A.K.Energygaps and a zero-ﬁeld quantum 4 3 Hall effect in graphene by strain engineering. Nat. Phys. 6,30–33 (2009). ph 3 2 ð2Þ ω ω 4π a K00 0 pﬃﬃ Δ ¼ 16 25. Soluyanov, A. A. et al. Type-II Weyl semimetals. Nature 527, 495–498 (2015). ph K0n 2 2 L ω 3 3 0 ðω Þ ω ~ K0n 0 26. Deng, K. et al. Experimental observation of topological Fermi arcs in type-II 2 ph ∂ F ðω Þ ph ½ 2 K0n Weyl semimetal MoTe . Nat. Phys. 12, 1105–1110 (2016). F ðω Þ ph K0n ∂ω 27. Huang, L. et al. Spectroscopic evidence for a type II Weyl semimetallic state in K0n ð49Þ ph 2 2 ph 2 2 MoTe . Nat. Mater. 15, 1155–1160 (2016). 5ðω Þ ω ~ ðω Þ ω ~ 2 ½½ K0n 0 K0n 0 ph 3 28. Huang, H., Zhou, S. & Duan, W. Type-II Dirac fermions in the PtSe2. class of 8ðω Þ K0n 2 transition metal dichalcogenides. Phys. Rev. B 94, 121117(R) (2016). ph ph 2 2 2 2 ∂ F ðω Þ ðω Þ ω ~ ½½ K0n K0n 0 þ : 29. Yan, M. et al. Lorentz-violating type-II Dirac fermions in transition metal ph 2 ph 2 ð∂ω Þ 8ðω Þ K0n K0n dichalcogenide PtTe . Nat. Commun. 8, 257 (2017). 30. Xiao, M., Lin, Q. & Fan, S. Hyperbolic Weyl point in reciprocal chiral metamaterials. Phys. Rev. Lett. 117, 057401 (2016). Isofrequency contour plots. In ascending order, the normalized frequency values 31. Chen, W.-J., Xiao, M. & Chan, C. T. Photonic crystals possessing multiple corresponding to the isofrequency contours are 0.98300, 0.98800, 0.99022, and Weyl points and the experimental observation of robust surface states. Nat. 0.99080 for Fig. 5a, 0.98080, 0.98180, 0.98250, and 0.98300 for Fig. 5b, and ﬁnally Commun. 7, 13038 (2016). 0.97215, 0.97278, 0.97350, and 0.97395 for Fig. 5c. 32. Noh, J. et al. Experimental observation of optical Weyl points and Fermi arc- like surface states. Nat. Phys. 13, 611–618 (2017). 33. Yang, B. et al. Direct observation of topological surface-state arcs in photonic Data availability. All relevant data are available from the corresponding authors metamaterials. Nat. Commun. 8, 97 (2017). upon reasonable request. 34. Wang, H.-X., Chen, Y., Hang, Z. H., Kee, H.-Y. & Jiang, J.-H. Type-II Dirac photons. npj Quantum Mater. 2, 54 (2017). Received: 17 November 2017 Accepted: 27 March 2018 35. Lin, J. Y., Hu, N. C., Chen, Y. J., Lee, C. H. & Zhang, X. Line nodes, Dirac points, and Lifshitz transition in two-dimensional nonsymmorphic photonic crystals. Phys. Rev. B 96, 075438 (2017). 36. Pyrialakos, G. G., Nye, N. S., Kantartzis, N. V. & Christodoulides, D. N. References Emergence of type-II Dirac points in graphynelike photonic lattices. Phys. Rev. 1. Novoselov, K. S. et al. Electric ﬁeld effect in atomically thin carbon ﬁlms. Lett. 119, 113901 (2017). Science 306, 666–669 (2004). 37. Jacqmin, T. et al. Direct observation of Dirac cones and a ﬂatband in a 2. Polini, M., Guinea, F., Lewenstein, M., Manoharan, H. C. & Pellegrini, V. honeycomb lattice for polaritons. Phys. Rev. Lett. 112, 116402 (2014). Artiﬁcial honeycomb lattices for electrons, atoms and photons. Nat. 38. Nalitov, A. V., Solnyshkov, D. D. & Malpuech, G. Polariton Z topological Nanotechnol. 8, 625–633 (2013). insulator. Phys. Rev. Lett. 114, 116401 (2015). 3. Zhu, S.-L., Wang, B. & Duan, L.-M. Simulation and detection of Dirac fermions 39. Karzig, T., Bardyn, C. E., Lindner, N. H. & Refael, G. Topological polaritons. with cold atoms in an optical lattice. Phys. Rev. Lett. 98, 260402 (2007). Phys. Rev. X 5, 031001 (2015). 4. Peleg, O. et al. Conical diffraction and gap solitons in honeycomb photonic 40. Bardyn, C. E., Karzig, T., Refael, G. & Liew, T. C. H. Topological polaritons lattices. Phys. Rev. Lett. 98, 103901 (2007). and excitons in garden-variety systems. Phys. Rev. B 91, 161413(R) (2015). 5. Han, D., Lai, Y., Zi, J., Zhang, Z. Q. & Chan, C. T. Dirac spectra and edge 41. Yi, K. & Karzig, T. Topological polaritons from photonic Dirac cones coupled states in honeycomb plasmonic lattices. Phys. Rev. Lett. 102, 123904 (2009). to excitons in a magnetic ﬁeld. Phys. Rev. B 93, 104303 (2016). 6. Gibertini, M. et al. Engineering artiﬁcial graphene in a two-dimensional 42. Yuen-Zhou, J. et al. Plexciton Dirac points and topological modes. Nat. electron gas. Phys. Rev. B 79, 241406(R) (2009). Commun. 7, 11783 (2016). 7. Torrent, D. & Sánchez-Dehesa, J. Acoustic analogue of graphene: observation 43. Alù, A. First-principles homogenization theory for periodic metamaterials. of Dirac cones in acoustic surface waves. Phys. Rev. Lett. 108, 174301 (2012). Phys. Rev. B 84, 075153 (2011). 8. Gomes, K. K., Mar, W., Ko, W., Guinea, F. & Manoharan, H. C. Designer 44. Yves, S. et al. Crystalline metamaterials for topological properties at Dirac fermions and topological phases in molecular graphene. Nature 483, subwavelength scales. Nat. Commun. 8, 16023 (2017). 306–310 (2012). 45. Kaina, N., Lemoult, F., Fink, M. & Lerosey, G. Negative refractive index and 9. Bellec, M., Kuhl, U., Montambaux, G. & Mortessagne, F. Tight-binding acoustic superlens from multiple scattering in single negative metamaterials. couplings in microwave artiﬁcial graphene. Phys. Rev. B 88, 115437 (2013). Nature 525,77–81 (2015). 10. Weick, G., Woollacott, C., Barnes, W. L., Hess, O. & Mariani, E. Dirac-like 46. Craig, D. P. & Thirunamachandran, T. Molecular Quantum Electrodynamics: plasmons in honeycomb lattices of metallic nanoparticles. Phys. Rev. Lett. 110, An Introduction to Radiation-Molecule Interactions. (Academic Press, 106801 (2013). London, 1984). 11. Yu, S. et al. Surface phononic graphene. Nat. Mater. 15, 1243–1247 (2016). 47. Fuchs, J. N., Piéchon, F., Goerbig, M. O. & Montambaux, G. Topological Berry 12. Castro Neto, A. H., Peres, N. M. R., Novoselov, K. S., Geim, A. K. & Guinea, F. phase and semiclassical quantization of cyclotron orbits for two dimensional The electronic properties of graphene. Rev. Mod. Phys. 81, 109–162 (2009). electrons in coupled band models. Eur. Phys. J. B 77, 351–362 (2010). 13. Katsnelson, M. I., Novoselov, K. S. & Geim, A. K. Chiral tunnelling and the 48. Hopﬁeld, J. J. Theory of the contribution of excitons to the complex dielectric Klein paradox in graphene. Nat. Phys. 2, 620–625 (2006). constant of crystals. Phys. Rev. Lett. 112, 1555–1567 (1958). 14. Tarruell, L., Greif, D., Uehlinger, T., Jotzu, G. & Esslinger, T. Creating, moving 49. Low, T. et al. Polaritons in layered 2D materials. Nat. Mater. 16,182–194 (2017). and merging Dirac points with a Fermi gas in a tunable honeycomb lattice. 50. Schrieffer, J. R. & Wolff, P. A. Relation between the Anderson and Kondo Nature 483, 302–305 (2012). Hamiltonians. Phys. Rev. 149, 491–492 (1966). 15. Bellec, M., Kuhl, U., Montambaux, G. & Mortessagne, F. Topological 51. Goerbig, M. O., Fuchs, J. N., Montambaux, G. & Piéchon, F. Tilted anisotropic transition of Dirac points in a microwave experiment. Phys. Rev. Lett. 110, Dirac cones in quinoid-type graphene and α-(BEDT-TTF) I . Phys. Rev. B 78, 2 3 033902 (2013). 045415 (2008). 10 NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications | | | NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-03982-7 ARTICLE C.-R.M acknowledges the EPSRC Centre for Doctoral Training in Metamaterials (Grant 52. McCann, E. & Fal’ko, V. I. Landau-level degeneracy and quantum Hall effect No. EP/L015331/1) and QinetiQ for additional funding. G.W. acknowledges ﬁnancial in a graphite bilayer. Phys. Rev. Lett. 96, 086805 (2006). support from Agence Nationale de la Recherche (Project ANR-14-CE26-0005 Q-Meta- 53. McCann, E. & Koshino, M. The electronic properties of bilayer graphene. Rep. Mat) and the Centre National de la Recherche Scientiﬁque through the Projet Interna- Prog. Phys. 76, 056503 (2013). tional de Coopération Scientiﬁque program (Contract No. 6384 APAG). W.L.B. 54. Mucha-Kruczyński, M., Aleiner, I. L. & Fal’ko, V. I. Strained bilayer graphene: acknowledges the ﬁnancial support from EPSRC (Grant No. EP/K041150/1) and EU Band structure topology and Landau level spectrum. Phys. Rev. B 84, ERC project Photmat (ERC-2016-ADG-742222). E.M. acknowledges ﬁnancial support 041404(R) (2011). from the Leverhulme Trust (Research Project Grant RPG-2015-101), and the Royal 55. Mariani, E., Pearce, A. J. & Von Oppen, F. Fictitious gauge ﬁelds in bilayer Society (International Exchange Grant No. IE140367, Newton Mobility Grants graphene. Phys. Rev. B 86, 165448 (2012). NI160073, and Theo Murphy Award TM160190). 56. Verberck, B., Partoens, B., Peeters, F. M. & Trauzettel, B. Strain-induced band gaps in bilayer graphene. Phys. Rev. B 85, 125403 (2012). 57. Pearce, A. J., Cavaliere, F. & Mariani, E. Conductance and shot noise in Author contributions strained bilayer graphene. J. Phys. Condens. Matter 25, 375301 (2013). C.-R.M. performed the theoretical calculations and wrote the manuscript with E.M.; T.J. 58. Bena, C. & Simon, L. Dirac point metamorphosis from third-neighbor S. contributed to the theoretical calculations; E.M. and G.W. conceived the project; and E. couplings in graphene and related materials. Phys. Rev. B 83, 115404 (2011). M. and W.L.B. supervised the project. All authors commented on the manuscript. 59. Montambaux, G. An equivalence between monolayer and bilayer honeycomb lattices. Eur. Phys. J. B 85, 375 (2012). Additional information 60. Sticlet, D. & Piéchon, F. Distant-neighbor hopping in graphene and Haldane models. Phys. Rev. B 87, 115402 (2013). Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- 61. Bhattacharya, U., Hutchinson, J. & Dutta, A. Quenching in Chern insulators 018-03982-7. with satellite Dirac points: the fate of edge states. Phys. Rev. B 95, 144304 Competing interests: The authors declare no competing interests. (2017). 62. Maier, S. A. Plasmonics: Fundamentals and Applications. (Springer-Verlag, Reprints and permission information is available online at http://npg.nature.com/ Berlin, 2007). reprintsandpermissions/ 63. Renger, J., Quidant, R., Van Hulst, N., Palomba, S. & Novotny, L. Free-space excitation of propagating surface plasmon polaritons by nonlinear four-wave mixing. Phys. Rev. Lett. 103, 266802 (2009). Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional afﬁliations. 64. Cheianov, V. V., Fal’ko, V. & Altshuler, B. L. The focusing of electron ﬂow and a Veselago lens in graphene p-n junctions. Science 315, 1252–1255 (2007). 65. Lee, G.-H., Park, G. & Lee, H. Observation of negative refraction of Dirac fermions in graphene. Nat. Phys. 11, 925–929 (2015). Open Access This article is licensed under a Creative Commons 66. Raoux, A. et al. Velocity-modulation control of electron-wave propagation in graphene. Phys. Rev. B 81, 073407 (2010). Attribution 4.0 International License, which permits use, sharing, 67. Downing, C. A. & Portnoi, M. E. Localization of massless Dirac particles via adaptation, distribution and reproduction in any medium or format, as long as you give spatial modulations of the Fermi velocity. J. Phys. Condens. Matter 29, 315301 appropriate credit to the original author(s) and the source, provide a link to the Creative (2017). Commons license, and indicate if changes were made. The images or other third party 68. Jackson, J. D. Classical Electrodynamics, 3rd edn (Wiley: New York, 1999). material in this article are included in the article’s Creative Commons license, unless 69. Power, E. A. & Thirunamachandran, T. Quantum electrodynamics in a cavity. indicated otherwise in a credit line to the material. If material is not included in the Phys. Rev. A 25, 2473–2484 (1982). article’s Creative Commons license and your intended use is not permitted by statutory 70. Colpa, J. H. P. Diagonalization of the quadratic boson Hamiltonian. Physica A regulation or exceeds the permitted use, you will need to obtain permission directly from 93, 327–353 (1978). the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Acknowledgements We would like to thank Ian Hooper for useful discussions and help with the numerical © The Author(s) 2018 simulations. C.-R.M. and T.J.S. acknowledge ﬁnancial support from the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom. In addition, NATURE COMMUNICATIONS (2018) 9:2194 DOI: 10.1038/s41467-018-03982-7 www.nature.com/naturecommunications 11 | | |
Nature Communications – Springer Journals
Published: Jun 6, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.
Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.
All the latest content is available, no embargo periods.
“Hi guys, I cannot tell you how much I love this resource. Incredible. I really believe you've hit the nail on the head with this site in regards to solving the research-purchase issue.”
Daniel C.
“Whoa! It’s like Spotify but for academic articles.”
@Phil_Robichaud
“I must say, @deepdyve is a fabulous solution to the independent researcher's problem of #access to #information.”
@deepthiw
“My last article couldn't be possible without the platform @deepdyve that makes journal papers cheaper.”
@JoseServera
DeepDyve Freelancer | DeepDyve Pro | |
---|---|---|
Price | FREE | $49/month |
Save searches from | ||
Create lists to | ||
Export lists, citations | ||
Read DeepDyve articles | Abstract access only | Unlimited access to over |
20 pages / month | ||
PDF Discount | 20% off | |
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.
ok to continue