Remote near-field spectroscopy of vibrational strong coupling between organic molecules and phononic nanoresonators

Article https://doi.org/10.1038/s41467-022-34393-4 Remote near-field spectroscopy of vibrational strong coupling between organic molecules and phononic nanoresonators Irene Dolado 1,8 , Carlos Maciel-Escudero 1,2,8 , Elizaveta Nikulina 1 , Evgenii Modin 1 ,FrancescoCalavalle 1 ,ShuChen 1 ,AndreiBylinkin 1,3 , Francisco Javier Alfaro-Mozaz 3 ,JiahanLi 4 ,JamesH.Edgar 4 , Fèlix Casanova 1,5 ,Sau l Vélez 6 ,LuisE.Hueso 1,5 ,RubenEsteban 2,3 , Javier Aizpurua 2,3 & Rainer Hillenbrand 5,7 Phonon polariton (PhP) nanoresonators can dramatically enhance the coupling of molecular vibrations and infrared light, enabling ultrasensitive spectroscopies and strong coupling with minute amounts of matter. So far, this coupling and the resulting localized hybrid polariton modes have been studied only by far-field spectroscopy, preventing access to modal near-field patterns and dark modes, which could further our fundamental understanding of nanoscale vibrational strong coupling (VSC). Here we use infrared near-field spectroscopy to study the coupling between the localized modes of PhP nanoresonators made of h-BN and molecular vibrations. For a most direct probing of the resonator-molecule coupling, we avoid the direct near-field interaction between tip and molecules by probing the molecule-free part of partially molecule-covered nanoresonators, which we refer to as remote nearfield probing. We obtain spatially and spectrally resolved maps of the hybrid polariton modes, as well as the corresponding coupling strengths, demonstrating VSC on a single PhP nanoresonator level. Our work paves the way for near-field spectroscopy of VSC phenomena not accessible by conventional techniques. Strong coupling between molecular vibrations and infrared photons (vibrationalstrong coupling, VSC) leads to hybrid light-matter states1–6. They offer intriguing possibilities for ultra-sensitive vibrational spectroscopy7–10 and for modifying chemical reactions3,4. Typically, VSC is achieved with molecules embedded into microcavities, implying large mode volumes and large amounts of molecules, which limits access to quantum phenomena that may be accesible only for nanoscale amounts of molecules or at the level of a few molecules. In this regard, plasmonic infrared resonators are a promising route to achieve VSC at the nanoscale7,8owing to their dramatically reduced mode volumes as compared to microcavities. Alternatively, phonon polaritons (PhP) can be employed for VSC experiments, offering stronger polariton confinement and larger quality factors9–11. Unfortunately, the far-field extinction cross-section of individual PhP Received: 18 July 2022 Accepted: 21 October 2022 Check for updates 1 CIC nanoGUNE BRTA, 20018 Donostia-San Sebastián, Spain. 2 Materials Physics Center, CSIC-UPV/EHU, 20018 Donostia-San Sebastián, Spain. 3 Donostia International Physics Center (DIPC), Donostia-San Sebastián, Spain. 4 Tim Taylor Department of Chemical Engineering, Kansas State University, Manhattan, KS 66506, USA. 5 IKERBASQUE, Basque Foundation for Science, 48009 Bilbao, Spain. 6 IFIMAC-Condensed Matter Physics Center, Instituto Nicolás Cabrera, and Departamento de Física de la Materia Condensada, Universidad Autónoma de Madrid, Madrid E-28049, Spain. 7 CIC nanoGUNE BRTA and Department of Electricity and Electronics, UPV/EHU, 20018 Donostia-San Sebastián, Spain. 8 These authors contributed equally: Irene Dolado, Carlos Maciel-Escudero. e-mail: r.hi[email protected] Nature Communications | (2022) 13:6850 1 1234567890():,; 1234567890():,; nanoresonators12–20 is extremely small (due to their small size compared to the infrared wavelength), challenging infrared far-field spectroscopy. Further, subradiant dark modes—offering the advantage of longer lifetimes—are difficult to probe by far-field spectroscopy. These problems can be circumvented by nanoscale Fourier Transform Infrared (nano-FTIR) spectroscopy, which employs the strong field concentration at the apex of a metallic scanning probe tip (the nearfield probe)21,22 to enable near-field spectroscopy and spatial mapping of both bright and dark modes of individual phononic nanoresonators17,18. However, the near-field mapping of VSC employing PhP nanoresonators has been elusive so far. Nano-FTIR spectroscopy has been employed to study the coupling between molecular vibrations and plasmonic resonators, but the near-field probe itself can also couple with the plasmonic resonator and the molecules7,23,24, eventually reaching strong tip-resonator and tip-molecule coupling. Although this coupling may be exploited for on-demand control of VSC, it may challenge the probing of the hybrid polariton modes that are exclusively formed by the moleculenanoresonator coupling. In ref. 17, nano-FTIR spectroscopy was applied to study strong coupling between molecular vibrations and propagating PhPs on extended, unstructured h-BN layers. In this experiment, polariton interferometry had to be applied to obtain the polariton momenta and thus the polariton dispersion. This measurement principle and its associated data analysis offers the advantage that the tip-molecule and tip-polariton coupling does not affect the polariton dispersion and analysis of VSC. However, such experiments do not allow for analyzing the coupling between molecular vibrations and localized PhP modes, which will be of utmost importance to experimentally explore, test and verify future VSC concepts employing individual PhP nanoresonators. Here we demonstrate that hybrid polariton modes formed by vibrational strong coupling between a single PhP nanoresonator and molecular vibrations can be studied and imaged in real space by nanoFTIR spectroscopy. We minimize the influence of the tip by probing the molecule-free part of partially molecule-covered PhP nanoresonators with a non-resonant metallic tip, which we refer to asremote near-field probing. We also verify our experimental results via comparative numerical simulations, where the near-field probe is modeled either by a point-dipole source (representing a non-disturbing nearfield probe) or by an oscillating metal tip (representing a more realistic and potentially disturbing near-field probe). Results Nano-FTIR spectroscopy of half-covered phononic nanoresonators Our experiment is illustrated in Fig. 1a. Half of a hexagonal boron nitride (h-BN) nanorod is covered with CBP molecules (4,4′-bis(N-carbazolyl)−1,1′-biphenyl; organic semiconductor) that exhibit a vibrational resonance at ω CBP = 1450 cm−1(Supplementary Note 1), as can be recognized from the near-field spectrum of a thin CBP layer shown by the green curve in Fig. 1d, e. The near-field probe, the non-resonant metallic tip25 of an atomic force microscope, is placed remotely with respect to the molecules at the opposite extremity of the h-BN rod. The tip (oscillating normal to the sample at frequency Ω) concentrates an illuminating broadband infrared laser beam at its apex to a nanoscale near-field spot, which excites phonon polaritons26–28 exhibiting FabryPerot (FP) resonances9,18,19 in the h-BN nanorod17.Figure1bshowsa simulation of the second-order FP mode excited by the near fields of a point dipole (note that this mode is often referred to as a dark mode, since it cannot be excited by far-field illumination due to its zero net electric dipole moment; see Supplementary Note 2). The coupling between the resonator mode and the layer of molecules—resulting from the strong overlap between the near field of the resonator mode and the molecules—is probed by recording of the tip-scattered field with an asymmetric Fourier transform spectrometer, yielding both near-field amplitude and phase spectra, s 3 (ω)andφ 3 (ω), respectively, 2.5 +-++ - CBP 1400 14801440 0.9 0.8 1.0 1.2 1.1 1.3 L (Pm) CBP SiO 2 s  /s 3,Au SiO 2 CBP h-BN Tip E sca E inc 500 nm = 1437.4 cm -1 E CBP h-BN z 1.0 0.9 0.8 0.6 0.7 0.5 Ez/Ez,Au L (Pm) 1400 14801440 CBP 3 0 0 (cm -1 ) Z (cm -1 ) Z Z 0.5 1.0 1.5 2.0 1 2 CBP Z CBP Z y x y cd e b x y z a Fig. 1 | Tip-enhanced near-field probing of half molecule-covered h-BN nanoresonators. a Illustration of the experiment. bNumerical simulation of the electric field distribution at ω=1437.4cm −1around an h-BN nanorod of 1000 nm length, 250 nm width and 87 nm height, whose right half is covered by a 50 nm thick CBP layer. The vertical arrow indicates a point-dipole source mimicking the tip. cTopography image of h-BN nanorods (length L, 87 nm height and 250 nm width) that are half covered by a 50 nm thick CBP layer. dExperimental nano-FTIR amplitude spectra recorded at the positions marked in panel cby dots of the respective color. For better visualization, spectra for L= 1.0 and 1.3 μmarescaledby a factor of 0.7. eSimulated nano-FTIR amplitude spectra, as explained in Fig. 3.For the width and thickness of the h-BN rods we used the nominal experimental values. The lengths Lwere chosen such that the second-order FP mode tunes across the molecular vibrational resonance of the CBP molecules at ω= 1450 cm−1. We attribute the differences between Lin the experiment and the simulation to fabrication uncertainties such as trapezoid-like rod cross sections.9,10 d,eGreen curves show experimental and simulated nano-FTIR amplitude spectra of a 50 nm thick bare CBP layer on a 250 nm thick SiO 2 on Si substrate. Red lines are guides to the eye and mark peak positions. Gray-dashed lines indicate the CBP vibrational resonance at ω CBP = 1450 cm−1whose line width is 6.4 cm−1(ref. 11). All spectra are normalized to that obtained on a Au reference surface, and are offset. We note that in the experimental spectra (Fig. 1d) we observe several smaller peaks, e.g., in the frequency range between 1420 and 1440 cm−1. They can be attributed to higher-order PhP modes10. In the simulated spectra, where the near-field probe is described by a point-dipole source (Fig. 1e), these peaks are only partially seen. A much better reproduction of these peaks is obtained in numerical simulations, where the nearfield probe is described by a conical metallic tip (for further details and discussion see Supplementary Note 4). Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 2 where the index indicates that the detector signal was demodulated at 3Ωto suppress background signals (nano-FTIR spectroscopy; see Methods). A topography image of the set of h-BN nanorods of 250 nm width, 87 nm height and various lengths LisshowninFig.1c, where one can also see the thin, homogeneous CBP layer of 50 nm thickness covering the right half of all resonators. The experimental nano-FTIR amplitude spectra (Fig. 1d; phase spectra are shown in Supplementary Fig. 3) recorded on the left rod extremity (measurement position marked by dots in Fig. 1c) clearly show the resonance peaks of the first and second-order FP resonances (illustrated by the schematic above the diagram in Fig. 1dandverified by the experimental mode pattern shown in Fig. 2), which shift to higher frequencies when the nanoresonator length Lis reduced from 1300 to 800 nm. However, the peak of the second order mode does not cross the CBP vibrational resonance at ω CBP =1450cm −1. Tracing the peak positions (marked by red lines) indeed reveals anti-crossing behavior, indicating that the nanoresonator near field couples with the molecular vibrations of the CBP layer. Most important, the nanoresonator-molecule coupling (occurring on the right half of the h-BN nanorods) can be well probed when the tip is placed on the left nanorod extremity, that is several 100 nm away from the molecules, where a direct near-field interaction between tip and molecules can be neglected (note that significant near-field interaction between tip and sample occurs only for distances smaller than the tip apex radius, here about 25 nm). The experimental near-field spectra can be well reproduced by numerical simulations (Fig. 1e), where the tip is modeled as a point-dipole source, which are described in more detail in Fig. 3. Spatio-spectral near-field mapping of phononic nanoresonators Remote near-field spectroscopy can be applied for nanoscale spatial mapping of the resonator-molecule coupling via hyperspectral nanoimaging—even for dark modes that are not accessible by far-field spectroscopy. In the future, such a possibility could be applied, for example, to study advanced resonator structures where a variety of different resonator modes may coexist and couple with the molecular vibrations. In Fig. 2we apply this capability for direct experimental identification of the phononic resonator mode that couples with the molecular layer and to verify that the peak splitting is caused by the presence of the molecules. To that end, we recorded nano-FTIR amplitude spectra along the principal axis of a bare (Fig. 2a, b) and a half-covered (Fig. 2c, d) h-BN nanoresonator. The lengths Lof the nanoresonators were chosen such that their second-order FP resonance occurs at the molecular vibrational resonance frequency, ω CBP =1450cm −1. For each spectral peak we observe strong near-field oscillations along the principal nanoresonator axis. The number of oscillations increases steadily with increasing frequency, revealing a series of longitudinal FP modes. The three near-field maxima (in the center and at the two extremities of the nanoresonator) at 1450 cm−1 clearly reveal the second-order (dark) FP mode. More important, the second-order FP modal near-field pattern of the half-covered nanoresonator exhibits a small spectral dip at the molecular vibrational resonance of CBP at 1450 cm−1(dark stripe in Fig. 2d) all along the principal resonator axis, which is absent in the hyperspectral linescan of the bare nanoresonator (Fig. 2b). For a better quantitative comparison, we show in Fig. 2e the nano-FTIR amplitude spectra recorded at the left extremity of the bare and molecule-covered h-BN nanoresonators (positions marked by blue and red dots in Fig. 2a, c, respectively). We clearly see that the peak of the first-order FP mode (far away from the CBP resonance at ω CBP =1450cm −1)isnearlythe same for both nanoresonators, whereas the peak of the second-order FP mode of the molecule-covered h-BN nanoresonator exhibits a spectral dip occurring at ω CBP , as typical for the coupling between molecular vibrations and nanoresonator modes7,9,10,23. The spatiospectral observation presented in Fig. 2demonstrates that the nanoresonator-molecule coupling can in principle be probed at any location where the nanoresonator mode can be activated by the nearfield probe. On the other hand, Fig. 2d reveals that the near-field signal on the molecule-covered part of the nanoresonator is strongly reduced due to the increased tip-nanoresonator distance and that the spectral line shape is modified due to the direct near-field interaction between tip and molecules, highlighting the advantages of probing the molecule-free part of the resonator, that is, probing larger near-field signals and reduction of undesirable near-field interaction between tip and molecules. The modification of the line shape is more clearly seen in Fig. 2f, which compares spectra recorded on the molecule-free nanoresonator part (B and C) with spectrum D that was recorded on the molecule-covered part. Specifically, we see that the spectrum D exhibits an asymmetric line shape and that the dip is shifted to slightly higher frequencies. This asymmetric line shape can be explained as CBP a b B A c d Height (nm) 0 140 0 0.5 1.0 1.5 min max 1500 1400 1450 s 3 /s 3,Au Position (Pm) 0.5 1.0 1.5 SiO2SiO2CBP CD CBP e B A 1390 14701430 B C D 1430 14701450 0 0.2 0.4 f 0 0.5 1.0 Norm. s3 Norm. s3 (cm -1 ) Z (cm -1 ) Z (cm -1 ) Z Position (Pm) Z Z 0 Fig. 2 | Nano-FTIR line scans of bare and half molecule-covered h-BN nanoresonators. a Topography image of a bare h-BN nanorod of 1.25μm length, 110 nm height and 250 nm width. bNano-FTIR amplitude spectra recorded along the horizontal dashed black line in panel a.cTopography image of an h-BN nanorod of 1.1 μm length, 87 nm height and 250 nm width, which is half covered with a 50 nm thick CBP layer. dNano-FTIR amplitude spectra recorded along the horizontal dashed black line in panel c.eNano-FTIR amplitude spectra of the bare (blue) and molecule-covered (red) h-BN nanorod at positions marked by blue and red dots in panels aand c, respectively. Both spectra are normalized to the peak maximum at 1408 cm−1.fNano-FTIR amplitude spectra of the molecule-covered h-BN nanorod at positions marked by red, orange and black dots in panel c. The spectra are normalized to the peak maximum at 1445 cm−1. The vertical green dashed line marks the frequency of the molecular vibrational resonance of CBP. Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 3 superposition of the symmetric nano-FTIR spectrum of the moleculecovered nanoresonator (B) and the asymmetric nano-FTIR spectrum of a bare CBP layer (such as the green spectrum in Fig. 1d) caused by the direct near-field interaction between tip and molecules. The later exhibits a derivative-like spectral line shape, that is, a peak and a dip to the left and right of the molecular resonance, respectively, which is a typical characteristic of nano-FTIR amplitude spectra of molecular layers29. Theoretical description of remote near-field probing of VSC To understand how the hybrid nanoresonator-molecule modes manifest in the near-field spectra, we first discuss in Fig. 3numerical simulations where the tip is modeled as a point-dipole source17,30 (red arrow in Fig. 3a, d) located above the h-BN resonator. We obtain complex-valued near-field spectrum by evaluating the vertical (z-) component of the electric field below the dipole (evaluation position marked by a cross in Fig. 3a, d) as a function of frequency, E z (ω). Further simulation details are described in the Methods section. Importantly, the dipole moment of the source is kept constant and consequently is not modified by the fields of the nanoresonator. This simplified modeling of the tip allows for excluding a potential coupling of the tip with the nanoresonator and the molecules, and thus exclusively reveals the spectral near-field signature of the coupling between the nanoresonator and the molecular vibrations. Figure 3b shows the simulated near-field amplitude spectrum |E z (ω)| (open symbols), of a h-BN nanorod that is half covered with a layer of permittivity ε= 2.8 (corresponding to the permittivity of CBP without the molecular resonance at 1450 cm−1). Its length was chosen such that the second-order FP resonance is at ω=1450cm −1,which manifests as a single peak in the near-field amplitude spectrum. Repeating the simulation when the h-BN nanorod is half covered by a CBP layer (Fig. 3e), the nanoresonators’near-field amplitude peak splits into two peaks (open symbols), which is a direct consequence of the coupling between the nanoresonator mode and molecular vibrations. For quantifying the coupling between the h-BN nanoresonator and the molecular vibrations, we fitthesimulatednear-field spectra by a model of two coupled harmonic oscillators (Methods), which represent the nanoresonator mode and the molecular vibration with eigenfrequencies ω i and damping parameters γ i (the index idenotes either PhP or CBP). The two oscillators are coupled to each other through the coupling strength g. For a most reliable analysis, we performed complex-valued fitting of the simulated near-field spectra EzωðÞ=∣EzωðÞ∣eiφωðÞ , where |E z (ω)| and φ(ω) are the amplitude and phase spectra (Fig. 3f). Plotting the fits by red solid lines in Fig. 3e, f, we find an excellent agreement with the simulated spectra (red open symbols), which confirms the validity of the two coupled harmonic oscillator model to fitthesimulatednear-field spectra. Before analysing the experimental nano-FTIR spectra, we note that the simulated near-field spectrum of the h-BN nanoresonator in absence of molecular vibrations yields a nearly circular trajectory in thecomplexplane(Fig.3c). Interestingly, when the nanoresonator is half covered with CBP molecules, the topology of the trajectory of the complex-valued near-field spectra changes. We find that a small loop appears within the main circular trajectory (Fig. 3f), which reveals the coupling between the resonator mode and the molecular vibrations, i.e., that not only a mere superposition of two modes is observed, (Supplementary Note 5), similar to a recent observation in far-field ellipsometry of excitonic coupling in classical microresonators31.The topology change of complex-valued trajectories may become an interesting means for characterizing coupling phenomena. Quantitative analysis of experimental nano-FTIR spectra For analysing the experimental near-field spectra (Fig. 4a, b), we first assume that the tip is solely illuminating the nanoresonator (i.e., that it acts like a dipole source with a constant dipole moment and for that reason does not need to be modeled as another coupled oscillator). We thus apply the same model of two coupled harmonic oscillators as for the simulated near-field spectra to fit the experimental spectra. We further assume that the PhP mode is not excited by far-field illumination because it is a dark mode. Figure 4b shows the experimental nearfield spectra of the half-covered nanoresonators in the complex plane (black lines), as well as the excellent fitting results (orange lines). The corresponding amplitude and phase spectra are shown in Fig. 4aand Supplementary Fig. 3, respectively. The small loops in the complex plane clearly reveal the coupling between the nanoresonator mode and the molecular vibrations, as predicted by the simulations in Fig. 3f. We note that the size of the small loop slightly varies in the experimental spectra of the different nanoresonators, which we attribute to variations of the resonators´ quality factor due to fabrication uncertainties. From the fits we determined for each experimental near-field spectrum (i.e. for each nanoresonator) the coupling strength g,thePhP eigenfrequency ω PhP and damping γ PhP , the molecular vibrational eigenfrequency ω CBP and damping γ CBP (see Methods and Supplementary Table 2). From these parameters we calculated the eigenfrequencies of the new hybrid modes (see Methods) for each nanoresonator, ω ± , which are shown in Fig. 4c (blue symbols) as a function of the nanoresonators´ eigenfrequencies ω PhP . We observe a clear anti-crossing of the hybrid modes, which indicates strong coupling. To determine the coupling regime for each individual resonator, we mark in Fig. 4e the transition from weak to strong coupling— defined by the standard criterium g=|γ PhP +γ CBP |/4(ref.32)—by black dashed horizontal lines. Interestingly, we find that gis above the black dashed line for all individual resonators, indicating that all of them are strongly coupled with the molecular vibrations. Importantly, applying the same fitting procedure to the simulated complex-valued near-field spectra (where the tip is modeled as a pointdipole source; see Methods and Supplementary Note 3) yields hybrid eigenmodes (Fig. 4d) and coupling strengths (Fig. 4f) that match well CBP h-BN Dipole 1440 14601450 1440 14601450 0.60 0.90 0.6 1.0 0.3 0 0.50 0.5 1 0.6 0.9 0 b Ez /Ez,Au ef a Re[Ez /Ez,Au] 0.7 0.8 0.1 0.2 0.75 0.25 0.8 0.75 Im[Ez /Ez,Au ] Re[Ez /Ez,Au] c d Z CBP Z Z Z CBP (cm -1 ) Z (cm -1 ) Z H  = 2.8 h-BN Dipole Ez /Ez,Au Im[Ez /Ez,Au ] Fig. 3 | Theoretical description of near-field probing of VSC. a Simulation geometry, showing a point-dipole source above an h-BN nanorod of 785 nm length, 87 nm height and 250 nm width, which is half covered with a 50 nm thick layer of a permittivity ε=2.8.bAmplitude of the z-component of the electric field below the dipole at the position marked by a cross in panel a(normalized to that obtained on a Au reference surface, ∣Ez=Ez,Au∣), as function of frequency ω.cz-component of the electric field below the dipole source at the position marked by a cross in panel a, plotted in the complex plane. The arrow indicates increasing frequency ω. dSimulation geometry, showing a point-dipole source above the same h-BN nanorod as in panel abut half covered with CBP molecules. eSame as panel b,but the h-BN nanorod is half covered by a 50 nm thick CBP layer. fSame as panel c,but for a half CBP-covered h-BN nanorod. b,c,e,fOpen symbols show simulation results. Solid lines show fits obtained with b,ca single harmonic oscillator model and e,fa coupled harmonic oscillator model describing the coupling between the nanoresonator modes and the molecular vibrations. Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 4 the experimental results. Since the tip is not involved in the nanoresonator-molecule coupling in the simulations, this quantitative agreement (for further discussion see Supplementary Note 6) lets us further assume that the coupling between the PhP mode and the molecular vibrations in our experiment can be well described with a simple two coupled harmonic oscillator model, without the need of considering the tip as a third oscillator. Our analysis thus demonstrates the capability of remote near-field spectroscopy to probe the hybrid modes arising from the coupling between a single PhP nanoresonator and nanoscale amounts of organic molecules. We note that the damping γ PhP obtained from the simulated nearfield spectra is a factor of three lower than that obtained from the experimental spectra (see Supplementary Note 3). We explain this finding by PhP scattering and absorption at fabrication-induced irregularities and damage of the h-BN nanorod edges9,10, which is not considered in the simulation. Although the damping has only a minor influence on the coupling strength gand on the determination of the hybrid eigenmodes (see Eq. (10) in Methods), it is a crucial parameter that determines whether a coupled system is in the weak or in the strong coupling regime. Marking the transition from weak to strong coupling in Fig. 4f by black horizontal dashed lines (analogue to Fig. 4e), we find that it occurs at much lower gfor the simulation as compared to experiment, which shows that in simulations we are deeper inside the strong coupling regime than in the experiment. Influence of the oscillating tip and signal demodulation So far, we have not explicitly considered in the simulations that the tip is a long metallic cone, that the tip is oscillating, and that the detector signal is demodulated at higher harmonics of the tip oscillation frequency. To elucidate the influence and impact of these key features on the spectral mode positions and linewidths, we performed numerical simulations where the tip is modeled as a metal cone and the tip oscillation and signal demodulation are considered. The results are shownintheSupplementaryNote6and7andsummarizedasfollows. We first applied the coupled harmonic oscillator model (same as used in Fig. 3)tofit the simulated near-field spectra where the tip is modeled as a metal cone but without considering the tip oscillation and signal demodulation. Compared to the fitting parameters obtained from the unperturbed molecule-covered resonator spectra (obtained from the point-dipole source simulations), we find that the presence of the metallic tip yields negligible spectral shifts of the bare h-BN nanoresonator mode and the hybrid polariton modes (<3 cm−1), a negligible change of the sum of damping parameters, γ PhP +γ CBP ,andonlyaslight reduction of the coupling strength by a factor of about 1.26. In a second simulation, we additionally implemented tip oscillation and signal demodulation. Fitting of the simulated spectra yields nearly the same coupling strengths and mode positions as before. However, for the sum of the damping parameters, γ PhP +γ CBP , we obtain values that are reduced by a factor of about 1.56 as compared to the results obtained from the simulations where the tip is modelled as a point-dipole source. Since the presence of the tip, its oscillation and signal demodulation can yield significantly different values for gand γ PhP +γ CBP ,we re-evaluate the coupling regime obtained in the experiment according to a modified condition, 1.26 g/[1.56(γ PhP +γ CBP )] > 0.25, where gand γ PhP +γ CBP are the parameters (listed in Table 1) obtained by fitting the experimental near-field spectra. We find that two of the nanoresonators (length L=1.2µmandL=0.9µm) satisfy the strong coupling criteria, whereas the other nanoresonators are at the onset of strong coupling. The transition from weak to strong couplingaccording to this re-evaluation—is marked in Fig. 4e by red dashed horizontal lines. Since the modified condition for evaluating the coupling regime is more conservative than the original one (g/[γ PhP +γ CBP ]>0.25),the Simulations Experiment c Re[ 3] 1440 1450 1460 1440 1450 1460 1440 1450 1460 e 1440 1450 1460 1440 1450 1460 5 10 0 d g(cm-1) (cm -1 ) f V Z sc sc wc + Z - Z PhP Z CBP Z PhP (cm-1) Z + Z - Z PhP Z CBP Z (cm -1 ) Z Z CBP Z s 3 /s 3,Au ab CBP 1 2 0.9 0.8 1.0 1.1 L= 1.2 Pm 1420 1450 1480 0.2 0.3 0.1 0.25 0.15 -0.3 0.3 0 0.2 0.3 0.1 0 0.10-0.1 0.10 0-0.6 -0.3 000-0.1-0.1-0.1-0.2-0.2-0.2 0.10.10.1 wc PhP (cm-1) Z PhP (cm-1) Z PhP (cm-1) Z Im[ 3] V Fig. 4 | Quantitative analysis of nano-FTIR spectra by coupled harmonic oscillator fitting. a Black curves show the experimental nano-FTIR amplitude spectra of h-BN nanorods of length L, which are half covered with CBP molecules (same data as in Fig. 1d). Orange curves show fits obtained using the coupled harmonic oscillator model. For better visualization, the spectrum for L=1.0μmis scaled by a factor of 0.7. bBlack curves show the complex-valued experimental nano-FTIR spectra, σ3ωðÞ=s3ωðÞeiφ3ωðÞ , plotted in the complex plane. Orange curves show fits obtained using the coupled harmonic oscillator model. cEigenfrequencies ω ± of the hybrid modes, the nanoresonators´ bare eigenfrequencies ω PhP (gray squares), and bare molecular vibrational eigenfrequency ω CBP (green squares), all of them obtained by fitting of the complex-valued experimental nano-FTIR spectra. dSame as panel c,butobtainedbyfitting complex-valued simulated near-field spectra of h-BN nanoresonators of length L=1.2to0.8μm, which are half covered by CBP molecules (the corresponding amplitude spectra are shown in Fig. 1e; the complex-valued spectra are shown in Supplementary Fig. 3). The tip is modeled by a point-dipole source as in Fig. 3.e,fCoupling strength gobtained from the fitting of the experimental and simulated near-field spectra, respectively. Black and red dashed horizontal lines indicate the transition from weak (WC) to strong (SC) coupling defined by g=|γ PhP +γ CBP |/4ref.32 and g=|γ PhP +γ CBP | /3.23 (see discussion below), respectively. Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 5 transition occurs in Fig. 4e at larger gvalues. We note that the factors used to re-evaluate the coupling regime are specific for our study and may change for different tip oscillation amplitudes or signal demodulation orders. Discussion In summary, we demonstrated that infrared near-field spectroscopy can be applied for nanoscale spatial mapping of strong coupling between molecular vibrations and resonating dark PhP modes. To that end, we introduced the concept of remote near-field probing with a non-resonant tip, where the tip and molecules are separated such that the direct near-field interaction between tip and molecules is avoided. Such minimal-invasive probing offers the advantage that the hybrid nanoresonator-molecule modes and coupling strengths can be determined within the model of two coupled harmonic oscillators without considering the tip as a third oscillator, which significantly simplifies the coupling analysis and increases its robustness. On the other hand, numerical simulations of remote near-field spectroscopy including the metal tip, its oscillation and signal demodulation reveal that damping parameters may be underestimated when tip oscillation and signal demodulation are not taken into account for fitting the near-field spectra with a coupled harmonic oscillator model. Since near-field spectroscopy is available in the wide spectral range between visible and terahertz frequencies, we envision remote near-field probing studies of strong coupling of different plasmonic and phononic resonators (which apart from h-BN could be made, for example, from VO 5 or MoO 3 that cover various infrared frequency ranges) with molecular vibrations or excitons in various resonator geometries. The possibility to probe dark modes and spatially control the selective excitation and probing of coexisting modes could pave the way to explore strong coupling configurations that are not accessible by far-field spectroscopy. Methods Growth of monoisotopic 10B h-BN crystals The growth of 10B enriched h-BN crystals was done by the metal flux method33. Powders of 40 g Fe, 25 g Cr, 1 g 10B were mixed in an alumina crucible and placed in a single-zone furnace. After a dwell time of 24 h at 1550 °C, the furnace was cooled at a rate of 4 °C/h to 1400 °C, to precipitate the h-BN crystals on the metal surface. This was subsequently followed by quick quenching to room temperature. Fabrication of h-BN nanoresonators We first mechanically exfoliated 10B enriched h-BN crystals using blue Nitto tape (Nitto Denko Co., SPV 224 P). Then the flakes were exfoliated from the tape onto a transparent polydimethyl-siloxane (PDMS) stamp. Flakes of appropriate thickness and size were then transferred onto a Si/SiO 2 (250 nm) substrate using the deterministic dry transfer technique34. For patterning the h-BN nanoresonators we used highresolution electron beam lithography using a double-layer poly(- methyl methacrylate) (PMMA) resist (495 A4/950 A2). After development of the resist, a 3 nm-thick Cr layer was deposited onto the sample by e-beam evaporation, followed by thermal evaporation of 40 nm of Al. After lift-off in acetone, chemical etching of the h-BN flake was performed using a SF6/Ar 1:1 plasma mixture at 20 sccm flow, 100 mTorr pressure and 100 W power for 60 s (RIE OXFORD PLASMALAB 80 PLUS reactive ion etcher). For removing the Cr-Al mask from the h-BN structures, the sample was immersed in a chromium etchant (Sigma-Aldrich Co., 651826) for 20 min. After rinsing in deionized water, the sample was dried using a N 2 gun. Thermal evaporation of CBP 4,4′-bis(N-carbazolyl)-1,1′-biphenyl with sublimed quality (99.9%) (Sigma Aldrich) was thermally evaporated in an ultra-high vacuum evaporator chamber (base pressure <10−9mbar), at a rate of 0.1 nm s−1 using a Knudsen cell. Preparation of a diamond mask for structured CBP deposition A diamond mask was used for depositing the thermally evaporated CBP molecules only on one half of the h-BN nanoresonators. The mask was cut from a polycrystalline diamond film (500 nm thick) by focused Ga-ion beam (FIB) milling in a Helios 450S (FEI) dual beam system (milling conditions: 30 kV and 9.3 nA). Using an in situ platinum deposition, the diamond mask was attached to the tungsten tip of an Omniprobe micromanipulator and physically placed on top of the nanoresonators such that it covered half of each. The mask was fixed on the sample surface by in situ Pt-deposition. After CBP deposition by thermal evaporation (see above), the Omniprobe micromanipulator was approached to the mask. The mask was fixed to the Omniprobe by in situ Pt deposition. By FIB milling, the mask was detached from the sample and could be removed via the Omniprobe micromanipulator. Nano-FTIR spectroscopy We used a commercial scattering-type scanning near-field optical microscope (s-SNOM) setup comprising a nano-FTIR module (Neaspec/Attocube, Germany), in which the oscillating (at a frequency Ω≅ 270 kHz) Pt/Ir-coated AFM tip (Arrow-NCPt-50, Nanoworld, NanoWorld AG, Neuchâtel, Switzerland) was illuminated by p-polarized midIR broadband radiation generated by a supercontinuum laser (average power of about 0.5 mW; frequency range 1200–1700 cm−1). The detector signal was demodulated at a frequency 3 Ωfor effective background suppression. Interferograms were measured by recording the demodulated detector signal as a function of the position of the reference mirror, d,atafixed tip position. For apodization of the interferograms a Planck-taper window function with ε=0.2was applied. After zero-filling (4× padding) we Fourier transformed the interferograms to obtain complex-valued near-field point spectra, E s (ω). Each spectrum of Fig. 1d is an average of three spectra recorded in steps of 40 nm along the y-axis. We normalized the obtained point spectra to a reference spectrum recorded on gold, E s,Au (ω). The spectral resolution was 6.25 cm−1. Numerical simulations where the tip is modeled as a pointdipole source We performed the numerical simulations (Figs. 1e, 3,4d, f) using the Radio Frequency Module of COMSOL Multiphysics software. This module solves Maxwell’s equations in the frequency domain based on the Finite Element Method (FEM). The h-BN nanoresonators were modeled as a rectangular structure of w= 250 nm width (along the xdirection), d= 87 nm thickness (along the z-direction) and variable length L(along the y-direction). The resonators are on top of a 250 nm thick layer with the permittivity of SiO 2 , which is on top of a Si substrate. The CBP layer was modeled as a 50 nm thick layer covering half the length of the h-BN nanorod. The material permittivities are provided below. To ensure numerical convergence of the simulated nearfield spectra, the complete structure (point-dipole source, h-BN Table 1 | Parameters obtained by fitting the experimental complex-valued spectra shown in Fig. 4bviathecoupled harmonic oscillator model L(µm) ω PhP γ PhP ω CBP γ CBP gγ PhP +γ CBP 1:26g 1:56ðγPhP +γCBP Þ 1.2 1441.2 9.4 1450.8 6.0 5.4 15.4 0.28 1.1 1447.8 8.9 1452.3 6.0 4.3 14.9 0.23 1.0 1451.0 8.8 1451.4 7.0 4.6 15.8 0.23 0.9 1454.4 10.2 1451.8 6.6 6.0 16.8 0.29 0.8 1462.4 10.7 1452.0 7.0 4.7 17.7 0.21 Allparameters from the second to the seventh column are expressed in cm−1. Errors according to the standard deviation obtained by fitting the individual spectra are presented in Supplementary Table 2. Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 6 nanorod and CBP layer) was located in a homogeneous rectangular box (filled with air) of 8 × wwidth, 25 × ddepth and 4 × Llength. We use perfectly matched layers (PML) for the boundaries of the simulation box and free triangular elements for the nanorod mesh and free tetrahedralelementsforallotherstructures. To obtain the numerical results shown in Figs. 1e, 3,4d, f, we modeled the tip as a point-dipole source oriented perpendicularly to the substrate (along the z-direction). This model assumes that the elongated tip in the experiment is oriented perpendicularly (z-direction) to the h-BN nanorod and is illuminated by p-polarized light. The experimental signal detected in the far field (E s ) was approximated by the vertical component of the electric near field (E z ) at an evaluation point below the point-dipole source. In all the calculations, the pointdipole and the evaluation points were located at coordinates (x=0, y=−L/2 + 50 nm,z= 350 nm) and (x=0,y=−L/2+50nm,z =65nm), respectively. We set the origin (x=0,y=0,z= 0) at the middle of the top surface of the nanorod. The CBP layer covered the nanorod for y> 0. For comparison with the experimental near-field spectra, the simulated near-field spectra were normalized to E z obtained when the dipole source is located above a gold (i.e. reference) substrate. Permittivity of h-BN The permittivity of the isotopically (10B) enriched h-BN (in-plane εhBN,?, out-of-plane εhBN,k) is described by a Drude-Lorentz model35 εhBN,j=ε1,j1+ ω2 LO,jω2 TO,j ω2 TO,jω2iωΓj ! ,ð1Þ where jindicates the in-plane (⊥) or out-of-plane (||) component. In Eq. (1), ω TO,j and ω LO,j are the TO and LO phonon frequencies, Γjis the damping constant and ε1,jis the high-frequency permittivity. The values used for each constant are presented in Table 2. Permittivity of CBP We model the permittivity of CBP as the sum of a non-resonant background plus four oscillators which describe the molecular vibrations of CBP within the frequency range of interest. Thus, the permittivity of CBP is given by11 εCBP =ε1,CBP +X 4 k=1 ω2 p,k ω2 0,kω2iωΓCBP,k ,ð2Þ where ε1,CBP =2:8 is the high-frequency permittivity and ω p,k ,ω 0,k and Γ CBP,k are the strength, natural frequency, and damping constant of the k-th oscillator, respectively. The values used for each constant are presented in Table 3. Permittivity of Si, SiO 2 and Au For the numerical calculations shown in the main text, we considered the h-BN nanorod to be on top of a Si/SiO 2 substrate. For the permittivity of Si, we used a constant value ε si =12.TheSiO 2 permittivity was approximated by the following analytical function εSiO2=a0+a1λ0+a2λ2 0+a3λ3 0+a4λ4 0+a5λ5 0+a6 λ0a7  2+a2 8 +ia9 λ0a7  2+a2 10 , ð3Þ with λ 0 being the free space wavelength expressed in microns, a 0 =942.883,a 1 =−712.749 μm−1,a 2 =215.563μm−2,a 3 =−32.5483 μm−3, a 4 = 2.45449 μm−4,a 5 =−0.0740283 μm−5,a 6 =−0.1404 μm2, a 7 =8.078μm, a 8 = 0.1404 μm, a 9 =−0.1426 μm2and a 10 = 0.1426 μm. The permittivity of Au was taken from ref. 36. Coupled harmonic oscillator model To fit the data shown in Figs. 3and 4,wedescribedthecoupling between the Fabry-Perot phonon polariton resonance and the molecular vibrations by a classical model of two coupled harmonic oscillators. The equations of motion that describe the coupled system are determined by the following relations37–40 € xPhP tðÞ+γPhP _ xPhP tðÞ+ω2 PhPxPhP tðÞ2g_ xCBP tðÞ=FPhP tðÞ,ð4Þ € xCBP tðÞ+γCBP _ xCBP tðÞ+ω2 CBPxCBP tðÞ+2g_ xPhP tðÞ=FCBP tðÞ,ð5Þ where the dots denote time derivatives, x PhP (t) represents the PhP mode with resonance frequency ω PhP and damping constant γ PhP ,and x CBP (t) represents the molecular vibration of CBP with resonance frequency ω CBP and damping constant γ CBP .gis the coupling strength between the two resonators and F PhP (t), F CBP (t) represent the effective forces that drive each resonator. The effective forces are proportional to the near fields provided by the tip. In particular, we set F CBP (t)=0as in our experiment the near fields of the tip do not act directly on the CBP molecules. Through a time-to-frequency Fourier transform of Eqs. (4)and(5), we find the following steady-state solutions ^ xPhP ωðÞ=ω2 CBP ω2iγCBPω ω2 PhP ω2iγPhPω  ω2 CBP ω2iγCBPω  2gωðÞ 2 ^ FPhP ωðÞ, ð6Þ ^ xCBPðωÞ=i2gω ω2 PhP ω2iγPhPω  ω2 CBP ω2iγCBPω  2gωðÞ 2 ^ FPhP ωðÞ, ð7Þ with ^ xPhP ωðÞ=F½xPhP tðÞ,^ xCBP ωðÞ=F½xCBP tðÞ, ^ FPhP ωðÞ=F½FPhP tðÞ where Fis the time-to-frequency Fourier transform. Fitting of harmonic oscillator model to near-field spectra To fit the simulated spectra (Figs. 3and 4d), we used Eq. (6)plusa complex-valued offset term (x 0 +y 0 i), which accounts for the nonpolaritonic near-field interaction between the tip and the Table 2 | Parameters for calculating the in-plane and out-ofplane permittivity components of h-BN according to Eq. (1) In-plane (ε ⊥ )Out-of-plane(ε || ) ε ∞,j 32.8 ω TO,j 1395 cm−1785 cm−1 ω LO,j 1630 cm−1845 cm−1 Γ j 2cm −11cm −1 Table 3 | Parameters for calculating the permittivity of CBP according to Eq. (2) kω 0,k [cm−1]ω p,k [cm−1]Γ CBP,k [cm−1] 1 1450.0 128 6.4 21478.6 47 4.4 3 1500.1 91 9.4 41507.4 99 6.1 Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 7 nanoresonator: ^ xsim ωðÞ=^ xPhP ωðÞ+x0+y0i:ð8Þ To fit the experimental spectra (Fig. 4a, b), we multiply Eq. (6)bya complex-valued factor (eiϕ) that accounts for the drift of the interferometer between the measurement of the nanoresonator spectra and the reference measurement on a flat gold surface: ^ xexp ωðÞ=^ xPhP ωðÞeiϕ+x0+y0i:ð9Þ In the fitting procedure, we used as free parameters ω PhP ,γ PhP ,g, ^ FPhP ωðÞ,x 0 ,y 0 and ϕ. The parameters ω CBP and γ CBP are fitting parameters limited to the ranges (1450, 1452.5) cm−1and (6,7) cm−1, respectively. From the fit parameters we obtain the eigenfrequencies ω ± of the coupled system using the following relation ω±=1 2ωPhP +ωCBP  ±1 2Re ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4∣g∣2+Δωi 2Δγ  2 s 2 43 5,ð10Þ with Δω=ω PhP -ω CBP and Δγ=γ PhP -γ CBP .Equation(10) can be obtained by solving for the eigenvalues of Eqs. (4)and(5), together with the secular approximations ω2 PhP ω2 ±  =ðωPhP +ω±ÞðωPhP ω±Þ≈2ω±ðωPhP  ω±Þand ω2 CBP ω2 ±  =ðωCBP +ω±ÞðωCBP ω±Þ≈2ω±ðωCBP ω±Þ. Data availability Data that support the results of this work are available upon reasonable request from the corresponding author. References 1. Shalabney, A. et al. Coherent coupling of molecular resonators with a microcavity mode. Nat. Commun. 6,5981(2015). 2. George, J., Shalabney, A., Hutchison, J. A., Genet, C. & Ebbesen, T. W. Liquid-phase vibrational strong coupling. J. Phys. Chem. Lett. 6, 1027–1031 (2015). 3. Thomas, A. et al. Ground-state chemical reactivity under vibrational coupling to the vacuum electromagnetic field. Angew. Chem. - Int. Ed. 55, 11462–11466 (2016). 4. Thomas, A. et al. Tilting a ground state reactivity landscape by vibrational strong coupling. Science 363,615–619 (2019). 5. Vergauwe,R.M.A.etal.Modification of enzyme activity by vibrational strong coupling of water. Angew. Chem. - Int. Ed. 58, 15324–15328 (2019). 6. Garcia-Vidal, F. J., Ciuti, C. & Ebbesen, T. W. Manipulating matter by strong coupling to vacuum fields. Science 373, eabd0336 (2021). 7. O’Callahan, B. T., Hentschel, M., Raschke, M. B., El-Khoury, P. Z. & Lea, A. S. Ultrasensitive tipand antenna-enhanced infrared nanoscopy of protein complexes. J. Phys. Chem. C 123, 17505–17509 (2019). 8. Dayal, G., Morichika, I. & Ashihara, S. Vibrational strong coupling in subwavelength nanogap patch antenna at the single resonator level. J. Phys. Chem. Lett. 12,3171–3175 (2021). 9. Autore, M. et al. Boron nitride nanoresonators for phonon-enhanced molecular vibrational spectroscopy at the strong coupling limit. Light Sci. Appl. 7,17172–17178 (2018). 10. Autore, M. et al. Enhanced light–matter interaction in 10B monoisotopic boron nitride infrared nanoresonators. Adv. Opt. Mater. 9, 1–9 (2021). 11. Bylinkin, A. et al. Real-space observation of vibrational strong coupling between propagating phonon polaritons and organic molecules. Nat. Photonics https://doi.org/10.1038/s41566-020-007253(2021). 12. Wang, T., Li, P., Hauer, B., Chigrin, D. N. & Taubner, T. Optical properties of single infrared resonant circular microcavities for surface phonon polaritons. Nano Lett. https://doi.org/10.1021/ nl4020342 (2013). 13. Caldwell, J. D. et al. Low-loss, extreme subdiffraction photon confinement via silicon carbide localized surface phonon polariton resonators. Nano Lett.https://doi.org/10.1021/nl401590g (2013). 14. Caldwell, J. D. et al. Sub-diffractional volume-confined polaritons in the natural hyperbolic material hexagonal boron nitride. Nat. Commun. https://doi.org/10.1038/ncomms6221 (2014). 15. Caldwell, J. D. et al. Low-loss, infrared and terahertz nanophotonics using surface phonon polaritons. Nanophotonics 4,44–68 (2015). 16. Wang, T. et al. Phonon-polaritonic Bowtie nanoantennas: controlling infrared thermal radiation at the nanoscale. ACS Photonics 4, 1753–1760 (2017). 17. Alfaro-Mozaz, F. J. et al. Nanoimaging of resonating hyperbolic polaritons in linear boron nitride antennas. Nat. Commun.8, 15624 (2017). 18. Tamagnone, M. et al. Ultra-confined mid-infrared resonant phonon polaritons in van der Waals nanostructures. Sci. Adv.https://doi. org/10.1126/sciadv.aat7189 (2018). 19. Duan, J. et al. Active and Passive Tuning of Ultranarrow Resonances in Polaritonic Nanoantennas. Adv. Mater. e2104954 https://doi.org/ 10.1002/adma.202104954 (2021). 20. Yu, S.-J. et al. Ultrahigh-quality infrared polaritonic resonators based on bottom-up-synthesized van der Waals nanoribbons. ACS Nano 16,3027–3035 (2022). 21. Keilmann, F. & Hillenbrand, R. Near-fieldmicroscopybyelasticlight scattering from a tip. Philos.Trans.R.Soc.AMath.Phys.Eng.Sci. 362,787–805 (2004). 22. Chen, X. et al. Modern scattering-type scanning near-field optical microscopy for advanced material research. Adv. Mater. 31, 1–24 (2019). 23. Muller, E. A. et al. Nanoimaging and control of molecular vibrations through electromagnetically induced scattering reaching the strong coupling regime. ACS Photonics 5,3594–3600 (2018). 24. Metzger, B. et al. Purcell-enhanced spontaneous emission of molecular vibrations. Phys.Rev.Lett.https://doi.org/10.1103/ PhysRevLett.123.153001 (2019). 25. Huth, F. et al. Resonant antenna probes for tip-enhanced infrared near-field microscopy. Nano Lett. 13,1065–1072 (2013). 26. Dai, S. et al. Tunable phonon polaritons in atomically thin van der Waals crystals of boron nitride. Science https://doi.org/10.1126/ science.1246833 (2014). 27. Yoxall, E. et al. Direct observation of ultraslow hyperbolic polariton propagation with negative phase velocity. Nat. Photonics https:// doi.org/10.1038/nphoton.2015.166 (2015). 28. Giles, A. J. et al. Ultralow-loss polaritons in isotopically pure boron nitride. Nat. Mater. https://doi.org/10.1038/NMAT5047 (2018). 29. Mastel,S.,Govyadinov,A.A.,DeOliveira,T.V.A.G.,Amenabar,I.& Hillenbrand, R. Nanoscale-resolved chemical identification of thin organic films using infrared near-field spectroscopy and standard Fourier transform infrared references. Appl. Phys. Lett.106, 023113 (2015). 30. Nikitin, A. Y. et al. Real-space mapping of tailored sheet and edge plasmons in graphene nanoresonators. Nat. Photonics 10, 239–243 (2016). 31. Thomas,P.A.,Tan,W.J.,Fernandez,H.A.&Barnes,W.L.Anew signature for strong light-matter coupling using spectroscopic ellipsometry. Nano Lett. 20,6412–6419 (2020). 32. Törmä, P. & Barnes, W. L. Strong coupling between surface plasmon polaritons and emitters: a review. Rep. Prog. Phys. 78, 13901 (2014). 33. Liu, S. et al. Single crystal growth of millimeter-sized monoisotopic hexagonal boron nitride. Chem. Mater. 30, 6222–6225 (2018). 34. Castellanos-Gomez, A. et al. Deterministic transfer of twodimensional materials by all-dry viscoelastic stamping. 2D Mater. https://doi.org/10.1088/2053-1583/1/1/011002 (2014). Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 8 35. Li, P. et al. Infrared hyperbolic metasurface based on nanostructured van der Waals materials. Science https://doi.org/10.1126/ science.aaq1704 (2018). 36. Olmon, R. L. et al. Optical dielectric function of gold. Phys.Rev.BCondens. Matter Mater. Phys. 86,1–9(2012). 37. Novotny, L. Strong coupling, energy splitting, and level crossings: a classical perspective. Am. J. Phys. 78, 1199–1202 (2010). 38. Wu, X., Gray, S. K. & Pelton, M. Quantum-dot-induced transparency in a nanoscale plasmonic resonator. Opt. Express 18, 23633 (2010). 39. Li, T. E., Subotnik, J. E. & Nitzan, A. Cavity molecular dynamics simulations of liquid water under vibrational ultrastrong coupling. Proc.NatlAcad.Sci.USA117,18324–18331 (2020). 40. Yoo, D. et al. Ultrastrong plasmon–phonon coupling via epsilonnear-zero nanocavities. Nat. Photonics 15,125–130 (2021). Acknowledgements This work was supported by the MCIN/AEI/10.13039/501100011033 under the María de Maeztu Units of Excellence Program (CEX2020001038-M) and the Projects RTI2018-094830-B-100, PID2021123949OB-I00, PID2019-107432GB-I00 and PID2021-122511OB-I00, as well as by the Graphene Flagship (GrapheneCore3, No. 881603). J.L. and J.H.E. are grateful for support from the Office of Naval Research (Award No. N00014-20-1-2474), for the BN crystal growth. S.V. acknowledges financial support by the Comunidad de Madrid through the Atraccionde Talento program (grant no. 2020-T1/IND-20041). C.M.-E., R.E., and J.A. received funding from grant no. IT 1526-22 from the Basque Government for consolidated groups of the Basque University. Author contributions R.H., I.D., and C.M.-E conceived the study. h-BN crystals were grown by J.L., supervised by J.H.E. I.D. fabricated the h-BN nanoresonators and performed the experiments and data analysis. C.M.-E. performed the calculations and data analysis. E.N. and E.M. fabricated the diamond mask. F. Calavalle deposited the molecular layer, supervised by F. Casanova and L.E.H. S.C. and F.J.A.-M. contributed to simulations. A.B. and R.E. participated in the data analysis. S.V. supervised the fabrication of samples. R.H., I.D., and C.M.-E. wrote the manuscript with input from all authors. All authors contributed to scientific discussion and manuscript revisions. R.H. and J.A. supervised the work. Competing interests R.H. is co-founder of Neaspec GmbH, a company producing scatteringtype scanning near-field optical microscope systems, such as the one used in this study. The remaining authors declare no competing interests. Additional information Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-34393-4. Correspondence and requests for materials should be addressed to Rainer Hillenbrand. Peer review information Nature Communications thanks the other anonymous reviewer(s) for their contribution to the peer review of this work. Reprints and permissions information is available at http://www.nature.com/reprints Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. © The Author(s) 2022 Article https://doi.org/10.1038/s41467-022-34393-4 Nature Communications | (2022) 13:6850 9