Plasmonic nanocavity-enabled universal detection of layer-breathing vibrations in two-dimensional materials

TL;DR

This paper introduces a new method using plasmonic nanocavities to detect hidden layer-breathing vibrations in 2D materials, enabling better understanding of their interfacial interactions.

Contribution

A universal plasmon-enhanced Raman spectroscopy strategy using nanocavities to detect weak or inactive layer-breathing vibrations in 2D materials.

Findings

Plasmonic nanocavities enhance and detect layer-breathing modes in multilayer graphene and hBN.

The method modifies polarization selection rules of layer-breathing vibrations.

An electric-field-modulated model explains intensity profiles and interfacial polarizability modulation.

Abstract

Conventional Raman spectroscopy faces inherent limitations in detecting interlayer layer-breathing (LB) vibrations with inherently weak electron-phonon coupling or Raman inactivity in two-dimensional materials, hindering insights into interfacial coupling and stacking dynamics. Here, we demonstrate a universal plasmon-enhanced Raman spectroscopy strategy using gold or silver nanocavities to strongly enhance and detect LB modes in multilayer graphene, hBN, and their van der Waals heterostructures. Plasmonic nanocavities even modify the linear and circular polarization selection rules of the LB vibrations. By developing an electric-field-modulated interlayer bond polarizability model, we quantitatively explain the observed intensity profiles and reveal the synergistic roles of localized plasmonic field enhancement and interfacial polarizability modulation. This model successfully describes the behavior of plasmon-enhanced LB vibrations across different material systems and nanocavity geometries. This work not only overcomes traditional detection barriers but also provides a quantitative framework for probing interlayer interactions, offering a versatile platform for investigating hidden interfacial phonons and advancing the characterization of layered quantum materials.

Plasmonic nanocavities enable universal detection of layer-breathing modes in two-dimensional materials via plasmon-enhanced Raman spectroscopy, overcoming the limit of weak electron-phonon coupling or Raman inactivity, and delivering a quantitative framework for interlayer coupling analysis.

Figures

Full text

Introduction

The strategic design of metallic nanocavities enables localized surface plasmon resonance to confine the optical field into sub-wavelength volumes, achieving extraordinary local electromagnetic field enhancement by several orders of magnitude^1–4^. This provides unprecedented opportunities for plasmon-enhanced Raman spectroscopy (PERS) to uncover light-matter interactions at the nanoscale^5,6^. The technique leverages both electromagnetic and chemical enhancement mechanisms to boost the Raman signal of molecular vibration modes, achieving single-molecule detection sensitivity^7–13^. In addition, plasmonic nanocavities spatially reconfigure the incident optical field, modifying Raman selection rules and activating symmetry-forbidden modes in conventional Raman spectroscopy^14–16^. These unique characteristics render PERS an indispensable tool for surface and interface analysis at the molecular level, as well as for in situ characterization of nanomaterials^4,11,12,17^.

PERS further advances studies on the fundamental properties of two-dimensional materials (2DMs) by leveraging their atomically thin nature, which enables exceptionally efficient coupling with plasmonic fields. Among them, the remarkable enhancement and selective excitation of intralayer vibration modes in various 2DMs^15,16,18^ have enabled precise characterization of localized lattice dynamics^19^, edge-related properties^20^, and quantitative measurements of the plasmonic field enhancement in nanocavities^15,16^. A promising application of PERS lies in the investigation of low-frequency interlayer modes^21^ in 2DMs. The inherent Raman inactivity or inherently weak electron-phonon coupling (EPC) of some interlayer modes prevents their observation via conventional Raman spectroscopy^21,22^. Remarkably, the out-of-plane displacement fields of interlayer phonon modes exhibit a spatial extent comparable to subwavelength plasmonic near-fields, in stark contrast to the atomic-scale confinement of intralayer modes. This unique overlap suggests significant modifications in the PERS response, a currently underexplored phenomenon requiring systematic study.

In this study, we demonstrated significant enhancement of interlayer layer-breathing (LB) modes in multilayer graphene, hBN and their van der Waals heterostructures (vdWHs) coupled with plasmonic gold (AuNCs) or silver (AgNCs) nanocavities. Upon resonant excitation, the AuNCs and AgNCs generate strong localized plasmonic fields and induce polarizability modifications in 2DMs. These synergistic effects enhance the LB mode signals and modulate the intensity profile of the observed LB modes in AuNC- or AgNC-coupled multilayer graphene, hBN and related vdWHs. We developed an electric-field-modulated interlayer bond polarizability model (E-IBPM) that comprehensively accounts for both variations in interlayer bond polarizability and the inhomogeneity of the plasmonic field within 2DMs. This framework successfully explains the observed Raman intensity profile of LB modes in multilayer graphene, hBN and the related vdWHs coupled with AuNCs or AgNCs, a finding that bridges the gap between the localized plasmonic enhancement and the Raman intensity of interlayer modes. The proposed methodology promises to be applicable to other elusive quasiparticles, such as interlayer excitons and certain plasmonic resonances in vdWHs.

Results

Plasmon-enhanced Raman modes in AuNCs/NLG

In conventional Raman scattering of phonons, the Raman intensity is critically governed by the EPC strength. This renders interlayer phonons in 2DMs (e.g., the LB modes in multilayer graphene) with weak EPC strength particularly challenging to detect. PERS can overcome this longstanding limitation, as it can amplify weak Raman signals by simultaneously enhancing the electric field of incident/scattered light and modulating the polarizability of the materials. Motivated by these advantages, we employed PERS to characterize the LB modes of multilayer graphene.

Figure 1a illustrates the PERS schematic for detecting LB modes in N-layer graphene (NLG) through coupling with AuNCs. Under the plasmon resonance condition, AuNCs generate a strongly confined plasmonic field localized at the nanoscale around their surfaces to enhance the LB modes in NLG with inherently weak electron-phonon coupling or Raman inactivity. In our experiments, we fabricated AuNCs/NLG structures by depositing an 8 nm Au film on NLG. Notably, the N of the NLG flake was pre-determined using Raman spectroscopy^23^ prior to Au deposition. The Au nanofilm comprises nanoislands which are varied in size and irregular in shape, separated by numerous nanogaps, as elucidated by the scanning electron microscopy (SEM) and transmission electron microscopy (TEM) images (Fig. 1b), which have been proven to serve as AuNCs with good plasmonic activity^24,25^. Indeed, the dark-field scattering spectroscopy of AuNCs/18LG (Fig. 1c) suggests a plasmon resonance peak at around 633 nm. When performing low-frequency Raman measurements on AuNCs/18LG under both plasmon on-resonance (λL = 633 nm) and off-resonance (λL = 532 nm) excitation conditions, we found significant enhancement of Raman modes under on-resonance excitation compared to the off-resonance case (Fig. 1d), unambiguously confirming the critical role of plasmon resonance in amplifying the Raman signals.Fig. 1.Plasmon-enhanced Raman modes in AuNCs/NLG. a Schematic illustration of the AuNCs/NLG and the PERS of LB modes. b Optical (left panel, scale bar is 20 μm), SEM (middle panel, scale bar is 40 nm), and TEM (right panel, scale bar is 25 nm) images of AuNCs/NLG. c Dark-field scattering spectra of AuNCs/18LG and 18LG. d Low-frequency Raman spectra of AuNCs/18LG under on-resonance (λL = 633 nm) and off-resonance (λL = 532 nm) conditions. e Raman spectra of AuNCs/NLG and the corresponding pristine NLG with λL of 633 nm (N = 1, 2, 6, and 18). The “*" and “+” symbols indicate the LB modes and new low-frequency modes, respectively

To confirm the robustness of the enhanced Raman modes in AuNCs/NLG, we further acquired the Raman spectra of AuNCs/NLG and the corresponding pristine NLG (N = 1, 2, 6, and 18) excited by λL = 633 nm, as illustrated in Fig. 1e. The G mode (~1580 cm^−1^) of AuNCs/NLG exhibits significant enhancement compared with that of NLG when N = 1, 2, and exhibits obvious strain-induced mode splitting when N = 1. As N increases, the enhancement effect weakens and the mode splitting vanishes, which is consistent with previous results and can be attributed to the weaker interaction between AuNCs and NLG for large N^24^. The lack of a D mode (~1360 cm^−1^) in the Raman spectra (Fig. 1e) indicates minimal defects^26^ in the AuNCs/NLG fabrication process. This validates the use of the pre-determined N for constructing physical models to reliably interpret the measured plasmon-enhanced Raman spectra in AuNCs/NLG.

Beyond the high-frequency intralayer modes, the distinct low-frequency interlayer modes of NLG also warranted particular attention^21^. Theoretically, NLG has N-1 pairs of doubly degenerate C modes and N-1 LB modes, which can be denoted by C_N,N−j_ and LB_N,N−j_ (j = 1, 2,⋯, N − 1) modes, respectively. The C_N,1_ (LB_N,1_) mode represents the C (LB) mode with the highest frequency. After subtracting the background signals from AuNCs or substrates, clear low-frequency Raman features emerge in both AuNCs/NLG and NLG (see Supplementary Note 1 and Fig. S1). Figure 1e demonstrates that all LB modes remain absent in NLG with N = 2, 6, and 18 due to their exceptionally weak EPC, while weak C_N,1_ modes can be observed. Nonetheless, the C_N,1_ modes are not evident in AuNCs/NLG. Instead, a series of layer-number-dependent Raman features emerged below 130 cm^−1^, as marked by asterisks and crosses in Fig. 1e. The absence of C_N,1_ modes in AuNCs/NLG can be attributed to the strong spectral overlap between these new modes and the inherently weak C_N,1_ modes, combined with the significantly higher intensity of the new modes. Notably, the frequencies of certain new modes in AuNCs/NLG surpass those of the C_N,1_ modes, i.e., the highest-frequency C modes in pristine NLG^21^. Based on this evidence, these new modes are preliminarily assigned to plasmon-enhanced LB modes of AuNCs/NLG. Raman mapping of a typical AuNCs/16LG sample in Fig.S2 further confirmed the scalability and reproducibility of the PERS-based method for enhancing LB modes in NLG. Notably, there are subtle differences between the plasmon-enhanced LB modes of AuNCs/NLG and the expected LB modes in pristine NLG. For instance, a new low-frequency mode is observable in AuNCs/1LG, whereas no LB mode is expected in pristine 1LG. These distinctions will be further elaborated in subsequent sections.

Plasmon-enhanced LB modes in AuNCs/NLG

To elucidate the physical origins of plasmon-enhanced LB modes in AuNCs/NLG, systematic low-frequency Raman spectroscopic investigations (up to 140 cm^−1^) were conducted on AuNCs/NLG systems (1≤ N ≤ 30) using λL = 633 nm excitation. The background-subtracted spectra are depicted in Fig. 2a, while Fig. 2b summarizes the frequencies of the observed LB modes, revealing a pronounced N-dependence.Fig. 2.Plasmon-enhanced LB modes in AuNCs/NLG. a Low-frequency Raman spectra (up to 140 cm^−1^) of AuNCs/NLG, where N varies from 1 to 30. The IMs are marked by “+". b Extracted LB mode frequencies from (a) as a function of N (circles), along with the fitting results (lines) by the LCM. The inset illustrates the LCM for AuNCs/NLG and the corresponding LB force constants. c Linearly polarized (including VV and HV configurations) and d circularly polarized (including σ + σ + and σ + σ − configurations) Raman spectra of AuNCs/6LG and t(3 + 3)LG (θ = 10.5^∘^). The C and LB modes are labeled. All spectra were obtained with λL = 633 nm

To better understand these LB mode frequencies, the linear chain model (LCM)^21^ was employed for AuNCs/NLG, which treats N-layer graphene as a system of rigid masses connected by springs^21,27^. The rigid mass represents the mass per unit area of each graphene layer, and the force constants of springs correspond to the interlayer coupling strength between adjacent graphene layers. However, we observed N LB modes in AuNCs/NLG (N = 1, 2, 3), deviating from the established theoretical LCM framework that predicts (N − 1) LB modes in pristine NLG (N > 1). We note that the robust interfacial coupling at both the graphene/substrate and AuNCs/graphene interfaces would perturb the out-of-plane acoustic phonon, causing it to emerge at frequencies lower than the expected LB modes. This perturbed mode, referred to as the interface mode (IM) of the hybrid structure, has been observed in various 2DMs like CVD-grown Bi_2_Te_3_ and WS_2_^28,29^. A detailed comparison of frequencies and interlayer displacements for the out-of-plane acoustic modes and LB modes in NLG, along with the LB modes and IM in AuNCs/NLG, is provided in Fig.S3. Thus, we modified the LCM to incorporate the interfacial LB force constants (i.e., kGr/Sub at the graphene/substrate interface and kAu/Gr at the AuNCs/graphene interface) at the two interfaces, in addition to the nearest-neighbor interlayer LB force constant in NLG ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\rm{Gr}}}^{{\rm{1st}}}=10.7\times 1{0}^{19}$$\end{document}  N m^−3^) and the next-nearest-neighbor LB force constant ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\rm{Gr}}}^{{\rm{2nd}}}=0.9\times 1{0}^{19}$$\end{document}  N m^−3^)^30^. This extended model successfully describes both IM and LB mode frequencies. By fitting the extracted LB mode frequencies in Fig. 2b with the modified LCM, we obtained \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{{\rm{Gr/Sub}}}=0.2{k}_{{\rm{Gr}}}^{{\rm{1st}}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{{\rm{Au/Gr}}}=0.3{k}_{{\rm{Gr}}}^{{\rm{1st}}}$$\end{document} (See Supplementary Note 2 in SI for the full description of LCM). The excellent agreement between theoretical and experimental results in Fig. 2b underscores the essential role of interfacial coupling in properly understanding the LB mode frequency of AuNCs/NLG probed by PERS.

PERS follows fundamentally different selection rules compared to conventional Raman scattering, providing a clear spectroscopic feature to distinguish these two mechanisms^12,31^. This distinction is exemplified by comparing the polarization-dependent responses of twisted (folded) double trilayer graphene (t(3 + 3)LG) and AuNCs/6LG systems, as depicted in Fig. 2c and 2d. In t(3 + 3)LG with a twist angle of 10.5^∘^, the C_3,1_ and C_3,2_ modes of the 3LG constituent and the LB_6,1_ and LB_6,3_ modes of 6LG are observed. This distinct difference between the C and LB modes stems from the contrasting interfacial coupling in twisted multilayer graphene (tMLG): while interfacial shear coupling is very weak, the interfacial LB coupling strength approaches that of interlayer coupling in graphite and AB-stacked multilayer graphenes^30^. Consequently, the C modes’ layer displacements are localized within the corresponding constituents, whereas the LB modes involve collective vibrations across all stacked layers in twisted (folded) multilayer graphenes^30^. This behavior is a general characteristic of C and LB modes in vdWHs^32–35^, providing a diagnostic criterion to distinguish them in experimental observations once the thicknesses of all the constituents are known. The two modes of t(3 + 3)LG exhibit strong enhancement under 633 nm excitation, in resonance with the twist-angle-dependent van Hove singularity (VHS) transition^30,36^. However, they display markedly different polarized responses^30^: the C mode exhibits non-zero intensity in linearly parallel-polarized (VV), linearly crossed (HV), and circularly cross-polarized (σ + σ −) configurations, while the LB modes show remarkable intensity in VV and circular co-polarized (σ + σ +) polarizations. These observations perfectly align with conventional Raman selection rules^30,37^. In striking contrast, the LB mode intensity in AuNCs/6LG remains non-zero across all four polarization configurations while showing negligible polarization dependence. This suggests that the AuNCs do not just passively enhance the signal by plasmonic field enhancement, but actively modify the Raman tensor properties associated with the bond polarizability of 6LG. These observations unambiguously demonstrate that AuNCs substantially altered the interlayer bond polarizability properties and the spatial distribution of the incident electric field in the Raman spectroscopy of adjacent multilayer graphene.

Electric-field-modulated interlayer bond polarizability model

The intensity analysis of Raman modes in 2DMs provides deeper insights into their underlying scattering mechanisms. It is noteworthy that quantitatively comparing the absolute Raman intensity of a specific mode across different samples is inherently challenging, as it is influenced by a multitude of factors, including plasmonic field distribution, optical interference effects, and experimental conditions. Therefore, the relative intensity between the LB modes within the same spectrum often provides a more robust and physically insightful metric. For the interlayer modes in pristine 2DMs, the interlayer bond polarizability model (IBPM)^38^ well explains their relative intensity by computing the Raman tensors from polarizability changes induced by interlayer bonds between adjacent layers. Within the electric dipole approximation, the IBPM neglects spatial variations of the incident electric field, linking the Raman intensity of interlayer modes to the derivative of the polarizability tensor with respect to layer displacements. In contrast, the PERS of LB modes in AuNCs/2DM systems violates the electric dipole approximation due to a strongly confined plasmonic field with nanometer-scale spatial localization at the AuNCs/2DM interface. In addition, the interfacial interaction between 2DMs and AuNCs^6,11,12^ can modify the polarizability of the 2DMs. Thus, a new theoretical framework incorporating both local field enhancement and polarizability modulation is required to accurately model the plasmon-enhanced intensity profile of the LB modes in AuNCs/2DM systems.

In PERS, the Raman response of a plasmon-enhanced j-th mode can be characterized by the induced Raman dipole moment p, closely linked to its Raman polarizability tensor Δαj and the localized plasmonic field (ELoc) generated by AuNCs, i.e., p = Δαj ⋅ ELoc. Unlike free-space dipole emission, the radiation from this dipole near the metal nanostructure is strongly modified by the plasmonic effects. Both its emission properties and the local excitation field are enhanced through analogous mechanisms. The far-field detectable Raman signal is thus determined by this modified radiation. The radiation process is quantitatively described by the Green’s function G, which relates the self-reaction radiation field (ER) to the Raman dipole moment, ER = G ⋅ p = G ⋅ Δαj ⋅ ELoc^31^. Thus, the PERS intensity of the j-th Raman mode in 2DMs can be expressed as^31,38^

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{j}\propto \frac{{n}_{j}+1}{{\omega }_{j}}{\left|{\bf{G}}\cdot \Delta {{\boldsymbol{\alpha }}}_{j}\cdot {{\bf{E}}}_{{\rm{Loc}}}\right|}^{2}$$\end{document}

where ωj is the frequency of the j-th mode. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${n}_{j}=1/[\exp (\hslash {\omega }_{j}/{k}_{b}T)-1]$$\end{document} is the Bose-Einstein factor. G is always described by a radiation enhancement factor, depending on several factors, including the substrate geometry and optical properties, the dipole position, orientation, and its emission frequency. Estimating G is an intractable task. For practical calculations in parallel-polarized backscattering configurations, the nondiagonal elements in G can be ignored and the diagonal element in G can be estimated as G ≈∣ELoc∣/∣E0∣, assuming that the amplitude and polarization of the local field remain similar at both the laser and Raman frequencies.

According to the empirical bond polarizability model, the Raman polarizability change Δαj can be approximated as a summation of polarizability variations induced by individual bonds^39,40^ (see Supplementary Note 3). For interlayer modes in 2DMs, only interlayer bonds contribute to the overall polarizability changes^38^. Thus, for the LB_N,N−j_ mode, Δαj is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {{\boldsymbol{\alpha }}}_{j}={\sum }_{l}{{\boldsymbol{\alpha }}}_{l}^{{\prime} }\Delta {z}_{j,l}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\boldsymbol{\alpha }}}_{l}^{{\prime} }={\left(\partial {{\boldsymbol{\alpha }}}_{j}/\partial {z}_{j,l}\right)}_{0}$$\end{document} denotes the polarizability derivative associated with the l-th layer, and Δzj,l is the layer displacement of the l-th layer in LB_N,N−j_ modes, obtained from LCM calculations. To account for plasmonic enhancement, the localized plasmonic field at each layer (ELoc,l) and the corresponding Green’s function (Gl) should be considered. Consequently, the PERS intensity of a given LB mode can be expressed as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{j}\propto \frac{{n}_{j}+1}{{\omega }_{j}}{\left|\mathop{\sum }\limits_{l}{{\bf{G}}}_{l}\cdot ({{\boldsymbol{\alpha }}}_{l}^{{\prime} }\Delta {z}_{j,l})\cdot {{\bf{E}}}_{{\rm{Loc}},l}\right|}^{2}$$\end{document}

This expression highlights that the Raman dipoles at each layer contribute to the LB mode intensity, similar to the case of molecules in an inhomogeneous local field^41^. The above framework constitutes the core idea of IBPM^34,38,42,43^, and this plasmonic field-adapted version is termed the electric-field-modulated IBPM (E-IBPM).

Although the out-of-plane LB modes can, in principle, be more efficiently enhanced by the z-component of plasmonic field (ELoc,l,z)^18,44^, our experimental configuration employs a backscattering geometry \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z(xx)\bar{z}$$\end{document} , which selectively probes the x**x component of Δαj (Δαj,x**x) and couples to the x component of the local field (ELoc,l,x) in Eq.(2). While the high numerical aperture (NA = 0.9) of the objective may permit partial collection of signals enhanced by ELoc,l,z, the observed intensity profiles are dominated by enhancements arising from ELoc,l,x. Moreover, the intensity profiles of the ELoc,l,z-enhanced signals exhibit spatial distributions that closely resemble those enhanced by ELoc,l,x (see Supplementary Note 4). As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{l,xx}^{{\prime} }$$\end{document} is only associated with the interlayer bond on both sides of the l-th layer, the interlayer bond parameter ηl,x related to the x component of the interlayer bond is introduced to express \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{l,xx}^{{\prime} }$$\end{document} (see Supplementary Note 3 for detailed derivation). Thus, Δαj,x**x can be expressed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta {\alpha }_{j,xx}=\mathop{\sum }\limits_{l}{\alpha }_{l,xx}^{{\prime} }\Delta {z}_{j,l}=\mathop{\sum }\limits_{l}({\eta }_{l,x}-{\eta }_{l+1,x})\Delta {z}_{j,l}$$\end{document}

Since ∣Gl∣≈∣ELoc,l,x∣/∣E0∣, the intensity of LB modes and IMs is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{j}\propto \frac{{n}_{j}+1}{{\omega }_{j}}{\left|\mathop{\sum }\limits_{l}({\eta }_{l,x}-{\eta }_{l+1,x})\Delta {z}_{j,l}{E}_{{\rm{Loc}},l,x}^{2}\right|}^{2}$$\end{document}

The relative intensity of the LB modes is dictated by the relative magnitude of the plasmonic field at each layer (i.e., the decay of the plasmonic field across layers) while the absolute Raman intensity scales with the field’s absolute strength. This core physical insight remains perfectly valid regardless of the plasmonic cavity shapes and densities. Once the local field generated by the plasmonic cavities is known, the relative Raman intensity profile of plasmon-enhanced LB modes can be obtained. The local field generated by the AuNCs with λL = 633 nm was simulated by the finite-difference time-domain (FDTD) method. Figure 3a presents the simulated distribution of ∣ELoc,l,x∣/∣E0∣ corresponding to the arrangement of AuNCs in Fig. 1b. Figure 3b demonstrates a rapid decay of the maximum ∣ELoc,l,x∣/∣E0∣ with increasing layer index l.Fig. 3.Electric-field-modulated interlayer bond polarizability model for plasmon-enhanced LB modes in AuNCs/NLG. a FDTD-simulated distribution of normalized x-component of local field (∣ELoc,l,x∣/∣E0∣) in different graphene layers near AuNCs. b Simulated maximum∣ELoc,l,x∣/∣E0∣at graphene layers near AuNCs and a schematic of AuNCs/NLG with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{l,xx}^{{\prime} }$$\end{document} . c Measured PERS spectra with λL= 633 nm and d calculated IM and LB modes of AuNCs/NLG (N = 5, 10, 15, 20) by E-IBPM. The IM and LB modes are labeled. The lineshapes of all LB modes were modeled as Lorentzian functions, with FWHM of 10, 7, 4, and 3 cm^−1^ for N = 5, 10, 15, and 20, respectively. e SEM images and f the corresponding Raman spectra of as-prepared (As-prep.) AuNCs/8LG, AuNCs/8LG annealed (Ann.) at 600 °C, as-prepared AgNCs/8LG and AgNCs/8LG annealed at 400 °C. Scale bars in (e): 250 nm. The IMs and LB modes are marked with dashed lines in (f)

For AuNCs/NLG, only the 1st layer (adjacent to the AuNCs) and the N-th layer (adjacent to the substrate) that show asymmetric interlayer bonds exhibit a non-zero polarizability derivative^38^ (Fig. 3b), i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{1,xx}^{{\prime} }={\eta }_{{\rm{Au/Gr}}}-{\eta }_{{\rm{Gr}}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{N,xx}^{{\prime} }={\eta }_{{\rm{Gr}}}$$\end{document} , where ηAu/Gr and ηGr are interlayer bond parameters at AuNCs/graphene interface and within NLG, respectively. Interior layers (l = 2 to N − 1) show \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{l,xx}^{{\prime} }={\eta }_{{\rm{Gr}}}-{\eta }_{{\rm{Gr}}}=0$$\end{document} . Moreover, the strong interaction between AuNCs and NLG results in ηAu/Gr ≫ ηGr under the plasmon resonance condition, similar to the case in which the inherent interlayer bond parameters of pristine hBN layers is considerably smaller than the interfacial bond parameter in hBN/WS_2_ vdWHs under exciton resonance^34^. Therefore, the Eq.(2) for the LB_N,N−j_ mode and IM in AuNCs/NLG can be simplified as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{j}\propto \frac{{n}_{j}+1}{{\omega }_{j}}{\left|\Delta {z}_{j,1}\right|}^{2}$$\end{document}

Using the interlayer displacements from LCM calculations, we simulated the intensities of the LB_N,N−j_ modes and IM in AuNCs/NLG using ηAu/Gr = 10ηGr. Figure 3c, d compare the experimental low-frequency spectra with the simulated PERS spectra (generated via the E-IBPM) for the LB modes in AuNCs/NLG (N = 5, 10, 15, 20), in which the LB modes were modeled with Lorentzian lineshapes, with full width at half maximum (FWHM) values of 4, 3, 2, and 1 cm^−1^ for N = 5, 10, 15, and 20, respectively. Notably, the different FWHMs for NLG with varied N were determined according to the corresponding Raman spectra from experiments, as widely adopted in previous studies^38,43,45^. The close match between the simulated and experimental LB mode intensities confirms that the E-IBPM accurately describes their underlying physics. The applicability of E-IBPM to NLG coupled with AuNCs with different shapes or densities is further confirmed, as indicated in Fig. 3e, f and Supplementary Note 5. We also prepared plasmonic cavities with different shapes and densities by annealing AuNCs and AgNCs at different temperatures. The AuNCs annealed at 600 °C exhibit nanoislands separated by larger nanogaps than those of the as-prepared AuNCs. Similar to the as-prepared AuNCs, the as-prepared AgNCs comprised nanoislands with varying sizes and irregular shapes. In contrast, the AgNCs annealed at 400 °C show well-separated nanohemispheres (Figs. 3e and S4). The PERS spectra of 8LG coupled to these different plasmonic nanocavities exhibited pronounced plasmon-enhanced LB modes, where the peak positions show a slight difference due to the different interfacial coupling constants at AuNCs/8LG ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\mathrm{Au/Gr}}}={0.3}\,{k}_{\,{\mathrm{Gr}}}^{{1}\mathrm{st}}$$\end{document} ) and AgNCs/8LG ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\mathrm{Ag/Gr}}}={0.1}\,{k}_{\,{\mathrm{Gr}}}^{{1}\mathrm{st}}$$\end{document} ). This difference also leads to distinct relative Raman intensity profiles between AuNCs/8LG and AgNCs/8LG, as discussed in detail in Supplementary Note 5. This indicates that the E-IBPM effectively captures the essential physics occurring within the active plasmonic regions, providing a versatile and efficient framework for understanding the fundamental PERS mechanism in our system.

Plasmon-enhanced LB modes in various systems coupled to AuNCs or AgNCs

The generality of the plasmon-enhanced LB mode detection strategy has been further demonstrated in other 2DMs such as hBN, where the extremely weak EPC typically hinders the observation of LB modes^22,34^. While the Raman spectrum of pristine 16L-hBN shows only a faint C_16,1_ mode, AuNCs/16L-hBN under resonant excitation (λL = 633 nm) exhibits pronounced plasmon-enhanced LB modes after background subtraction (Fig. 4a). The assignment of these peaks to LB modes is confirmed by the agreement between experimental frequencies and calculated frequencies using the modified LCM (see Supplementary Note 6). The interlayer LB force constant of multilayer hBN is khBN = 9.9 × 10^19^ N m^−3^^27^, and the derived interfacial LB force constants are khBN/Sub = 2khBN and kAu/hBN = 0.1khBN. The significantly different interfacial LB coupling constants compared with those in NLG can be ascribed to the different interlayer distance between 2DMs and the underlying substrate (or Au nanofilm) and the mean charge densities at the interfaces^30^.Fig. 4.Plasmon-enhanced LB modes in hBN, NLG, and tMLG coupled to AuNCs or AgNCs. Experimental (Exp.) and calculated (E-IBPM or IBPM) Raman spectra of the LB modes for a AuNCs/16L-hBN, b1 AuNCs/t(3 + 3)LG (twist angle θ = 10. 5^∘^) and the corresponding pristine counterparts under 633 nm excitation. b2 The schematic diagram for interlayer bond parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{l,xx}^{{\prime} }$$\end{document} of AuNCs/t(3 + 3)LG. c Experimental (Exp.) and calculated (E-IBPM or IBPM) Raman spectra of the LB modes for AgNCs/30LG and AgNCs/50L-hBN under 532 nm excitation. Asterisks (*) indicate residual C modes

Building on the E-IBPM framework established for AuNCs/NLG, we extended the analysis to hBN by considering two critical interfacial parameters, i.e., ηAu/hBN (at the AuNCs/hBN interface) and ηhBN (within hBN), with ηAu/hBN = 10ηhBN. Non-zero polarizability derivatives emerge exclusively in the boundary layers, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{1,xx}^{{\prime} }={\eta }_{{\rm{Au/hBN}}}-{\eta }_{{\rm{hBN}}}$$\end{document} for the top layer and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{N,xx}^{{\prime} }={\eta }_{{\rm{hBN}}}$$\end{document} for the bottom layer. The simulated LB mode intensities of AuNCs/NL-hBN show excellent agreement with experimental data, where LB modes are modeled with Lorentzian lineshapes featuring a FWHM of 3 cm^−1^ (Fig. 4a).

For a complex case, the heterointerface induces non-zero polarizability derivatives for the layers adjacent to the heterointerface, leading to an additional contribution to LB mode intensities. This effect is clearly demonstrated in t(3 + 3)LG (Fig. 4b), where the twist-induced electronic VHSs in the electronic band structure generate non-zero interior polarizability derivatives \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{l,xx}^{{\prime} }$$\end{document} ^35,43^. Figure 4b1 compares the Raman spectra of t(3 + 3)LG and AuNCs/t(3 + 3)LG (θ = 10. 5^∘^) in σ^+^σ^+^ configuration measured at λL = 633 nm, where the twist angle of t(3 + 3)LG was identified by the emergent R and R’ modes (Fig. S5) in the high-frequency G-band region^46–48^. In pristine t(3 + 3)LG, only two high-frequency LB modes (LB_6,1_ and LB_6,3_) are observed. In contrast, in AuNCs/t(3 + 3)LG, both low- and high-frequency LB modes as well as IM are observed. Remarkably, the frequencies of the observed LB modes can be accurately described using the same LB force constants \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\rm{Gr}}}^{{\rm{1st}}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\rm{Gr}}}^{{\rm{2nd}}}$$\end{document} and kAu/Gr derived from AuNCs/6LG due to the negligible impact of twist interfaces on interlayer LB couplings^30^. The interaction between the substrate and graphene layer is minimized (kGr/Sub ≈ 0) owing to the specific sample preparation method employed. Consequently, the 1 cm^−1^ blue shift of the LB_6,1_ mode in AuNCs/t(3 + 3)LG compared with that of t(3 + 3)LG provides direct evidence for the existence of interfacial coupling at the AuNCs/t(3+3)LG interfaces.

To simulate the LB mode intensities in AuNCs/t(3 + 3)LG using the E-IBPM, additional interlayer bond parameters are required. Figure 4b2 illustrates the schematic of AuNCs/t(3 + 3)LG with the relevant bond parameters: ηAu/Gr (AuNCs/graphene interface), ηt1 (twist interface), ηt2 (adjacent to twist interface), and ηGr (within NLG). To preserve the structural symmetry of t(3+3)LG, the polarizability derivatives \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{l,xx}^{{\prime} }$$\end{document} are parameterized as follows: \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{1,xx}^{{\prime} }={\eta }_{{\rm{Au/Gr}}}-{\eta }_{{\rm{Gr}}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{2,xx}^{{\prime} }={\eta }_{{\rm{Gr}}}-{\eta }_{{\rm{t2}}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{3,xx}^{{\prime} }={\eta }_{{\rm{t2}}}-{\eta }_{{\rm{t1}}}$$\end{document} , and symmetrically for the 4th–6th layers. As shown in Fig. 4b1, the experimental LB mode intensities are well reproduced using parameters ηAu/Gr = 1, ηt1 = − 1, ηt2 = 0.8, ηGr = 0.1, combined with the FDTD-simulated ELoc,l,x (Fig. 3b). The hierarchy ∣ηt1∣≈∣ηt2∣≫∣ηGr∣indicates strong polarizability modulation near the twisted interface under resonance excitation. The negative sign of ηt1 is attributed to the imaginary part of the polarizability under VHS resonance^38,43^. The validity of these parameters is further supported by the excellent agreement between simulation and experiment for both t(3 + 3)LG and t(2 + 2)LG by IBPM and AuNCs/t(2 + 2)LG by E-IBPM, using the same interlayer bond parameters between pristine graphene layers and at twisted interfaces as those in E-IBPM for AuNCs/t(3 + 3)LG (see Supplementary Note 7 and Fig. S7). Notably,∣ηAu/Gr∣≫∣ηGr∣highlights the dramatic plasmon-induced enhancement at the AuNCs interface. This quantitative match across the above distinct material systems (i.e., insulator hBN and tMLG) validates that our theoretical approach correctly captures the essential physics of the PERS of LB modes, implying that the proposed PERS mechanism and E-IBPM are broadly applicable across different AuNCs/2DMs and AuNCs/vdWHs structures, independent of the specific materials.

Furthermore, AgNCs synthesized under varying conditions with different morphologies and densities also exhibit a distinct plasmon resonance peak at 532 nm in their dark-field scattering spectra (see Fig.S11b), which suggests the potential to extend the excitation wavelength range of this approach. As examples, Fig. 4c presents the observed LB modes in AgNCs/30LG and AgNCs/50L-hBN under 532 nm excitation, along with the corresponding E-IBPM simulation results. For AgNCs/30LG, the frequencies and the relative Raman intensities of the observed plasmon-enhanced LB modes under 532 nm excitation can be accurately described using the same set of LB force constants ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\rm{Gr}}}^{{\rm{1st}}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{\,{\rm{Gr}}}^{{\rm{2nd}}}$$\end{document} , kAg/Gr and kGr/Sub) and relevant bond parameters (ηAg/Gr = 10ηGr) as those used under 633 nm excitation. The consistency of these parameters across the two excitation wavelengths indicates strong interactions between AgNCs and NLG under both plasmon resonance excitations. For AgNCs/50L-hBN, the frequencies of the plasmon-enhanced LB modes can be well fitted with khBN = 9.9 × 10^19^ N m^−3^^27^, khBN/Sub = 2khBN and kAg/hBN = 0.02khBN. Remarkably, the relative Raman intensities of these observed LB modes are also well reproduced by the E-IBPM using ηAg/hBN = 10ηhBN, which also suggests a strong interaction between AgNCs and hBN. These results fully validate the feasibility of a wavelength-tunable LB mode detection strategy achieved by utilizing different plasmonic nanocavities. This flexibility allows for the selection of appropriate metal plasmonic cavities for specific material systems, thereby effectively mitigating potential photoluminescence interference.

Discussion

In summary, we have developed a PERS approach using AuNCs or AgNCs to characterize LB modes in NLG, hBN, and related vdWHs. This technique successfully overcomes the detection limitations imposed by weak EPC in conventional Raman scattering. Combining experimental and theoretical analysis, we elucidate the fundamental enhancement mechanism: the synergistic interplay between interfacial polarization modulation at AuNCs/2DM (AgNCs/2DM) interface and the localized electric field enhancement of AuNCs (AgNCs) dominantly amplifies the Raman dipole moment of the 2DM layer at the interface, contributing to the Raman intensity of the LB modes. The proposed E-IBPM provides a quantitative understanding of the relative LB mode intensities across different Au(Ag)NCs/2DMs and Au(Ag)NCs/vdWHs structures.

Our findings establish: (1) A highly sensitive PERS methodology for probing interlayer phonons and interfacial coupling in 2DMs; (2) A PERS mechanism in 2DMs coupled to AuNCs (AgNCs) and a universal theoretical framework (E-IBPM) for understanding plasmon-coupled interlayer vibrational modes. As interlayer vibration modes carry unique information about the interlayer or interfacial properties in 2DMs and related heterostructures, the plasmonic nanocavity-enabled LB modes in such systems provide a unique, non-destructive tool for characterizing the structure and interlayer coupling at hidden interfaces in complex, multilayer stacks, which are inaccessible to many other techniques. Our work highlights that in complex plasmonic cavity-coupled van der Waals systems, the relative Raman intensity serves as a more reliable probe of the underlying PERS mechanisms than the hard-to-quantify absolute enhancement factors, an insight that is broadly applicable to the field. In addition, the plasmonic nanocavities in this work act as a universal amplifier, capable of not only enhancing intrinsic dipole transitions but also converting weak processes that indirectly couple to the field via polarizability changes into strong, observable signals. The methodology holds potential for application to a broader range of elusive quasiparticles, with interlayer excitons and certain plasmonic resonances in vdWHs being prime candidates. This work opens avenues for precise characterization of interlayer interactions in emerging two-dimensional devices while providing fundamental insights into plasmon-phonon coupling phenomena.

Materials and methods

Sample preparation

Multilayer graphene and hBN flakes were mechanically exfoliated from their bulk crystals and transferred onto 90 nm SiO_2_/Si substrates with oxygen plasma-treated surfaces to improve adhesion. The layer number of graphene was determined via Raman spectroscopy by comparing the intensity ratio of the Si peak (~520.7 cm^−1^) from the substrate beneath the sample to that from an adjacent bare SiO_2_/Si substrate under the same focusing condition^23^. The thicknesses of hBN flakes were measured by atomic force microscopy and further verified using the frequencies of the plasmon-enhanced LB modes. The tMLG and hBN/tMLG vdWHs were assembled via a standard polycarbonate-assisted dry-transfer technique. Finally, an 8 nm-thick Au film was thermally evaporated onto the 2DMs and vdWHs at a deposition rate of 0.5 Å/s.

Structural characterization

The surface morphology of the plasmonic nanocavities was characterized by SEM (Helios 5 CX DualBeam, Thermo Scientific) with an accelerating voltage of 3 kV. To investigate the cross-sectional interfacial structure, the samples were prepared by the focused ion beam (FIB) lift-out technique and subsequently analyzed by high-resolution TEM (Talos F200X, Thermo Scientific) at an accelerating voltage of 200 kV.

Dark-field scattering spectroscopy

The dark-field scattering measurements were performed using a home-made SmartRaman confocal micro-Raman module (Institute of Semiconductors, Chinese Academy of Sciences) coupled with a dark-field microscope and a Horiba iHR320 spectrometer. The white light from a tungsten lamp was directed to a 100 × objective lens (numerical aperture, N.A. = 0.9) to form uniform dark-field illumination. The scattered signal of the samples was collected by the same objective, which was then directed to the iHR320 spectrometer.

Raman spectroscopy

Raman measurements were performed using a Jobin-Yvon HR800 system equipped with a charge-coupled device (CCD). The Raman scattering signal was collected through a 100 × objective lens (N.A. = 0.9). Two laser wavelengths were used in the experiments, including 633 nm from a He-Ne laser and 532 nm from a solid-state laser. The laser plasma lines were removed using a BragGrate^TM^ bandpass filter (OptiGrate Corp.), and low-frequency Raman measurements down to 5 cm^−1^ were achieved with three BragGrate^TM^ notch filters (OptiGrate Corp.). The typical spectral resolution is 0.34 cm^−1^ per CCD pixel when using an 1800 lines/mm grating with 633 nm excitation. The laser power was kept below 1 mW to prevent sample heating.

Simulations

The local field generated by AuNCs was investigated using the three-dimensional finite-difference time-domain (3D-FDTD) simulation in Lumerical software. The incident light was modeled as a plane wave with x-polarization. The Bloch boundary condition in the x**y-plane was applied to emulate infinite periodicity, and perfectly matched layers along the z-direction eliminated artificial boundary reflections. A critical uniform grid resolution of 0.1 × 0.1 × 0.1 nm^3^ was implemented to minimize numerical dispersion errors. The dielectric function of Au incorporated a modified Drude-Lorentz formula to account for nanoscale surface effects, and the optical response of graphene was governed by the Falkovsky model embedded in the software. The simulations extended to 1500 fs with dynamic convergence monitoring, ensuring field residuals stabilized within machine precision thresholds. The components of the local electric field were extracted through Fourier analysis of time-domain data at the characteristic wavelength of 633 nm.

Supplementary information

SUPPLEMENTAL MATERIAL for Plasmonic nanocavity-enabled universal detection of layer-breathing vibrations in two-dimensional materials