# Fileset

[sciadv.adp3487.pdf](https://mdr.nims.go.jp/filesets/881508e3-009a-4694-8318-627de5032e7c/download)

## Creator

Lukas Wehmeier, Suheng Xu, Rafael A. Mayer, Brian Vermilyea, Makoto Tsuneto, Michael Dapolito, Rui Pu, Zengyi Du, Xinzhong Chen, Wenjun Zheng, Ran Jing, Zijian Zhou, [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), Adrian Gozar, Qiang Li, Alexey B. Kuzmenko, G. Lawrence Carr, Xu Du, Michael M. Fogler, D. N. Basov, Mengkun Liu

## Rights

[Creative Commons BY-NC Attribution-NonCommercial 4.0 International](https://creativecommons.org/licenses/by-nc/4.0/)

## Other metadata

[Landau-phonon polaritons in Dirac heterostructures](https://mdr.nims.go.jp/datasets/358f8bd7-f9bb-43de-aa4b-043de581b14e)

## Fulltext

Landau-phonon polaritons in Dirac heterostructuresWehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e1 of 16P H Y S I C SLandau-phonon polaritons in Dirac heterostructuresLukas Wehmeier1,2*, Suheng Xu3, Rafael A. Mayer1, Brian Vermilyea4, Makoto Tsuneto1,  Michael Dapolito1,3, Rui Pu1, Zengyi Du1, Xinzhong Chen1,3, Wenjun Zheng1,  Ran Jing1,5, Zijian Zhou1, Kenji Watanabe6, Takashi Taniguchi7, Adrian Gozar8,9,10,  Qiang Li1,5, Alexey B. Kuzmenko11, G. Lawrence Carr2, Xu Du1, Michael M. Fogler4*,  D. N. Basov3*, Mengkun Liu1,2*Polaritons are light-matter quasiparticles that govern the optical response of quantum materials at the nanoscale, enabling on-chip communication and local sensing. Here, we report Landau-phonon polaritons (LPPs) in magne-tized charge-neutral graphene encapsulated in hexagonal boron nitride (hBN). These quasiparticles emerge from the interaction of Dirac magnetoexciton modes in graphene with the hyperbolic phonon polariton modes in hBN. Using infrared magneto-nanoscopy, we reveal the ability to completely halt the LPP propagation in real space at quantized magnetic fields, defying the conventional optical selection rules. The LPP-based nanoscopy also tells apart two fundamental many-body phenomena: the Fermi velocity renormalization and field-dependent magne-toexciton binding energies. Our results highlight the potential of magnetically tuned Dirac heterostructures for precise nanoscale control and sensing of light-matter interaction.INTRODUCTIONPolaritons are light-matter quasiparticles that play a fundamental role in the optical response of polarizable materials (1–18). Phonon polaritons were studied historically first (18), and they are examples of modes demonstrating strong light-matter coupling. In complex materials, polaritons can involve several distinct matter excitations, yielding a rich variety of collective phenomena (3, 19–21). If the optical properties of a material are tunable, polaritons inherit this tunability. For example, the dispersion of plasmon-polaritons in two-dimensional (2D) conductors can be controlled by changing their charge carrier concentration (13–15) or applying an electric current (16, 17). However, attaining strong mode coupling with conducting materials is difficult because of their high electronic losses. Graphene is one of the promising polaritonic platforms be-cause of its low intrinsic electron scattering rate (22) and corre-sponding high quality factors (3, 8, 9).Here, we report the discovery of the Landau-phonon polaritons (LPPs) in a 2D graphene-hexagonal boron nitride (hBN) hetero-structure. The LPPs result from the hybridization (19–21) of pho-non polaritons of the hBN encapsulating layers (8–10) with Dirac magnetoexcitons (6, 7) [or “Landau polaritons” (5)] of charge-neutral graphene (6, 7). LPPs belong under a broader umbrella of magneto-phonon resonance (MPR) effects, resulting from a near coincidence of the energy spacing between a pair of Landau levels (LLs) and the energy of an optical phonon. We comment on other MPR effects (23–28), such as magneto polarons (23–25), in the dis-cussion. We also note that, in addition to graphene, Landau polaritons in semiconductor 2D electron gas systems have been demonstrated extensively at terahertz frequencies, which exhibit intriguing ultra-strong light-matter coupling phenomena in a solid-state cavity quan-tum electrodynamics system (29). Using the state-of-the-art magneto scanning near-field optical microscopy (m-SNOM) (6, 30–32), we have imaged real-space interference patterns created by the LPPs in a graphene-hBN heterostructure. We demonstrate that the LPP propa-gation can be switched on and off using magnetic fields. We have been able to detect as many as six different LPP branches. Several of them originate from optically dark transitions, suggesting that the usual se-lection rules (33–39) no longer apply in the extended momentum-frequency space accessible with the m-SNOM. Our high-precision mapping of the LPP dispersion has also enabled us to quantify many-body effects that yield the effective Fermi velocity in gra-phene (33–37, 40, 41).Our experimental setup is depicted in Fig. 1A. The experiments involved focusing infrared radiation onto the tip of an atomic force microscope that acted as a scannable nanoscale antenna. Light scat-tered from the tip carried near-field information to a far-field detec-tor. Another, stationary nanoantenna in the form of a metallic bar deposited on graphene played the dual role of an electrical contact and a polariton launcher. Both the sample and the m-SNOM resided in an optical cryostat allowing the control of temperature and mag-netic field applied in the out-of-plane direction [see Materials and Methods and (6)]. We present and discuss the results of these mea-surements below after we have introduced the necessary theoretical background.RESULTSHigh-momentum magneto-optics of grapheneIn a transverse magnetic field, the density of states in graphene splits into LLs of energy En = sgn(n)√2�n��ℏvF ∕ lB� , where n = 0, ±1, ±2, 1Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794, USA. 2National Synchrotron Light Source II, Brookhaven National Laborato-ry, Upton, NY 11973, USA. 3Department of Physics, Columbia University, New York, NY 10027, USA. 4Department of Physics, University of California, San Diego, La Jolla, CA 92093, USA. 5Condensed Matter Physics and Materials Science Division, Brookhaven National Laboratory, Upton, NY 11973, USA. 6Research Center for Elec-tronic and Optical Materials, National Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan. 7Research Center for Materials Nanoarchitectonics, Na-tional Institute for Materials Science, 1-1 Namiki, Tsukuba 305-0044, Japan. 8De-partment of Physics, Yale University, New Haven, CT 06520, Fairfield University, Department of Physics, Fairfield, CT 06824, USA. 9Energy Sciences Institute, Yale University, West Haven, CT 06516, USA. 10Fairfield University, Department of Phys-ics, Fairfield, CT 06824, USA. 11Department of Quantum Matter Physics, University of Geneva, 1211 Geneva, Switzerland.*Corresponding author. Email: mengkun.​liu@​stonybrook.​edu (M.L.); db3056@​columbia.​edu (D.N.B.); mfogler@​ucsd.​edu (M.M.F.); lwehmeier@​bnl.​gov (L.W.)Copyright © 2024 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC). Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024mailto:mengkun.​liu@​stonybrook.​edumailto:db3056@​columbia.​edumailto:db3056@​columbia.​edumailto:mfogler@​ucsd.​edumailto:lwehmeier@​bnl.​govhttp://crossmark.crossref.org/dialog/?doi=10.1126%2Fsciadv.adp3487&domain=pdf&date_stamp=2024-09-13Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e2 of 16… is the LL index, vF is the Fermi velocity, e is the elementary charge, lB =√ℏ∕e ∣B ∣ is the magnetic length, and B is the magnetic field (Fig. 1B). This characteristic square-root n- and B-dependence is a manifestation of the Dirac-like energy-momentum dispersion of graphene quasiparticles. In a charge-neutral graphene, the optical transitions can occur between LLs with indices of opposite sign, −n → n′, at frequencies ω ∝√�n� +√∣n� ∣ . The oscillator strength of each transition is a function of the in-plane momentum k. Conven-tional far-field infrared experiments excite graphene at very small k, with a nonnegligible oscillator strength only at |n| − ∣n� ∣ = ± 1. This selection rule at k → 0 is evident in the nonlocal optical con-ductivity σ(ω,k) of graphene shown in Fig.  1C (blue curve in Fig. 1D). Such peaks in the optical conductivity Re σxx have been observed in many magneto-optical absorption experiments (33, 35, 36, 38, 39, 42).The theory also predicts that at finite momenta transitions be-tween any pair of LLs become possible (43). Among these addition-al “forbidden” transitions, the first ones to become noticeable as k increases are the –n → n transitions; see Fig. 1C. Faint signatures of such forbidden modes have been seen in previous far-field experi-ments (33). They were attributed to mildly relaxed momentum conservation due to disorder scattering. As discussed below, our ex-periments have revealed much stronger evidence of the forbidden inter-LL transitions (ILTs), presumably because the requisite large in-plane momenta were created by scattering of light with the tip. The forbidden transitions become comparable in strength to the nearby allowed ones at momenta of the order of the inverse mag-netic length, e.g., l−1B= 71 μm−1 at B = 3.35 T. The momentum range important in the m-SNOM is illustrated by the bell-shaped curve in Fig. 1C. For the estimated tip radius of rtip = 30 nm, it is centered at k = 33 μm−1 marked by the vertical dashed line (14). At such k, the forbidden transitions are only slightly weaker than the allowed ones; see the orange line in Fig. 1D. In addition, the allowed transitions at nonzero k are diminished with respect to the k = 0 case (the blue line in Fig. 1D) to fulfill the optical sum rule.Modeling of polariton dispersionEach ILT gives rise to a collective excitation, which has been previously referred to as a Landau polariton (5) (the term we use here), Dirac magnetoexciton, or magnetoplasmon (6, 7). If the Landau polari-tons are tuned in resonance with the hyperbolic phonon-polaritons in hBN by changing the applied magnetic field, the hybrid modes, which are the aforementioned LPPs, can form. We have carried out numerical simulations to model the LPP dispersion expected un-der our experimental conditions. As customary in near-field stud-ies, we deduce the dispersion of the collective modes from the frequency- and momentum-dependent p-polarized reflection coef-ficient of the sample, rp = rp(k, ω). Figure 2 (A to C) demonstrates the imaginary part of rP calculated for three representative values of the magnetic field. The multiple branches of phonon polaritons in the up-per Reststrahlen band of hBN (~1360 to 1610 cm−1) are evident in all three cases (8–10). Without the magnetic field (Fig. 2A), the charge-neutral graphene influences the response of the heterostructure only weakly via its “universal” optical conductivity σ = e2/4ℏ (33).At 3.35 T (Fig. 2B), the frequency of the −1 → 2 ILT is inside the hBN upper Reststrahlen band, which generates avoided crossings in the polariton dispersion. These features manifest as a coupling and hybridization of the −1 → 2 inter-LL Landau polariton with the Fig. 1. High-momentum magneto-optics of graphene. (A) Schematics of our sample and m-SNOM setup. Gold contacts enable transport measurements and gating of graphene and also serve as polariton launchers. (B) LL energy as a function of magnetic field B and LL index n = 0, ±1, …, ±4. Black (red) arrows mark −n → n ± 1 and −n → n inter-LL transitions (ILTs) for a photon energy of ℏω = 188 meV (1519 cm−1). (C) Analytically calculated real part of the graphene conductivity (43) at B = 3.35 T as a function of frequency ω and in-plane momentum k calculated using Fermi velocity vF = 1.19 × 106 m/s and damping γ = 24.3 cm−1. The relevant −n → n ± 1 (−n → n) ILTs are labeled in black (red). The black bell-shaped curve illustrates the momenta accessible via m-SNOM (14), which peak at around k = 1/rtip = 33 μm−1 as marked by the vertical dashed line. (D) The line cuts at momenta k = 0 and k = 33 μm−1 extracted from (C).Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e3 of 16hBN phonon polaritons, i.e., the formation of the LPPs. Modeling the system as two coupled harmonic oscillators (44, 45) (see Materials and Methods), we extract the mode splitting Ω = 43.3 cm−1 at the largest avoided crossing and mode linewidths of ΓLandau = 50.6 cm−1 and ΓhBN = 4.0 cm−1 for the uncoupled Landau polariton and hBN phonon polariton, respectively. Hence, the strong coupling criterion C =2Ω2Γ2Landau+Γ2hBN= 1.5 > 1.0 is fulfilled at this magnetic field.At 6.0 T (Fig. 2C), our calculations indicate no ILTs inside the Reststrahlen band, so the phonon-polariton dispersion is again largely unaffected by graphene. Notably, these calculations show that the polariton damping at 6.0 T should be lower than that at 0 T because the LL quantization makes graphene more optically trans-parent away from the discrete ILT frequencies (33).Nano-imaging of LPPWe now turn to our experimental nano-imaging results that reveal field-tunable features of LPPs in real space. Figure 2 (D to F) shows m-SNOM images acquired at a temperature of 154 K and mag-netic fields of B = 0.0, 3.35, and 6.0 T, matching Fig. 2 (A to C, re-spectively). The incident photon energy is 188 meV (wavenumber, ω = 1519 cm−1). Note that our scan area contained three different regions: 1) a gold electrode on the left, which served as a polariton launcher; 2) hBN-encapsulated graphene on the top right, showing propagating LPP polaritons; 3) hBN without graphene in the bot-tom right, showing phonon-polariton modes only. At 0.0 T (Fig. 2D) and 6.0 T (Fig.  2F), we observed polariton fringes parallel to the gold electrode in both regions 2) and 3). At 0.0 T, the fringes in the region containing graphene exhibited a higher damping. At 6.0 T, the impact of graphene was minimal. These findings are consistent with our simulations (Fig. 2, A and C) and also previous work (33). On the other hand, at 3.35 T (Fig. 2E), there is a notable contrast between the regions with and without graphene. The polariton propagation in hBN-graphene ceases such that all but the first fringe is suppressed. This gives clear evidence for the existence of the hy-bridization gap in the LPP dispersion, i.e., the strong mode cou-pling, predicted by our theoretical calculations (Fig. 2B).To study the magnetic-field dependence of the LPP dispersion in detail, we obtained a field-tip position map of the m-SNOM signal (Fig. 3A) by sweeping B from −6.0 to +6.0 T. The maps were ac-quired by performing repeated scans with the tip along lines per-pendicular to the gold electrode, as marked by the black arrow in Fig. 2F. At our selected photon energy of 188 meV (ω = 1519 cm−1) within the hBN Reststrahlen band, we observe the suppression of the fringes for certain distinct field values, e.g., for the discussed case of B  =  3.35 T. When approaching such fields from a higher (lower) absolute magnetic field side, the polariton wavelength de-creases (increases) along with an overall reduction in near-field sig-nal and a decrease of the propagation length. Figure 3B shows line profiles that have been extracted at B=0.0, ±3.3, and ±5.8 T, re-spectively. While we observe oscillatory polariton fringes at 0.0 and ±5.8 T, at the −1 → 2 ILT at B = ±3.3 T, the fringes are strongly damped, consistent with the predicted observations in Fig. 2. These features are observed for both directions of the magnetic field, B > 0 and B < 0. We note that the ILTs can also be suppressed by doping graphene off charge neutrality [via the Pauli blocking (33, 35, 39)], Fig. 2. Hybridization of hBN phonon polaritons with graphene Landau polaritons, resulting in LPPs. (A to C) Calculated LPP dispersion at magnetic fields of 0.0, 3.35, and 6.0 T, respectively. The false color represents Im rp(k, ω), the analytically calculated imaginary part of the reflection coefficient for p-polarized light. Graphene is as-sumed to be charge neutral with a constant LL broadening (33) γ = 24.3 cm−1 and Fermi velocity vF = 1.19 × 106 m/s, the latter being the value extracted from Fig. 4C. Inset in (B): An enlarged view of the region exhibiting strong coupling and an avoided crossing between the Landau and the hBN phonon polaritons; the arrow marks ω = 1519 cm−1 corresponding to the data in (D) to (F). (D to F) Nano-imaging data collected from the region marked by the red rectangle in Fig. 1A at T = 154 K and B = 0.0, 3.35, and 6.0 T, respectively. The near-field signal S3 (demodulated at the third harmonic of the tip frequency; refer to Materials and Methods) shows relative differences between regions with and without graphene that strongly depend on the magnetic field. The enhanced signal-to-noise ratio in (E) and (F), compared to (D), is due to a slightly longer integration time. We also note that the mechanical stability of our system is slightly better at higher magnetic fields. The double-headed arrow in (F) marks the location of the line scan analyzed further in Fig. 3. a.u., arbitrary unit.Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e4 of 16which provides additional opportunities for controlling LPPs (see Materials and Methods).Averaging the data in Fig. 3A over the tip position removes the spatial oscillations, which allows us to focus on the B dependence of the signal (Fig. 3C). We can compare the trace in Fig. 3C with the theoretical simulations in Fig. 3D. For simplicity, in these simula-tions, we did not include a separate polariton launcher. Instead, we modeled a more common m-SNOM setup where the tip is the only nanoantenna interacting with the sample, which is uniform and in-finite in size. It is fair to compare the predictions of this model with Fig. 3C since in both cases, the signal depends only on the magnetic field not the tip position. We see that Fig. 3 (C and D) are in good agreement. The dips in the theoretical curve come from distinct ILTs. Assuming that this is also the case in the data, we can label them accordingly. Furthermore, fitting the data to the theory allows us to determine the graphene Fermi velocity, which we discuss in detail below. Unassigned dips in the experimental spectrum are mi-nor and might be attributed to the experimental conditions, such as mechanical vibrations.We find that the dips corresponding to the −1 → 2 transitions are the strongest in our frequency window, testifying to a strong mode coupling regime (see also Materials and Methods). In addition, we observe clear signatures of several other transitions. They include allowed transitions −2 → 3 and −3 → 4, as well as transitions −1 → 1, −2 → 2, and −3 → 3 forbidden by the standard selection rules (33, 36, 38, 39). Any of these ILTs also induces a clear modification of the polariton wavelength as well as a reduction of the quality factor, which further supports our assignment and will be used for a quan-titative analysis later. In total, we can resolve six different ILTs in our data. Notably, the forbidden transitions (26, 33, 46) show up much stronger compared to what was previously seen in far-field infrared spectroscopy (33). As hypothesized above, this massive breakdown of the selection rules originates from the greater role of high-momentum field components k ∼ l−1B in our m-SNOM measure-ments (Fig. 1, C and D).Tunability of LPPsWe have fitted the polaritonic fringes in Fig. 3 to exponentially de-caying sine waves ∼eikx, where k = Re k + i Im k is the complex po-lariton momentum (Materials and Methods). From this fitting, we deduced the LPP wavelength λP = 2 π/ Re k (Fig. 4A) and the qual-ity factor Q = Re k/ Im k (Fig. 4B). For example, at B = 0 T, we found λP = 647 nm and Q = 12. As B increases, the crossing of each ILT results in a deep minimum of the Q factor as well as a characteristic change in λp. Within the studied magnetic field range, we have ob-served a modulation depth of λmax/λmin ∼ 2 (Fig. 4A) and Qmax/Qmin ∼ 10 (Fig. 4B). The latter is much larger than the values Qmax/Qmin ∼ 2 reported for gate-tuning of doped graphene-hBN polaritons (19). In particular, near the −1 → 2 transition, altering the magnetic field by only 10% changes the Q factor by a factor of five. The great tunability of our system can be explained by the large optical con-ductivity associated with the LL transitions at infrared frequen-cies and the resultant externally controllable gap in the polariton dispersion (Fig.  2B). Therefore, the magnetic field provides a Fig. 3. Magnetic-field dependence of the polariton dispersion. (A) Near-field signal S3 acquired via a repeated line scan while sweeping the magnetic field from −6.0 to 6.0 T at a rate of 0.4 mT/s; measurement was taken at ω = 1519 cm−1 and T = 154 K. The direction of this line scan was perpendicular to the gold contact to the left, which served as a polariton launcher; see Figs. 1A and 2F. (B) Line profiles extracted from (A) at magnetic fields B = 0.0, ±3.3, and ±5.8 T. (C) The near-field signal in (A) averaged over the distance. Minima of the averaged signal are assigned to the ILTs shown by the labels. The assignment is based on the calculation shown in (D). (D) Calculated near-field signal (Materials and Methods) as a function of magnetic field. Parameter values are chosen to be the same as in Figs. 1 (C and D) and 2 (A to C).Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e5 of 16feasible path toward on/off switching of polariton propagation in 2D systems.Many-body effectsLast, we discuss many-body effects manifested in deviations of the ILT frequencies from the √B - law valid for free Dirac fermions. An alternative description of these deviations is the renormalization of the effective Fermi velocity veffF defined by Eq. 1 below. From the minima in the Q factor, we can read off the magnetic fields B associ-ated with each ILT and obtain the corresponding field-dependent veffF . We find that veffF decreases with B for all transition types follow-ing a non-logarithmic B dependence (squares in Fig.  4C). These values agree unexpectedly well with previous far-field infrared (33) and Raman (34) spectroscopy results. We extract veffF associated with both allowed and forbidden ILTs with the same measurement. In this regard, the m-SNOM provides a unified approach for LL spectroscopy that lifts many previous limitations.Our theoretical calculations (Materials and Methods) of the ef-fective Fermi velocity (diamonds in Fig. 4C) show a good agreement with the data using one adjustable parameter, the value of veffF at one specific ILT (here, −2 → 3 gives the best agreement). These calcula-tions also corroborate our experimental observation that the effec-tive Fermi velocity of the −n → n ILTs is consistently below the trend followed by the −n → n ± 1 ones.Our explanation for the above observation of veffF is as follows: De-spite its common usage, the term “effective” Fermi velocity is some-what misleading in the present context. A more accurate statement is that the interaction corrections to the observed ILT energy ℏω, resulting in veffF , include contributions from both the Fermi velocity renormalization (a polaronic effect) and excitonic effects. Namely, ℏω is given by the LL energy difference, |En| + ∣En�∣, minus the magneto-exciton binding energy Δnn�The LLs En in this expression obey the quantization rule |En| = E(qn) where E = ℏvrenFq is the renormalized quasiparticle dispersion and qn = l−1B√2�n� is the quantized momentum of a Dirac fermion residing at the nth LL (inset Fig. 4C). The effective Fermi velocity veffF is (approximately) equal to the renormalized vrenF only if the magne-toexciton binding energy Δnn� is neglected. In that case, a logarith-mic dependence of veffF on E (at fixed n and n′) follows from the perturbation theory formula vF(E) ≈ vF(Λ)[1+ 14αln|Λ∕E|+…] where α = e2/(κℏvF) ≪ 1 is the Coulomb coupling constant, Λ is the high-energy cutoff, and κ is the effective dielectric constant of the graphene environment (36). This formula has been derived for ℏω ≡ ℏveffFlB�√2�n� +√2�n���=���En�� + ��En� ���− Δnn� (1)Fig. 4. Magnetic-field dependence of LPP properties and Fermi velocity renormalization. (A) Polariton wavelength λP and (B) polariton quality factor Q = Re k/ Im k as a function of the magnetic field B. Solid lines show experimental values extracted from Fig. 2; shaded regions show the SD of the measurement. (C) Effective Fermi velocity veffF as a function of the logarithmic magnetic field ln(B) derived for different ILTs (see section Many-body effects): Squares show experimental values derived from (B). Diamonds represent calculated values of veffF (see Materials and Methods) (37). We observe a nonlogarithmic trend. Inset: The red (black) points show veffF for the −1 → 1 (−1 → 2) ILT measured via Raman spectroscopy (34) [far-field infrared spectroscopy (33)]. The tapering shape of the Dirac cone illustrates the Fermi velocity renormaliza-tion (33, 34, 40), resulting in a logarithmic B dependence of the far-field data (33, 34). (D) Squares and diamonds show the exciton binding energy Δnn� of the Landau po-laritons derived from the experiment and theory, respectively. The exciton binding energy is larger for the ILTs with n′ = −n compared to those with n′ = −(n ± 1) and generally increases with magnetic field. Inset: The dependence of the exciton binding energy on the magnetic field and type of the ILT can be explained within a semiclas-sical model where quantized electronic orbitals of the LLs are shaped as narrow rings of radius rj = lB√2�j� , j = n, or n′. The magnetoexciton binding energy Δnn� (see text) is given by the Coulomb attraction energy of these rings. For a fixed n, this binding energy is the largest when the ring radii are equal, at n′ = −n.Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e6 of 16graphene in zero magnetic field; however, it remains approximately correct at nonzero B (Materials and Methods), meaning that the re-normalized Fermi velocity is first and foremost a function of energy, vrenF= vrenF(E) . Since En and En′ are B dependent (33, 34), vrenF usually acquires a logarithmic B dependence for a given ILT, as found in previous far-field spectroscopy studies (inset in Fig. 4C) (31, 32, 34, 35, 38). On the other hand, here, we have studied ILTs at a fixed laser frequency so that the transition energy ℏω ≈ |En| + ∣En�∣ re-mained the same, being split roughly equally between |En| and ∣En� ∣. Therefore, in our experiments, the renormalized Fermi velocity vrenF≈ vrenF(ℏω∕2) should have little B-field dependence and the ob-served variation of veffF (Fig.  4C) should mostly come from the change of magnetoexciton binding energy Δnn�, which follows non-logarithmic trend with changing magnetic field. The data analysis in some of the previous works was oversimplified by emphasizing the logarithmic trend and not considering the additional field-dependent contribution due to the magnetoexciton binding energy.Our theoretical calculation of the two competing terms, |En| + ∣En�∣ and Δnn�, in Eq. 1 confirms that at fixed ℏω, the former gives a nearly constant contribution to veffF for all measured ILTs so that veffF variation comes from the latter, with characteristic dips occurring at n′ = − n points (Materials and Methods). This allows us to extract the binding energy Δnn� (Fig. 4D) from veffF . The absolute value of Δnn� generally increases with magnetic field and is larger for the ILTs with n′ = − n compared to those with n′ = − (n ± 1). A simple way to think about the magnetoexciton binding energy Δnn� is to imagine that it is equal to the Coulomb attraction energy of two LL orbitals shaped as concentric rings, one with charge +e, the other with charge −e (inset Fig. 4D). The ring radii are given by the formula rj =∣Ej ∕evrenFB ∣ , where j = n or n′, which is the semiclassical cyclo-tron radius of a Dirac particle with energy Ej. (Note additional rela-tions rj = lB√2||j|| = l2Bqj .) For a fixed n, this attraction energy is the largest when the ring radii are equal, i.e., at n′ = −n, yielding the lowest veffF at such ILTs.DISCUSSIONOur study has shown that the physics of LPPs is very rich, and it involves simultaneously three types of effects: polaritonic, excitonic, and polaronic. These effects have distinct characteristics: 1) The po-laritonic effects change collective mode properties in the hetero-structure. Forbidden optical transitions are now accessible in the momentum space offered by m-SNOM. The mode coupling be-tween Landau polaritons (magnetoexcitons) in graphene and pho-non polaritons in hBN generates a tunable avoided crossing, which could potentially be further tailored by using other ILTs (e.g., 0 → 1 ILT) or multilayer engineering (e.g., adding additional layers of gra-phene). 2) The excitonic effects are manifestations of the electron-electron interactions. They lead to a finite binding energy, which also modifies the LPP dispersion. This binding energy can be fur-ther tuned via dielectric screening engineering. 3) The polaron effect is another term for the renormalization of the quasiparticle dispersion. Although above we emphasized the role of electron-electron interactions as the reason for the renormalization of the Fermi velocity vrenF , this interaction is screened by hBN. Hence, the interaction of electrons in graphene with phonons in hBN is includ-ed implicitly. In our case, vrenF does not change much with magnetic field since we keep the incident photon energy the same throughout the experiments. In addition, in our calculation of the renormal-ized Fermi velocity, we approximated the hBN dielectric function by its dc (ω = 0) value. Goals for future work can be: i) incorporat-ing more sophisticated theoretical approaches into our model to properly treat electron-phonon coupling and ii) experimentally studying in more detail the momentum dependence (e.g., via differ-ent m-SNOM tip radii) and frequency dependence of LPPs.As mentioned in the introduction to this paper, LPPs are specific examples of MPR effects. Other known MPR effects include magneto-polarons (23–25), dc magneto-transport oscillations (27, 28, 47), and mode splitting in magneto-Raman spectroscopy (26). Most of them have been studied in bulk crystals or a single material system. It would be interesting to investigate whether these phenomena are af-fected by finite-momenta LPPs in a 2D heterostructure. Last, it would be desirable to explore a variety of other nano-magneto-optics phe-nomena using m-SNOM, including chiral edge magnetoplasmons (48–50), cavity magneto optics (51), magnon polaritons and magnon-phonon polaritons (5, 52), the polaritonic Hofstadter butterfly (53), magnetoexcitons of fractional quantum Hall states (54), and collec-tive modes of stripe phases in partially filled LLs (55).MATERIALS AND METHODSExperimental setup: Magneto infrared nanoscopyDetails of our experimental setup have been described in references (6, 56). Infrared nano-imaging in magnetic fields up to 7 T is dem-onstrated using a home-built scattering-type scanning near-field optical microscope (57) that is placed within a closed-cycle cryostat (OptiCool, Quantum Design). The accessible sample temperature with our scanning probe system is ~10 to 350 K. For infrared near-field imaging, a tunable QCL laser (Hedgehog mid-IR laser by DRS Daylight Solutions) is focused onto an atomic force microscopy tip via a parabolic mirror. As the resulting light spot is diffraction lim-ited, it is much larger than the atomic force microscopy tip. Hence, it also illuminates the gold electrode that here serves as a launcher for the polariton waves (refer to the sketch of the experimental setup in Fig. 1A). The light scattered from the tip contains information about both the sample material and its polaritonic excitation. We detect the scattered light with a mercury-cadmium-telluride detec-tor via a self-homodyne detection scheme (11, 57, 58), which yields a stable scattering signal and a good signal contrast. To separate the near-field signal from far-field contributions to the detected optical signal, we use lock-in demodulation at higher harmonics nΩtip of the tip tapping frequency Ωtip, with n = 3 for the data shown in the manuscript. We use an Akiyama-type scanning probe with a resonance frequency of about Ωtip ≈ 65 kHz (56), and scanning stages and positioners from Attocube (ANSxyz100/LT/UHV and ANPxyz101/LT/UHV, respectively).Fabrication of graphene-hBN heterostructureThe sample studied in this work is hBN-encapsulated monolayer graphene. We fabricated the hBN/graphene/hBN stack using the standard van der Waals dry transfer technique (59). Graphene and hBN flakes were exfoliated onto Si/SiO2 (285 nm) substrates. The flakes were examined using an atomic force microscope and those with pristine surfaces and suitable thicknesses were selected. A polycarbonate (PC)/ polydimethylsiloxane (PDMS) was prepared and used to assemble the stack. We used the stamp to pick up the top hBN at ~120°C; then, monolayer graphene was picked up by Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e7 of 16the top hBN at ~100°C; and then, the bottom hBN was picked up by the graphene layer at ~120°C. The stack on the stamp was trans-ferred onto a clean Si/SiO2 (285 nm) substrate at 180°C, where the PC melted and detached from the PDMS. The finished stack was annealed in forming gas (95% N2 and 5% H2) at 450°C for 4 hours to minimize the polymer residue on the surface. Electrical contacts to the graphene channel were defined via standard electron beam lithography. The top hBN and graphene at the contact area were etched in CHF3/O2 (10:1) plasma (40 mTorr, 60 W), exposing the side of the graphene layer. Cr (5 nm)/Au (40 nm) was thermally evaporated to form side contacts to graphene (60).Pauli-blocking and gate-switchable dispersionA transition n → n′ between two LLs can be suppressed (Pauli-blocked) by emptying the initial LL or filling the final LL via doping graphene off charge neutrality (33, 35, 39). Hence, Pauli-blocking is another clear signature of the ILTs. Here, we show that Pauli-blocking can be used to tune the LPP dispersion, requiring lower carrier doping than nonmagnetic tuning. This is in stark contrast to the tuning mechanism of plasmons in doped graphene that de-pends on modifying the plasma frequency via the charge carrier density (13, 15). Studying the same line as in Fig. 3A, Fig. 5A shows a measurement at constant magnetic field strength B = 3.3 T, with filling factor ν = 2πℏN/eB varied from 0 to 19 (N is the charge car-rier concentration). Figure 5B shows corresponding line profiles at distinct filling factors ν = 0,5,10,15. At low ν, we find a strongly damped polariton, matching well to our magnetic-field dependent measurements of charge-neutral graphene above. When increasing the filling factor, due to the Pauli-blocking, we regain a propagating polariton mode with clear oscillatory behavior.A quantitative analysis of the line profiles yields the filling-factor dependence of the polariton wavelength λP (Fig.  5C) and quality factor Q (Fig. 5D). At a constant magnetic field of 3.3 T, for charge-neutral graphene with ν = 0, λP is larger than its value without a magnetic field due to the coupling with the −1 → 2 inter-LL reso-nance. However, the polaritons are strongly damped, preventing a Fig. 5. Gate dependence of the −1➔➔ 2 graphene LPP at B = 3.3 T. (A) Near-field signal acquired via a repeated scan of the same line as in Fig. 3A while sweeping the gate voltage and, thus, the filling factor from ν = 0 to 19; measurements were done at ω = 1519 cm−1, and T = 154 K. (B) Line profiles extracted from (A) at four different filling factors of ν = 0, 5, 10, and 15, respectively. (C) Polariton wavelength λP and (D) polariton Q factor as a function of ν and charge carrier concentration N at 3.3 T; the shaded regions show the SD of the measurement.Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e8 of 16meaningful estimate of λP at ν = 0. When increasing the filling factor to ν ≈ 10, λP decreases to a local minimum and approximately re-gains the wavelength observed at 0 T (compare to Fig. 3A). At the same time, the Q factor increases by more than a factor of 5. This is readily attributed to the expected Pauli-blocking of the −1 → 2 transition: One expects the second LL to be filled at a filling factor of ∣ν∣ = 4 ∣n�∣ + 2 = 10 with n′ referring to the index of the highest LL contributing to the transition (33), i.e., here, n′ = 2. When increasing ν beyond 10, the Q factor remains almost constant while λP slowly increases again. This increase of λP is not directly related to the Pauli-blocking of the inter-LL excitation and can rather be assigned to the hybridization between hBN phonon polaritons and plasmon polariton modes in doped graphene (19). Overall, combining gate tuning with a constant magnetic field opens unique possibilities for polariton tunability in the low-voltage regime, with a charge carrier concentration difference of less than 1011 cm−2, resulting in a polari-ton wavelength change of several percent.Additional insight into the coupling between graphene Landau polaritons and hBN phonon polaritons can be attained by studying the frequency dependence of the LPP dispersion along with its gate dependence. Like Fig. 5A, Fig. 6 shows near-field images that were acquired by repeatedly scanning the same line perpendicular to the gold electrode while sweeping the back gate voltage. The images have been acquired at six different photon frequencies; in Fig.  6 from left to right, these are: ω = 1498,1509,1519,1529,1537, and 1547 cm−1. From top to bottom, the gate voltage was changed from 0 to −50 V for each scan (0 to −40 V for the scan at 1529 cm−1).From the symmetry of the near-field images, one can see that the charge neutrality point (CNP) of graphene for these scans is found approximately in the middle of the scans, i.e., around −25 V. All the measurements shown so far (Figs. 2 to 5) have been taken in a single measurement run at 154 K. For the measurements shown in Fig. 6, the sample had been heated up to room temperature (to exchange the m-SNOM tip) and then cooled to 76 K. While both the thermal cycling (33) as well as the lower temperature may affect the mea-surement, the results in Fig. 6 appear consistent with the results dis-cussed so far. As also observed in Fig. 5, Fig. 6 illustrates how the ILT gets Pauli-blocked when doping graphene away from charge neutrality. Well matching to the calculated dispersion in Fig. 2 (A to C) and literature results on the uncoupled hBN phonon polaritons(8), the LPP fringes show a clear positive dispersion, i.e., increasingpolariton momentum (decreasing wavelength) with increasing fre-quency. Notably, Fig. 6 also directly illustrates how the gap in theLPP dispersion (Fig. 2B) opens and closes with frequency: While,around charge neutrality, the LPP propagation is gapped out at thecentral frequencies in Fig. 5 (e.g., at 1519 and 1529 cm−1), one can see the width of the gap (along the voltage axis) decrease when mov-ing to either lower or higher photon frequencies. For example, both at 1498 and 1547 cm−1, the fringe suppression due to the coupling with the graphene ILT is less efficient, and some fringes can even be observed at charge neutrality. Again, this behavior agrees very well with the calculated dispersion shown in Fig. 2B.Numerical calculation of graphene magneto- optical conductivity, polariton dispersion, and near- field signal: (Figs. 1, C and D, 2D, and 3D)We calculate the longitudinal conductivity of graphene (plotted in Fig. 1, C and D of the main text) using the approximate Kubo formulaIn this formula, the effect of disorder and electron interactions on LL broadening is expressed using a single phenomenological param-eter γ. Equation 2 neglects contributions from intra-LL transitions, which are model dependent and difficult to compute. [For an example of such calculation, see (61).] These intra-LL contributions should be small at frequencies ω ≫ γ of interest and should be absent complete-ly under the integer quantum Hall effect conditions, where each of the Fermi factors fn is equal to either zero or unity (no partial LL filling occurs). In principle, the matrix elements Fnm(q) in Eq. 2 can be sepa-rated into two parts, Fnm(q) = F0nm(q) + δFnm , where the first one rep-resents the response of an ideal system without electron scattering and the second one is a correction due to interactions and disorder. If γ and δFnm originate from the same physical mechanism, such as the scattering of electrons by static disorder, then the matrix element δFnm are expected to scale as γ2/(Em − En)2. In the aforementioned paper, for short-range disorder and in the absence of electron-electron interactions, the following result has been derived (61)In our experiments, we typically have γ2 ∕ω2c∼ 10−5 … 10−4 so that δFnm should be negligible. We do not consider such corrections below; however, the question of whether electron interaction can produce larger δFnm for the same small γ is conceptually interesting and may warrant future theoretical study.σxx(q,ω)= ie2ω∞∑m=−∞∞∑n=m+1(Em − En)(fm − fn)(Em−En)2− ℏ2(ω+ iγ)2Fnm(q)(2)δFnm=γ24π(1+δ0,m)(1+δ0,n)(Em−En)2[2|m|+2|n|+δn,0(2|n|−2|m|+1)2+(m↔n)] (3)Fig. 6. Frequency dependence of the gate dependence at B = 3.3 T. Near-field signal acquired via a repeated scan of approximately the same line as in Fig. 3A while sweeping the back gate voltage from 0 to −50 V. The charge neutrality point (CNP) of graphene for these scans is found approximately in the middle of the scans, i.e., around −25 V. To investigate the dependence on the photon energy, the scans were performed at different frequencies, from left to right: ω = 1498,1509,1519,1529,1537,1547 cm−1. The measurements were done at temperature T = 76 K. Scale bar, 1 μm.Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e9 of 16The formulas for Fnm(q) in the ideal system, which are well known, can be written as (43)where Lαn(x) is the associated Laguerre polynomial. At small mo-menta q ≪ l−1B , function Fnm(q) scales as q2||m|−|n||−2 so that the con-ductivity σxx(q, ω) is dominated by the allowed transitions with |m|−|n| = ± 1, as discussed in the main text. Precisely at q = 0 (“far-field limit”), the conductivity is given byNote that because of the factors 1 + δ0,n in Eq. 4, the oscillator strength of the principal cyclotron transitions 0 → ± 1 is effectively doubled compared to other ILTs. In Eq. 5, this doubling is affected implicitly by summing over the band index s.In the case of zero temperature, which we consider from now on, the numerators in the series of Eq. 5 can have only three possible values: Nms = 0, −1, or −2. For example, at CNP, they are equal to Nm,−1 = − 2 if s = − 1 (inter-band transitions) and Nm,1 = 0 if s = 1 (intra-band transitions) for all m except m = 0, for which N0,−1 = N0,1 = − 2. The series for finite doping differs from that for the CNP only in a finite number of terms for which Nms deviates from this rule. Hence, the far-field limit of the conductivity of a doped system can be easily computed from that for CNP by adding suitable cor-rections due to these terms.Numerical evaluation of the series in Eq. 5 at CNP poses a prob-lem, however, because of its slow convergence. To expedite the con-vergence, we transformed this series into an integral using standard techniques of complex analysis. The final expression for σCNPxx is somewhat lengthy, so we split it into three partsThe first part C1 is a partial sum of the original series (cf. Eq. 5),The upper limit of summation M is largely arbitrary, and the minimal choice (where M = 0 or 1) is given below. The second part C2 is the integralWe ensure that this integral is free of singularities by taking M = 0 if Re mω > 1 where mω =14(ω−1ω)2and M = 1 otherwise. Parameter mω has the property that at m = mω, the summand in Eq. 6b has a pole in the complex plane of m. The final part of our formula comes from the residue of that polewhere Θ(x) is a unit step function. Numerical evaluation Eq. 6a to d is very fast, yielding the numerically exact results for conductivity σCNPxx(0,ω) at CNP and thus σxx(0, ω) at any desired doping.Next, we discuss our procedure for computing σxx(q, ω) at non-zero q. Here, both allowed and forbidden ILTs must be included. Once again, the problem is difficult because the terms in Eq. 2 de-crease slowly with m and n. For simplicity, we focus on the CNP where only the inter-band transitions are possible. Applying the Wentzel-Kramers-Brillouin approximation to solve the differential equation for Ln2−n1n1(x) , we derived the following asymptotic formula, valid for m ≥ n ≫ 1Substituting this formula into Eqs. 4 and 2, we see that the high-order terms in Eq. 2 scale as O(x−1/2 n−2) if n≫ ω and that the dom-inant contribution comes from m’s that satisfy the double inequality √n ≤√m ≤√n +√x . Let us define the series residual ΔσCNPxx as the difference between σCNPxx  and the series of Eq.  2 truncated at some Nc so thatFrom the above reasoning, we expect that ΔσCNPxx scales as N−1∕2c  if Nc ≫ ω. Hence, for strong enough magnetic fields where ω = (ω + iγ)∕ωc  is at most of the order of unity, we can achieve, say, 5% accuracy by simply summing over N2c∼ 0.05−4 ∼ 105 terms. As the magnetic field decreases, ω increases, and this naïve strategy be-comes unfeasible because of long computation times and, more im-portantly, the accumulation of roundoff errors in evaluating the Laguerre functions. As an alternative, we briefly experimented with various series acceleration methods (e.g., Levin’s method) but found them to be numerically unstable due to their own set of roundoff errors. Eventually, we settled on estimating the residual by using the asymptotic formula (Eq. 7) and replacing the summation with inte-gration. After some manipulations, our result, in the absence of damping, can be presented as followsFnm(q)=2πx|||sgn(mn) An<,n>(x)+An<−1,n>−1(x)|||2,n<=min(|m|, |n|), n>=max(|m|, |n|), x≡ q2l2B2,An1,n2(x)=(1+δ0,n121+δ0,n22n1!n2!xn2−n1)1∕2Ln2−n1n1(x)e−x∕2(4)(5)σCNPxx(0,ω) =e22πℏ(C1 + C2 + C3)(6a)C1 = −2iωM�m=01�√m+1+√m�2− ω21√m + 1 +√m(6b)(6c)C3 = iω4 − 12ω4πtanπmωΘ(|ω| − 1)Θ(Remω −M −12)(6d)(7)σCNPxx�q,ω�=e2ℏiωNc�n=0Nc�m=n√m+√n�√m+√n�2−ω2Fnm�q�+ΔσCNPxx�q,ω�(8)ΔσCNPxx�q,ω�=−e2π2ℏq2iω�q2v2F−ω2q∫0dk−�q2−k2−arccosC�k−�,C�k−�=k2+v2F−2q2v2F+ω2k2+v2F−ω2,k+ =max�q, 2qc−k−�, qc =√2NclB(9)Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e10 of 16In weak magnetic fields where x ≡ q2l2B∕2 > 4Nc , Eq. 9 can be evaluated analytically. Here, we have k+ = q, C(k−) = − 1 so thatwhich is a well-known expression for the conductivity of gra-phene at CNP in the absence of the magnetic field. Apparently, σCNPxx(q,ω)≃ ΔσCNPxx(q,ω) in this regime. At stronger fields where <4Nc, Eq. 9 can be easily evaluated by a numerical quadrature. Last, to include a small LL broadening γ ≪ ω, we replaced ω by ω + iγ in these equations. Including ΔσCNPxx in the calculation significantly im-proves the accuracy of the procedure so that reasonably accurate results can be obtained for modest Nc. We found that Nc = 70 is suf-ficient for the entire parameter space of momenta, frequencies, and magnetic fields we explored.The parameters we used for calculating σCNPxx(q,ω) according to Eqs. 4, and 8 to 10 were as follows. On the basis of our experimental results for the −1 → 2 transition (Fig. 4C), we choose a Fermi veloc-ity vF = 1.19 × 106 m/s. Our choice of the LL broadening γ is based on (33), which suggests an energy-dependent LL broadening of the form γ = γ0 + βE with γ0 = 0.75 meV and β = 0.012. For our photon energy E = ℏω = 188 meV, it gives γ = 3.01 meV (24.27 cm−1), which is the constant value that we used for our calculations.Having calculated graphene conductivity, next, we compute the reflection coefficient rp for the whole heterostructure using the trans-fer matrix formalism [see, e.g., (8, 13, 14)]. We assume that the sample has the following structure: hBN (16 nm), graphene, hBN (55 nm), SiO2 (285  nm), and Si (semi-infinite); see Fig.  1D. Representative plots for Im rp as a function of momentum, frequency, and magnetic field are shown in Fig. 2 (A to C) of the main text.The near-field signal in Fig. 3D is calculated as described in (62), assuming the tip radius of rtip = 30 nm, length of the long semi-axis of the spheroid of 160 nm, and tip tapping amplitude of A = 40 nm.Momentum-range accessible via m-SNOMAs m-SNOM is based on scattering-type SNOM (s-SNOM), the ac-cessible momentum range of both techniques can be described in the same manner, which originates in the simple point-dipole model (14, 57, 58, 63) of the tip-sample interaction. The bell-shaped curve in Fig.  1C represents the time-averaged function k2e−2kz(t). Here, the time-dependent effective dipole position z(t) is (14): z(t) = b + A(1 − cos Ωtipt), b is the shortest distance between the model dipole and the sample, A is the tip tapping amplitude, and Ωtip is the tip tapping fre-quency. The analytical formula for this bell-shaped curve is therefore 〈k2e−2kz〉t = k2e−2k(b+A)I0(kA), where I0(x) is the modified Bessel func-tion of the first kind. We choose b = frtip with f = 0.79 based on s-SNOM modeling literature (64, 65), and we assume rtip = 30 nm and A = 40 nm.Extracting the coupling strength via the coupled harmonic oscillator modelCoupled harmonic oscillators model for hBN phonon polaritons and graphene Landau polaritonsWe model the interaction between hBN hyperbolic phonon polari-tons and graphene Landau polaritons as coupling of two harmonic oscillators (44, 45, 66, 67). The oscillator resonant frequencies and damping parameters are ωi, and Γi, respectively, where i = {hBN, Landau}. The complex eigenfrequencies of this system are given bywhere the sign ± indicates the lower (−) and the upper (+) branch solu-tion (42). Here, ω =(ωhBN + ωLandau)∕2 , Γ =(ΓhBN + ΓLandau)∕2 , G is the coupling strength, and H = ωhBN − ωLandau −i2(ΓhBN − ΓLandau) . The real parts of the eigenfrequencies areBy calculating ω± and the parameters ωi and Γi of the uncoupled oscillators, we can determine the coupling strength G from Eq. 12. Because of the dispersive nature of polaritons, each parameter de-pends on the momentum k. The following sections describe the pro-cedure for extracting ωi(k), Γi(k), and ω±(k) from the simulated dispersions.Extraction of the oscillator parameters for uncoupled modesTo calculate the parameters of the oscillators, we simulate the imag-inary part of the reflection coefficient Im rp that yields the disper-sion for each uncoupled mode. For hBN phonon polaritons, we considered the Im rp map at B = 6 T shown in Fig. 2C in the main text (and Fig. 7A). The reason we preferred to use B = 6 T instead of B = 0 T is that at zero field, the charge-neutral graphene influences the polaritonic response of the heterostructure via its universal opti-cal conductivity σ = e2/4ℏ (33). We isolate the fundamental hyper-bolic branch by applying a window mask to the dispersion map. The mask values are 1 in the interest region and 0 elsewhere (red area in Fig. 7A). Figure 7B shows line profiles from the map extracted for different momenta. The oscillator parameters were obtained by fit-ting these profiles to the functionThe extracted dispersion ωhBN(q) (white line in Fig. 7A) matches well with the dispersion in the Im rp map. Figure 7C shows the ex-tracted values for fhBN and ΓhBN.To calculate the parameters for the Landau polaritons, we simu-lated Im rp for graphene on top of SiO2 (thickness, 285 nm). Figure 7D shows the map at B = 3.35 T. Using the same procedure described previously for hBN polaritons, we extracted the uncoupled oscillator parameters for the LLPs. Figure 7E shows the extracted ωLandau(B, k). Figure 7D (black line) shows ωLandau(k), and Fig. 7F shows ΓLandau and fLandau for B = 3.35 T.Dispersion of branches and coupling strength extractionWe extracted ω±(k, B) from Im rp by locating its maxima for every fixed k. The hBN phonon polariton dispersion was used as a refer-ence to separate the ω− and ω+ branches. The result of this extrac-tion is illustrated in Fig. 8 (A and D), for B = 3.15 T and B = 3.4 T, respectively. Figure  8 (B, C, E, and F) shows the spectral gap 2Δω±(k) = 2(ω± − ω) as a function of the complex parameter H(k). The coupling strength was obtained by fitting the spectral gap to Eq. 12.Coupling strength and criterium for strong couplingFigure  9A shows the calculated coupling strength G for different magnetic fields. We can observe that the coupling strength increases until it reaches a maximum at B = 3.3 T. The frequencies of hBN and ΔσCNPxx�q,ω�=e24ℏiω√q2v2 − ω2(10)ω± + iΓ±∕2 = ω − i Γ∕2 ±√G2 +14H2 (11)ω± = ω ± Re√G2 +14H2 (12)Im rp ∝ ωf 2iIm(1− ω2 − iωΓi + ω2i)(13)Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e11 of 16Landau modes in which dispersion crossing occurs are shown in Fig. 9B. By calculating the branches splitting Ωtip = ω+ − ω− where crossing occurs (Fig.  9C), we calculated the criterion C = 2Ω2 ∕(Γ2hBN+ Γ2Landau) (Fig. 9D) (44). Strong coupling is achieved for 3.15 T < B < 3.4 T, where C > 1.To verify the results, we show in Fig. 9 (E to H) the calculated frequency branches from Eq. 12, in which we used the parameters from the uncoupled hBN and Landau modes and the extracted cou-pling strength. Notice that we have calculated coupling values only in a narrow range of magnetic fields. This is because the coupled oscillator model is only well-defined at frequencies and momenta in which the dispersions of both modes overlap. Nevertheless, the hBN polariton dispersion is not fully recovered for B > 3.4 T (see Fig. 9H at B = 3.6 T), indicating that these polariton modes can still interact with each other.Weak polariton fringes parallel to graphene boundary in fig. 2 (D to F)Figure 2 (D to F) shows clear polariton fringes parallel to the gold electrode. The analysis of these fringes and their dispersion with magnetic field is a central part of this study. We note that weaker polariton fringes can be observed parallel to the graphene boundary marked by the dashed white line in Fig. 2 (D to F). This mode reflec-tion is attributed to the polariton momentum mismatch between hBN regions with and without graphene. Therefore, the interference fringes parallel to the graphene boundary are most pronounced in pure hBN at 3.35 T where the momentum mismatch between the regions 2) and 3) is the largest. In our analysis, we do not focus on these weak polariton fringes.Many-body effects: Fermi velocity and magnetoexciton binding energyFigure 2 (D to F): The effective Fermi velocity was defined by Eq. 1, copied below for convenience (assuming n, n′ ≥ 0)In the main text, we explained that the many-body correction ΔE−n→n� to the energy E−n→n� of a −n → n′ ILT (at zero momentum q = 0) has two parts: i) the correction due to the electron self-energies ΔEm of m = −n, n′ LLs and ii) the magnetoexciton binding energy Δ−n,n�We referred to the former as the Fermi velocity renormaliza-tion effect and to the latter as the excitonic effect. In this ap-proach, the renormalized Fermi velocity of −n → n′ ILT is given bywhere vF is the bare Fermi velocity, so that veffF,−n→n� and vF,−n→n� are related byveffF,−n→n�≡ lBℏE−n→n�√2n +√2n� (14)ΔE−n→n� =(ΔEn� − ΔE−n)− Δ−n,n� (15)vF,−n→n�≡ vF+lBℏΔEn�− ΔE−n√2n +√2n�(16)veffF,−n→n�= vF,−n→n�−lBℏΔ−n,n�√2n +√2n�(17)Fig. 7. Extraction of oscillator parameters from polariton dispersions. (A) Im rp at B = 6 T from the main text (Fig. 2C). The red-shaded region indicates the area ex-cluded by the window function. The white solid line is the phonon polariton dispersion extracted from fitting. (B) Line profiles from (A) extracted for different momenta values (black lines). The fitted curve is represented as dashed lines. (C) hBN phonon polaritons damping and oscillator strength extracted from fitting. (D) Im rp simulated for graphene on SiO2 (thickness, 285 nm) at B = 3.35 T. The black solid line represents the dispersion extracted from fitting. (E) Extracted Landau polariton resonance as a function of magnetic field and momenta. The white dashed lines represent isofrequencies. (F) Landau polariton damping and oscillator strength from fitting.Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e12 of 16Following previous work (37), we compute the quantities that enter these equations as followsHere, fm is again the Fermi occupation factor of mth LL [fm = Θ(−m) if m ≠ 0 and f0 = 1/2], gmn(q) is the matrix element defined in Eq. 4, and v(q)=2πe2ϵ(q)q is the screened Coulomb interaction. We ΔEn� −ΔE−n=− ∫d2q(2π)2v(q)I−n,n�(q),I−n,n�(q)=∑mfm(|||gm,n�(q)|||2−|||gm,−n(q)|||2) (18)Δ−n,n� = ∫d2q(2π)2v(q)g−n,−n(q)gn�n�(q)(19)Fig. 9. Coupling strength analysis and model validation across varying magnetic fields. (A) Coupling strength as a function of the magnetic field, extracted via fitting. (B) Ex-tracted spectral position where the hBN phonon polariton and Landau polariton dispersion cross. (C) Branch separation Ω at the crossing point. (D) Criterion C for strong coupling. The blue-shaded areas represent the uncertainty of the presented variables, assuming straight connections between the data points as a guide to the eye. (E to H) Calculated branch dispersion based on the uncoupled modes and the extracted coupling strength for (E) B = 3.2 T, (F) B = 3.3 T, and (G) B = 3.4 T. For (H) B = 3.6 T, we only show the hBN phonon polari-ton and the Landau polariton dispersions. The Im(rp) is presented on the background to show the expected dispersion from the simulation.Fig. 8. Determination of coupling strength from hybrid modes. (A and D) Extracted dispersion of lower (green) and upper (yellow) branches for (A) B = 3.15 T and (D) B = 3.4 T. The simulated Im rp and the dispersion of the uncoupled modes are shown for reference. Dependence of 2Δω±(q) on (B and E) the real part and (C and F) the imaginary part of H(k) for B = 3.15 T and B = 3.4 T.Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e13 of 16neglect dynamical screening effects; however, we include static screening through the momentum-dependent dielectric functionFunction κ(q) has the meaning of the effective dielectric constant of graphene environment. Parameter ϵhBN ≡ (ϵ⊥hBNϵ∥hBN)1∕2= 4.9 is the dc (i.e., ω = 0) dielectric constant of hBN defined as the geomet-ric average of its in-plane and out-of-plane dielectric constants, ξhBN ≡ (ϵ⊥hBN∕ϵ∥hBN)1∕2 is the hBN anisotropy factor, ϵSiO2 = 3.9 is the dc dielectric constant of SiO2, ϵSi = 11.7 is the dc dielectric con-stant of Si, and ϵ0 = 1 is the dielectric constant of vacuum. Parame-ters d1 = 16  nm and d2 = 55  nm are thicknesses of the top and bottom hBN layers, respectively, and d3 = 285  nm is the SiO2 thickness. Function κg(q) specified by Eq.  19c accounts for the (static) screening of the Coulomb potential by electrons in gra-phene and σxx(q, ω) is the longitudinal conductivity of graphene in the presence of magnetic field computed as described in section Numerical calculation of graphene magneto-optical conductivity. Representative plots of functions κ(q), κg(q), and ϵ(q) are shown in Fig. 10 (A and B).As one can see from these graphs, κ(q) (the middle curve) is equal to (ϵ0 + ϵSi)/2 at zero q, then has a small dip to approach (ϵ0 + ϵSiO2)/2 within a narrow range of relatively low momenta d−13≲ q ≲ d−12 , and then rises and tends to ϵhBN at q ≳ d−11= 6.7 × 105 cm−1 . The total dielectric function ϵ(q) (the top curve) shows the same small dip at low q, goes through a modest maximum, and then approaches the limiting valueat high momenta q ≫ l−1B.The calculation of the binding energies Δ−n,n� in Eq. 19 involves numerical evaluation of four integrals of the formwhich are well behaved. On the other hand, the integral for the self-energy in Eq. 18 diverges at large momenta qlB ≫√n,√n′ becauseas can be deduced from Eqs. 4 and 18. We regularize this divergence by renormalization, i.e., subtraction of vF,−n→n� evaluated at some reference field B0. Let l0 ≡ lB(B0) be the magnetic length at B0. Since the high momenta enter through the product qlB, we can rescale the ϵ(q)= κ(q)+ κg(q)(20a)κ(q)=ϵhBN2(ϵhBNtanhξhBNqd1+ϵ0ϵ0tanhξhBNqd1+ϵhBN+ϵhBNtanhξhBNqd2+ϵSOSϵSOStanhξhBNqd2+ϵhBN)ϵSOS=ϵSiO2ϵSiO2tanhqd3+ϵSiϵSitanhqd3+ϵSiO2(20b)κg(q)= 2πiq limω→ 0ω−1σxx(q,ω)(20c)ϵ∞ = ϵhBN +πα2, α =e2ℏvrenF≈ 2.2 (21)∫d2q(2π)2v(q)Ln1(q2l2B2)Ln2(q2l2B2)e−q2 l2B∕2 (22)I−n,n��q�≃ −√2n +√2n�4qlB(23)Fig. 10. Calculation of the Fermi velocities and many-body effects. (A and B) Effective dielectric function for electrons in graphene. Functions κ(q), κg(q), and ϵ(q) de-fined by Eq. 19 (A to C) for B = 3.35 T. (B) shows a magnification of the small dip at low q that is observed for κg(q) and ϵ(q). (C) Renormalized Fermi velocities given by Eq. 26. (D) Effective Fermi velocities defined by Eq. 17 that include excitonic corrections.Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e14 of 16integration variable in Eq. 18 by the ratio R = lB ∕ l0 =√B0 ∕B to can-cel the divergence. Performing the subtraction, we findwhereNote that if the dielectric function is replaced by a constant, e.g., ϵ∞, then v(q) = 2πe2/(ϵ∞q) and ΔvF,−n→n�(B, B0) vanishes identically. In this case, the conventional logarithmic-in-B rule (37) for the Fermi velocity renormalization is recovered.In our experiment, all the ILTs were measured at the same pho-ton energy ℏω. Therefore, it is convenient to select a particular ILT, e.g., −2 → 3, as a reference, so that B0 = B−2→3, where B−2→3 is the field at which this ILT occurs. Equation 24 entailsOur calculations using this formula show that vF,−n→n� is almost the same for all the ILTs, as shown in Fig. 10C. On the other hand, veffF,−n→n� , given by Eq. 17, exhibits characteristic dips at n = n′ be-cause of relatively larger excitonic corrections Δ−n,n� at such ILTs; see Fig. 10D.To obtain the experimental binding energies Δexp,−n,n�, we use Eqs. 17 and 26 but set the effective Fermi velocities equal to their ex-perimentally measured values veffF,exp,m→n . Solving for Δ−n,n�, we obtainThis gives Δexp,−n,n� in terms of veffF,exp,−n→n� and veffF,exp,−2→3 , the cal-culated electron self-energies, and the calculated binding energy Δ−2,3 for the reference ILT.Processing of experimental near-field dataThe following section describes the processing of the experimental near-field data that are displayed in Figs. 2 (D to F), 3A, and 5A. We want to emphasize that we tried to keep data and image processing to a minimum. The data processing that we did apply is described below for each (sub-)figure showing experimental near-field data:Figure  2 (D to F): The 2D scans at 0.0 T (Fig.  2D) and 3.35 T (Fig.  2E) show raw data. No filter was applied. The scan at 6.0 T (Fig. 2F) was corrected for horizontal strokes using the “correct hori-zontal scars (strokes)” function of gwyddion (program version 2.62).Figure  3A: The near-field data were corrected for slow spatial drift during this linescan using the topography information record-ed simultaneously with the near-field data. Each line was horizon-tally shifted such that the topography step of the gold electrode aligns for all lines (refer to Fig. 1D for a sketch of the sample). The topography information was also used to select the region without the gold electrode. Using this information derived from the topog-raphy data, the same procedure was then applied to the near-field data to obtain the drift-corrected version of the near-field response on hBN and hBN-graphene-hBN shown in Fig. 3A. In addition, we found that the average near-field signal of our measurement slowly decreased with time. This can be assigned to minor drift in the opti-cal alignment over the long measurement time of >8 hours. Assum-ing a linear dependence on time, we estimated a drift of 1.01% signal reduction per hour of measurement time relative to the initial near-field signal. We multiplied the near-field signal of each line with a constant compensating for this slow linear drift.Figure  5A: The near-field data were corrected for slow spatial drift during this line scan using the topography information record-ed simultaneously to the near-field data; this procedure is identical to the one applied for the near-field data of Fig. 3A. We did not com-pensate for any potential temporal drift of the near-field signal, as we assume that it will be negligible because of the shorter measure-ment time compared to Fig. 3A.Extracting the polariton wavelength and quality factor from near-field line profilesEach horizontal line of the near-field data displayed in Figs. 3A and 5A represents a line profile such as the examples displayed in Figs. 3C and 5B. We analytically describe each line profile via an exponentially damped sine wave (20, 68–70) with linear back-ground termThe vertical offset S0, the polariton wavelength λp, the polariton decay length Lp, and the phase offset φ0 are optimized for each line using a least-square fitting procedure. For this fitting procedure, the left-most 10% of each line are neglected. The linear background slope S1 and the amplitude A are constants that are optimized once for the whole dataset. However, S1 and A are not part of the fitting procedure and, thus, have the same value for all line profiles within a dataset.Applying this fitting procedure to the data displayed in Figs. 3A (Fig.  5A) directly yields the polariton wavelength λp displayed in Figs. 4A and 5C. The quality factor Q (Figs. 4B and 5D) is calculated as Q = 2πLp/λp, i.e., it directly follows from the ratio of the fit param-eters Lp and λp.For our analytical description, we assume that our observed fringe pattern is dominated by polariton waves that are launched by the gold electrode (refer to Figs. 1D, and 2 D to F). This is opposed to polaritons launched by the tip and then reflected by the gold electrode.Note that our analytical description above may be modified to reflect the geometrical spreading of the polariton waves, resulting in a function of the following formThe additional factor x−p reflects the amplitude reduction due to the geometrical spreading; the most common version of this factor uses a power p = 1/2, i.e., a factor x−1/2 corresponding to circularly spreading waves in a plane (20, 68–70). If we apply this modified (24)(25)vF,−n→n� = vF,−2→3−e28ℏϵ∞ln�BB−2→3�−lB√2ℏ ∫d2q(2π)2�v�q� I−n,n��q�√n+√n�−Rv�Rq� I−2,3�q�√2+√3� (26)(27)S(x) = S0 + S1x + Ae−x∕Lpsin(2πxλp− φ0)(28)S(x) = S0 + S1x + Ax−p e−x∕Lpsin(2πxλp− φ0)(29)Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e15 of 16description to our data, we still obtain a good match to our data, but the resulting polariton quality factors are about three times larger compared to the values shown in the main text. However, while it is common to apply the geometrical spreading factor for polaritons launched by the tip or small antennas (20, 68–70), in our case, the edge of a gold contact (refer to Fig. 1A for the sample design) repre-sents a linear launcher with extension larger than the analyzed polariton propagation length (70). This is why we chose not to add the fac-tor x−1/2 in our analysis. Nevertheless, this also suggests that the quality factors shown in Figs. 4B and 5D may be an underestimation.REFERENCES AND NOTES  1. T . Low, A. Chaves, J. D. Caldwell, A. Kumar, N. X. Fang, P. Avouris, T. F. Heinz, F. Guinea,  L. Martin-Moreno, F. Koppens, Polaritons in layered two-dimensional materials. Nat. Mater. 16, 182–194 (2017).  2.  Q. Zhang, G. Hu, W. Ma, P. Li, A. Krasnok, R. Hillenbrand, A. Alù, C.-W. Qiu, Interface nano-optics with van der Waals polaritons. Nature 597, 187–195 (2021).  3. D . N. Basov, M. M. Fogler, F. J. Garcia de Abajo, Polaritons in van der Waals materials. Science 354, aag1992 (2016).  4.  F. Dirnberger, J. Quan, R. Bushati, G. M. Diederich, M. Florian, J. Klein, K. Mosina, Z. Sofer,  X. Xu, A. Kamra, F. J. García-Vidal, A. Alù, V. M. Menon, Magneto-optics in a van der Waals magnet tuned by self-hybridized polaritons. Nature 620, 533–537 (2023).  5. D . N. Basov, A. Asenjo-Garcia, P. J. Schuck, X. Zhu, A. Rubio, Polariton panorama. Nanophotonics 10, 549–577 (2020).  6.  M. Dapolito, M. Tsuneto, W. Zheng, L. Wehmeier, S. Xu, X. Chen, J. Sun, Z. Du, Y. Shao,  R. Jing, S. Zhang, A. Bercher, Y. Dong, D. Halbertal, V. Ravindran, Z. Zhou, M. Petrovic,  A. Gozar, G. L. Carr, Q. Li, A. B. Kuzmenko, M. M. Fogler, D. N. Basov, X. Du, M. Liu, Infrared nano-imaging of Dirac magnetoexcitons in graphene. Nat. Nanotechnol. 18, 1409–1415 (2023).  7. T . M. Slipchenko, J.-M. Poumirol, A. B. Kuzmenko, A. Y. Nikitin, L. Martín-Moreno, Interband plasmon polaritons in magnetized charge-neutral graphene. Commun. Phys. 4, 110 (2021).  8.  S. Dai, Z. Fei, Q. Ma, A. S. Rodin, M. Wagner, A. S. McLeod, M. K. Liu, W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M. Thiemens, G. Dominguez, A. H. C. Neto, A. Zettl, F. Keilmann, P. Jarillo-Herrero, M. M. Fogler, D. N. Basov, Tunable phonon polaritons in atomically thin van der waals crystals of boron nitride. Science 343, 1125–1129 (2014).  9.  J. D. Caldwell, A. V. Kretinin, Y. Chen, V. Giannini, M. M. Fogler, Y. Francescato, C. T. Ellis,  J. G. Tischler, C. R. Woods, A. J. Giles, M. Hong, K. Watanabe, T. Taniguchi, S. A. Maier,  K. S. Novoselov, Sub-diffractional volume-confined polaritons in the natural hyperbolic material hexagonal boron nitride. Nat. Commun. 5, 5221 (2014).  10. E . Yoxall, M. Schnell, A. Y. Nikitin, O. Txoperena, A. Woessner, M. B. Lundeberg, F. Casanova, L. E. Hueso, F. H. L. Koppens, R. Hillenbrand, Direct observation of ultraslow hyperbolic polariton propagation with negative phase velocity. Nat. Photonics 9, 674–678 (2015).  11. T . V. A. G. Oliveira, T. Nörenberg, G. Álvarez-Pérez, L. Wehmeier, J. Taboada-Gutiérrez,  M. Obst, F. Hempel, E. J. H. Lee, J. M. Klopf, I. Errea, A. Y. Nikitin, S. C. Kehr,  P. Alonso-González, L. M. Eng, Nanoscale-confined terahertz polaritons in a van der Waals crystal. Adv. Mater. 33, e2005777 (2021).  12.  A. Fali, S. T. White, T. G. Folland, M. He, N. A. Aghamiri, S. Liu, J. H. Edgar, J. D. Caldwell,  R. F. Haglund, Y. Abate, Refractive index-based control of hyperbolic phonon-polariton propagation. Nano Lett. 19, 7725–7734 (2019).  13.  Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. McLeod, M. Wagner, L. M. Zhang, Z. Zhao,  M. Thiemens, G. Dominguez, M. M. Fogler, A. H. C. Neto, C. N. Lau, F. Keilmann, D. N. Basov, Gate-tuning of graphene plasmons revealed by infrared nano-imaging. Nature 487, 82–85 (2012).  14.  Z. Fei, G. O. Andreev, W. Bao, L. M. Zhang, A. S. McLeod, C. Wang, M. K. Stewart, Z. Zhao,  G. Dominguez, M. Thiemens, M. M. Fogler, M. J. Tauber, A. H. Castro-Neto, C. N. Lau,  F. Keilmann, D. N. Basov, Infrared nanoscopy of dirac plasmons at the graphene–SiO2 interface. Nano Lett. 11, 4701–4705 (2011).  15.  J. Chen, M. Badioli, P. Alonso-González, S. Thongrattanasiri, F. Huth, J. Osmond,  M. Spasenović, A. Centeno, A. Pesquera, P. Godignon, A. Z. Elorza, N. Camara,  F. J. G. de Abajo, R. Hillenbrand, F. H. L. Koppens, Optical nano-imaging of gate-tunable graphene plasmons. Nature 487, 77–81 (2012).  16.  Y. Dong, L. Xiong, I. Y. Phinney, Z. Sun, R. Jing, A. S. McLeod, S. Zhang, S. Liu, F. L. Ruta,  H. Gao, Z. Dong, R. Pan, J. H. Edgar, P. Jarillo-Herrero, L. S. Levitov, A. J. Millis, M. M. Fogler,  D. A. Bandurin, D. N. Basov, Fizeau drag in graphene plasmonics. Nature 594, 513–516 (2021).  17.  W. Zhao, S. Zhao, H. Li, S. Wang, S. Wang, M. I. B. Utama, S. Kahn, Y. Jiang, X. Xiao, S. Yoo,  K. Watanabe, T. Taniguchi, A. Zettl, F. Wang, Efficient fizeau drag from dirac electrons in monolayer graphene. Nature 594, 517–521 (2021).  18. C . H. Henry, J. J. Hopfield, Raman scattering by polaritons. Phys. Rev. Lett. 15, 964–966 (1965).  19.  S. Dai, Q. Ma, M. K. Liu, T. Andersen, Z. Fei, M. D. Goldflam, M. Wagner, K. Watanabe,  T. Taniguchi, M. Thiemens, F. Keilmann, G. C. A. M. Janssen, S.-E. Zhu, P. Jarillo-Herrero,  M. M. Fogler, D. N. Basov, Graphene on hexagonal boron nitride as a tunable hyperbolic metamaterial. Nat. Nanotechnol. 10, 682–686 (2015).  20.  A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega,  K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillenbrand, F. H. L. Koppens, Highly confined low-loss plasmons in graphene–boron nitride heterostructures. Nat. Mater. 14, 421–425 (2015).  21.  A. T. Costa, P. A. D. Gonçalves, D. N. Basov, F. H. L. Koppens, N. A. Mortensen, N. M. R. Peres, Harnessing ultraconfined graphene plasmons to probe the electrodynamics of superconductors. Proc. Natl. Acad. Sci. U.S.A. 118, e2012847118 (2021).  22.  G. X. Ni, A. S. McLeod, Z. Sun, L. Wang, L. Xiong, K. W. Post, S. S. Sunku, B.-Y. Jiang, J. Hone, C. R. Dean, M. M. Fogler, D. N. Basov, Fundamental limits to graphene plasmonics. Nature 557, 530–533 (2018).  23.  S. Das Sarma, Theory of two-dimensional magnetopolarons. Phys. Rev. Lett. 52, 859–862 (1984).  24.  J. Zhu, S. M. Badalyan, F. M. Peeters, Electron-phonon bound states in graphene in a perpendicular magnetic field. Phys. Rev. Lett. 109, 256602 (2012).  25. E . J. Johnson, D. M. Larsen, Polaron induced anomalies in the interband magnetoabsorption of InSb. Phys. Rev. Lett. 16, 655–659 (1966).  26. C . Faugeras, M. Orlita, M. Potemski, Raman scattering of graphene-based systems in high magnetic fields. J. Raman Spectrosc. 49, 146–156 (2018).  27. D . C. Tsui, T. Englert, A. Y. Cho, A. C. Gossard, Observation of magnetophonon resonances in a two-dimensional electronic system. Phys. Rev. Lett. 44, 341–344 (1980).  28. L . V. Butov, V. I. Grinev, V. D. Kulakovskii, T. G. Andersson, Direct observation of magnetoplasmon-phonon coupled modes in the magnetophotoluminescence spectra of the two-dimensional electron gas in InxGa1-xAs/GaAs quantum wells. Phys. Rev. B 46, 13627–13630 (1992).  29.  X. Li, M. Bamba, Q. Zhang, S. Fallahi, G. C. Gardner, W. Gao, M. Lou, K. Yoshioka,  M. J. Manfra, J. Kono, Vacuum Bloch–Siegert shift in Landau polaritons with ultra-high cooperativity. Nat. Photon. 12, 324–329 (2018).  30.  F. Keilmann, D. W. van der Weide, T. Eickelkamp, R. Merz, D. Stöckle, Extreme sub-wavelength resolution with a scanning radio-frequency transmission microscope. Opt. Commun. 129, 15–18 (1996).  31. H . U. Yang, E. Hebestreit, E. E. Josberger, M. B. Raschke, A cryogenic scattering-type scanning near-field optical microscope. Rev. Sci. Instrum. 84, 023701 (2013).  32.  R. H. J. Kim, J.-M. Park, S. J. Haeuser, L. Luo, J. Wang, Cryogenic magneto-terahertz scanning near-field optical microscope (cm-SNOM). arXiv:2210.07319 [cond-mat.mes-hall] (2022).  33. I . O. Nedoliuk, S. Hu, A. K. Geim, A. B. Kuzmenko, Colossal infrared and terahertz magneto-optical activity in a two-dimensional Dirac material. Nat. Nanotechnol. 14, 756–761 (2019).  34. C . Faugeras, S. Berciaud, P. Leszczynski, Y. Henni, K. Nogajewski, M. Orlita, T. Taniguchi,  K. Watanabe, C. Forsythe, P. Kim, R. Jalil, A. K. Geim, D. M. Basko, M. Potemski, Landau  level spectroscopy of electron-electron interactions in graphene. Phys. Rev. Lett. 114,  126804 (2015).  35.  B. J. Russell, B. Zhou, T. Taniguchi, K. Watanabe, E. A. Henriksen, Many-particle effects in the cyclotron resonance of encapsulated monolayer graphene. Phys. Rev. Lett. 120, 047401 (2018).  36. D . N. Basov, M. M. Fogler, A. Lanzara, F. Wang, Y. Zhang, Colloquium: Graphene spectroscopy. Rev. Mod. Phys. 86, 959–994 (2014).  37.  K. Shizuya, Many-body corrections to cyclotron resonance in monolayer and bilayer graphene. Phys. Rev. B 81, 075407 (2010).  38.  M. L. Sadowski, G. Martinez, M. Potemski, C. Berger, W. A. de Heer, Landau level spectroscopy of ultrathin graphite layers. Phys. Rev. Lett. 97, 266405 (2006).  39.  Z. Jiang, E. A. Henriksen, L. C. Tung, Y.-J. Wang, M. E. Schwartz, M. Y. Han, P. Kim,  H. L. Stormer, Infrared spectroscopy of landau levels of graphene. Phys. Rev. Lett. 98, 197403 (2007).  40. D . C. Elias, R. V. Gorbachev, A. S. Mayorov, S. V. Morozov, A. A. Zhukov, P. Blake,  L. A. Ponomarenko, I. V. Grigorieva, K. S. Novoselov, F. Guinea, A. K. Geim, Dirac cones reshaped by interaction effects in suspended graphene. Nat. Phys. 7, 701–704 (2011).  41.  J. González, F. Guinea, M. A. H. Vozmediano, Non-fermi liquid behavior of electrons in the half-filled honeycomb lattice (A renormalization group approach). Nucl. Phys. B. 424, 595–618 (1994).  42. I . Crassee, J. Levallois, A. L. Walter, M. Ostler, A. Bostwick, E. Rotenberg, T. Seyller,  D. van der Marel, A. B. Kuzmenko, Giant Faraday rotation in single- and multilayer graphene. Nat. Phys. 7, 48–51 (2011).  43.  R. Roldán, J.-N. Fuchs, M. O. Goerbig, Collective modes of doped graphene and a standard two-dimensional electron gas in a strong magnetic field: Linear magnetoplasmons versus magnetoexcitons. Phys. Rev. B. 80, 085408 (2009).  44.  A. Bylinkin, M. Schnell, M. Autore, F. Calavalle, P. Li, J. Taboada-Gutièrrez, S. Liu, J. H. Edgar, F. Casanova, L. E. Hueso, P. Alonso-Gonzalez, A. Y. Nikitin, R. Hillenbrand, Real-space observation of vibrational strong coupling between propagating phonon polaritons and organic molecules. Nat. Photonics 15, 197–202 (2021).Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024Wehmeier et al., Sci. Adv. 10, eadp3487 (2024)     13 September 2024S c i e n c e  A d v a n c e s  |  R e s e ar  c h  A r t i c l e16 of 1645. L . Novotny, Strong coupling, energy splitting, and level crossings: A classical perspective. Am. J. Phys. 78, 1199–1202 (2010).46. C . Neumann, S. Reichardt, M. Drögeler, B. Terrés, K. Watanabe, T. Taniguchi, B. Beschoten, S. V. Rotkin, C. Stampfer, Low B field magneto-phonon resonances in single-layer and bilayer graphene. Nano Lett. 15, 1547–1552 (2015).47. R. J. Nicholas, The magnetophonon effect. Prog. Quant. Electron. 10, 1–75 (1985).48. Y. Tokura, N. Nagaosa, Nonreciprocal responses from non-centrosymmetric quantum materials. Nat. Commun. 9, 3740 (2018).49. T . Van Mechelen, W. Sun, Z. Jacob, Optical N-invariant of graphene’s topological viscous Hall fluid. Nat. Commun. 12, 4729 (2021).50. A. Arora, M. S. Rudner, J. C. W. Song, Quantum plasmonic nonreciprocity in parity-violating magnets. Nano Lett. 22, 9351–9357 (2022).51. F. Appugliese, J. Enkner, G. L. Paravicini-Bagliani, M. Beck, C. Reichl, W. Wegscheider, G. Scalari, C. Ciuti, J. Faist, Breakdown of topological protection by cavity vacuum fields in the integer quantum Hall effect. Science 375, 1030–1034 (2022).52. P. Sivarajah, A. Steinbacher, B. Dastrup, J. Lu, M. Xiang, W. Ren, S. Kamba, S. Cao, K. A. Nelson, THz-frequency magnon-phonon-polaritons in the collective strong-coupling regime. J. Appl. Phys. 125, 213103 (2019).53. V . Rokaj, M. Penz, M. A. Sentef, M. Ruggenthaler, A. Rubio, Polaritonic Hofstadter butterfly and cavity control of the quantized hall conductance. Phys. Rev. B 105, 205424 (2022).54. L . Winter, O. Zilberberg, Fractional quantum Hall edge polaritons. arXiv:2308.12146 [cond-mat.mes-hall] (2023).55. M. M. Fogler, A. A. Koulakov, B. I. Shklovskii, Ground state of a two-dimensional electron liquid in a weak magnetic field. Phys. Rev. B 54, 1853–1871 (1996).56. M. Dapolito, X. Chen, C. Li, M. Tsuneto, S. Zhang, X. Du, M. Liu, A. Gozar, Scattering-type scanning near-field optical microscopy with Akiyama piezo-probes. Appl. Phys. Lett. 120, 013104 (2022).57. X. Chen, D. Hu, R. Mescall, G. You, D. N. Basov, Q. Dai, M. Liu, Modern scattering-type scanning near-field optical microscopy for advanced material research. Adv. Mater. 31, e1804774 (2019).58. B. Knoll, F. Keilmann, Enhanced dielectric contrast in scattering-type scanning near-field optical microscopy. Opt. Commun. 182, 321–328 (2000).59. A. K. Geim, I. V. Grigorieva, Van der Waals heterostructures. Nature 499, 419–425 (2013).60. L . Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe,L. M. Campos, D. A. Muller, J. Guo, P. Kim, J. Hone, K. L. Shepard, C. R. Dean, One-dimensionalelectrical contact to a two-dimensional material. Science 342, 614–617 (2013).61. U. Briskot, I. A. Dmitriev, A. D. Mirlin, Quantum magneto-oscillations in the ac conductivity of disordered graphene. Phys. Rev. B 87, 195432 (2013).62. B.-Y. Jiang, L. M. Zhang, A. H. Castro Neto, D. N. Basov, M. M. Fogler, Generalized spectral method for near-field optical microscopy. J. Appl. Phys. 119, 054305 (2016).63. J. Aizpurua, T. Taubner, F. J. García de Abajo, M. Brehm, R. Hillenbrand, Substrate-enhanced infrared near-field spectroscopy. Opt. Express 16, 1529–1545 (2008).64. J. Renger, S. Grafström, L. M. Eng, R. Hillenbrand, Resonant light scattering by near-field-induced phonon polaritons. Phys. Rev. B. 71, 075410 (2005).65. L . Wehmeier, D. Lang, Y. Liu, X. Zhang, S. Winnerl, L. M. Eng, S. C. Kehr, Polarization-dependent near-field phonon nanoscopy of oxides: SrTiO3, LiNbO3, and PbZr0.2Ti0.8O3. Phys. Rev. B 100, 035444 (2019).66. I . Dolado, C. Maciel-Escudero, E. Nikulina, E. Modin, F. Calavalle, S. Chen, A. Bylinkin, F. J. Alfaro-Mozaz, J. Li, J. H. Edgar, F. Casanova, S. Vélez, L. E. Hueso, R. Esteban, J. Aizpurua, R. Hillenbrand, Remote near-field spectroscopy of vibrational strong coupling between organic molecules and phononic nanoresonators. Nat. Commun. 13, 6850 (2022).67.  X. Wu, S. K. Gray, M. Pelton, Quantum-dot-induced transparency in a nanoscale plasmonic resonator. Opt. Express. OE 18, 23633–23645 (2010).68.  W. Ma, P. Alonso-González, S. Li, A. Y. Nikitin, J. Yuan, J. Martín-Sánchez, J. Taboada-Gutiérrez, I. Amenabar, P. Li, S. Vélez, C. Tollan, Z. Dai, Y. Zhang, S. Sriram, K. Kalantar-Zadeh, S.-T. Lee, R. Hillenbrand, Q. Bao, In-plane anisotropic and ultra-low-loss polaritons in a natural van der Waals crystal. Nature 562, 557–562 (2018).69.  P. Li, I. Dolado, F. J. Alfaro-Mozaz, F. Casanova, L. E. Hueso, S. Liu, J. H. Edgar, A. Y. Nikitin, S. Vélez, R. Hillenbrand, Infrared hyperbolic metasurface based on nanostructured van der Waals materials. Science 359, 892–896 (2018).70.  R. A. Mayer, F. H. Feres, F. C. B. Maia, I. D. Barcelos, A. S. McLeod, A. Rodin, R. O. Freitas, Guidelinesfor engineering directional polariton launchers. Phys. Rev. Appl. 18, 034089 (2022).Acknowledgments: We are grateful to R. Freitas (Brazilian Synchrotron light laboratory, lnlS),  A. Sternbach, and F. Ruta (both columbia University) for the helpful discussion, as well as for the helpful discussion and technical support from J. li, h. Wang, and W. Wang. this was produced under contract no. deSc0012704 with the US department of energy. the US Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this paper, or allow others to do so, for US Government purposes. dOe will provide public access (http://energy.gov/downloads/doe-public-access-plan). distributed under a creative commons Attribution noncommercial license 4.0 (cc BY-nc). Funding: Research on m- SnOM scanning probe platform is supported as part of Programmable Quantum Materials, an energy Frontier Research center funded by the U.S. department of energy (dOe), Office of Science, Basic energy Sciences (BeS), under award de- Sc0019443. l.W., X.c., M.l., and d.n.B. acknowledge U.S. department of energy, Office of Science, national Quantum information Science Research centers, co- design center for Quantum Advantage (c2QA) under contract number de- Sc0012704 for the support of the data analysis. M.l. acknowledges support from the nSF Faculty early career development Program under grant no. dMR- 2045425 for the development of the Akiyama probe–based SnOM. M.l. also acknowledges Gordon and Betty Moore Foundation dOi: 10.37807/gbmf12258 for supporting the development of polaritonic materials. M.d., R.J., M.l., and Q.l. acknowledge U.S. department of energy, Office of Basic energy Sciences, division of Materials Sciences and engineering, under contract no. de- Sc0012704 for supporting the sample characterization and theory development. R.A.M. thanks grant 2022/06709- 0, São Paulo Research Foundation (FAPeSP) for modeling development. K.W. and t.t. acknowledge support from the JSPS KAKenhi (grant numbers 20h00354, 21h05233. and 23h02052) and World Premier international Research center initiative (WPi), MeXt, Japan. A.B.K. is supported by the Swiss national Science Foundation under grant no. 200020_201096. X.d. acknowledges support from the nSF under awards dMR- 1808491 for sample fabrication. Author contributions: l.W., d.n.B., and M.l. conceived the project and designed the experiments. Q.l., G.l.c., X.d., M.M.F., d.n.B., and M.l. supervised the project. R.P., K.W., t.t., and X.d. prepared the devices. l.W., M.t., Z.d., and M.l. performed the experimental measurements with support from M.d., W.Z., and Z.Z. l.W. and M.l. analyzed the experimental data with input from S.X., R.A.M., R.J., A.G., A.B.K., G.l.c., M.M.F., and d.n.B. B.v. and M.M.F. developed the theoretical description. S.X. and R.A.M. developed the numerical simulation with input from l.W., X.c., R.J., M.M.F., and M.l. l.W., M.M.F., and M.l. cowrote the manuscript with input from all co- authors. Competing interests: M.d., X.c., A.G., and M.l. have filed a patent regarding magneto scanning near- field optical microscopy (publication number: WO/2023/049225; publication date: 30 March 2023; international application no.: Pct/US2022/044311; international filing date: 22 September 2022; organizations issuing patent: the Research Foundation For the State University Of new York, 35 State Street, Albany, new York 12207, USA; Yale University, two Whitney Ave, new haven, connecticut 06511, USA). the other authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper.Submitted 20 March 2024 Accepted 7 August 2024 Published 13 September 2024 10.1126/sciadv.adp3487Downloaded from https://www.science.org at National Institute for Materials Science on December 19, 2024http://energy.gov/downloads/doe-public-access-planhttp://dx.doi.org/10.37807/gbmf12258 Landau-phonon polaritons in Dirac heterostructures INTRODUCTION RESULTS High-momentum magneto-optics of graphene Modeling of polariton dispersion Nano-imaging of LPP Tunability of LPPs Many-body effects DISCUSSION MATERIALS AND METHODS Experimental setup: Magneto infrared nanoscopy Fabrication of graphene-hBN heterostructure Pauli-blocking and gate-switchable dispersion Numerical calculation of graphene magneto-optical conductivity, polariton dispersion, and near-field signal: (figs. 1, C and D, 2D, and 3D) Momentum-range accessible via m-SNOM Extracting the coupling strength via the coupled harmonic oscillator model Coupled harmonic oscillators model for hBN phonon polaritons and graphene Landau polaritons Extraction of the oscillator parameters for uncoupled modes Dispersion of branches and coupling strength extraction Coupling strength and criterium for strong coupling Weak polariton fringes parallel to graphene boundary in fig. 2 (D to F) Many-body effects: Fermi velocity and magnetoexciton binding energy Processing of experimental near-field data Extracting the polariton wavelength and quality factor from near-field line profiles REFERENCES AND NOTES Acknowledgments