Quantum Effects in the Acoustic Plasmons of Atomically-Thin Heterostructures

Recent advances in nanofabrication technology now enable unprecedented control over 2D heterostructures, in which single- or few-atom thick materials with synergetic opto-electronic properties can be combined to develop next-generation nanophotonic devices. Precise control of light can be achieved at the interface between 2D metal and dielectric layers, where surface plasmon polaritons strongly confine electromagnetic energy. Here we reveal quantum and finite-size effects in hybrid systems consisting of graphene and few-atomic-layer noble metals, based on a quantum description that captures the electronic band structure of these materials. These phenomena are found to play an important role in the metal screening of the plasmonic fields, determining the extent to which they propagate in the graphene layer. In particular, we find that a monoatomic metal layer is capable of pushing graphene plasmons toward the intraband transition region, rendering them acoustic, while the addition of more metal layers only produces minor changes in the dispersion but strongly affects the lifetime. We further find that a quantum approach is required to correctly account for the sizable Landau damping associated with single-particle excitations in the metal. We anticipate that these results will aid in the design of future platforms for extreme light-matter interaction on the nanoscale.


I. INTRODUCTION
The isolation of monolayer graphene [1] has stimulated extensive research efforts in two-dimensional (2D) materials, due in part to their unique electronic and optical properties, which are well-suited for use in compact, ultra-efficient photonic and opto-electronic devices [2]. Polaritons in 2D materials are particularly appealing in the field of nano-optics because they can confine external electromagnetic fields into extremely small volumes [3], enabling control of quantum and nonlinear optical phenomena on the nanoscale [4]. Additionally, 2D polaritons are extremely sensitive to their surrounding environment, a property that renders them as good candidates for optical sensing [5], but also as enablers of new electro-optical functionalities when different atomically-thin materials are combined to form heterostructures [6,7].
Surface-plasmon polaritons (SPPs) are formed when light hybridizes with the collective oscillations of charge carriers at the interface between dielectric and conducting media [8], offering particularly strong confinement of electromagnetic energy down to truly nanometer length scales [9]. Noble metals are the traditional material platform used in plasmonics research, although they are difficult to actively tune and suffer from large ohmic losses [10,11]. Graphene can help circumvent these limitations, as it supports highly-confined and long-lived plasmonic resonances that can be electrically modulated [12][13][14]. In particular, the encapsulation of exfoliated monolayer graphene (MG) in hexagonal boron nitride (hBN) has been shown to dramatically improve the quality factor * Electronic address: javier.garciadeabajo@nanophotonics.es of plasmon resonances, with measured lifetimes of propagating modes ∼ 0.5 ps at room temperature [15], and even beyond 1 ps at lower temperatures [16]. However, graphene plasmon studies have been so far limited to the terahertz and mid-infrared (mid-IR) spectral domains because the resonance energies of these excitations are severely constrained by the doping densities that can be sustained by the carbon layer, although some prospects have been formulated to extend their range of operation into the near-IR region [17].
Hybrid systems comprising noble metal layers and graphene potentially alleviate the limitations of plasmons in both of these materials by capitalizing on the appealing electrical tunability of graphene combined with the visible and near-IR plasmon resonances of noble metals. For example, tuning the damping of metal plasmons with electrostatically-gated graphene has been shown to enable a small degree of electrical tunability [18], while this approach has been predicted to reach order-unity active modulation if a graphene layer is deposited on a thin metallic film [19]. Still in the mid-IR, the proximity of a metal film in a heterostructure comprising hBNencapsulated graphene and optically-thick metal surfaces has been recently demonstrated to push plasmon confinement down to a few atomic lengths outside the graphene sheet by rendering the plasmon dispersion acoustic, running close to the continuum of graphene intraband excitations [20,21]. Indeed, acoustic plasmons can result from the hybridization of modes in two neighboring metal surfaces [22], or also in graphene and a metal surface, emerging as a low-energy branch in the energymomentum dispersion diagram [23]. In contrast to the intuition that metal screening quenches the graphene response, an acoustic plasmon branch has been predicted to manifest itself even when graphene is directly deposited arXiv:1901.07098v3 [cond-mat.mtrl-sci] 24 Jan 2019 on the metal without a spacer [24]. For atomic-scale graphene-metal separations, the nonlocal nature of the optical response in both the metal and graphene layers plays an important role, demanding a rigorous theoretical treatment to accurately describe the dispersion and lifetime of acoustic plasmons [25]. The hydrodynamic model [26][27][28][29] has been extensively employed in this context, consisting in describing thin metal films within the framework of classical electrodynamics and introducing conduction electrons as a classical plasma [30][31][32][33][34]. However, the hydrodynamical model neglects effects associated with the electronic band structure that can become important when optical confinement reaches the atomic length scale [35][36][37][38][39]. First principles simulations have also been employed based on time-dependent densityfunctional theory to characterize the plasmons of metal films [40][41][42][43][44][45][46], but they involve computationally expensive simulations that are difficult to extend to systems comprising heterostructures and lacking a single overall atomic periodicity. Nevertheless, a more rigorous study is still pending on the nonlocal effects that affect the plasmons supported by thin metal films and their hybridization with graphene modes.
Here, we explore the role played by quantum and finitesize effects in the acoustic plasmons resulting from hybridization of MG with few-atom-thick metallic films. We simulate the graphene and the metal within the random-phase approximation (RPA), incorporating in the metal phenomenological information on the main characteristics of the conduction electronic band structure, such as surface states, electron spill out, and a directional band gap associated with bulk atomic-layer corrugation. Comparison among different levels of approximation to the electron confinement potential, as well as to classical electromagnetic approaches, reveals that quantum and finite-size effects produce dramatic changes in the plasmon lifetime, which require a quantum description of the system to be properly accounted for. Remarkably, we find that a single metal monolayer can render the graphene plasmon dispersion acoustic, while the addition of more monolayers hardly modifies the dispersion but does change the lifetime. Also, the plasmon lifetime is strongly affected by Landau damping in the metal, which is automatically incorporated in the RPA description. Our study and methodology are of general use for the description of plasmons in thin metal films and their interaction with nearby 2D materials.

II. RESULTS AND DISCUSSION
Acoustic plasmons result from the repulsion between modes in closely spaced metallic films [8] and their name derives from the linear relation between their frequency ω and in-plane wave vector Q, similar to acoustic waves. A tutorial description of this concept is provided by considering metal films of small thickness d m and local permittivity m (ω), the response of which can be approx-imated as that of a zero-thickness layer of 2D conductivity [48,49] . When deposited on the planar surface of a dielectric of permittivity d , the reflection coefficient for p polarization (i.e., the one associated with surface plasmons) reduces to r = (1 − d + ξ)/(1 − d + ξ) (see Appendix), where ξ = 4πiσQ/ω. Now, the plasmons of two films separated by a dielectric of thickness d d and permittivity d are determined by the Fabry-Perot condition r 2 e −2Qd d = 1, where the exponential represents the round trip for signal propagation across the dielectric spacer. Putting these elements together, we obtain two plasmon branches described to the expression . In the Drude model for the metal conduction electrons, we can further approximate the metal permittivity as m ≈ b − ω 2 p /ω 2 , where ω p is the bulk metal plasma frequency and b is a background permittivity accounting for polarization of inner electrons. This leads to a dispersion relation ω ≈ ω p / b − 1 + F ± /(Qd m ), which we plot in Fig.  1 (blue-solid curves) for silver films (i.e., b = 4 and ω p = 9.17 eV, as obtained by fitting optical data [47]) of small thickness (N = 1 and 5 Ag(111) atomic layers, with 0.236 nm per layer) separated by a silica film ( d = 2.1) of varying thickness d d . We find an upper plasmon branch and an acoustic plasmon at lower energies. In the small Q limit, the above expressions can be further approximated as ω/ω p ≈ √ Qd m for the upper branch (i.e., with a characteristic ω ∝ √ Q dependence) and ω/ω p ≈ Q d d d m /(2 d ) for the acoustic plasmon branch (Fig. 1, blue-dashed curves). This simple tutorial model is in reasonable agreement with a local dielectric description of the system (Fig. 1, black-dashed curves) and it correctly captures the increasing degree of mode repulsion as d d is reduced. In this work, we produce a more realistic quantum-mechanical simulation (see below), which we anticipate to yield a dispersion relation in remarkably excellent agreement with the above picture ( Fig. 1, color density plots, representing the loss function Im{R}, where R is the reflection coefficient of the metaldielectric-metal structure for p polarization). The above discussion can be readily extended to metal-graphene and double-layer-graphene structures, yielding equally good agreement for the dispersion relations [23,50]. Nevertheless, as we show below, an accurate calculation of the acoustic plasmon lifetimes is not possible without incorporating quantum-mechanical elements in the description of the system.
We are interested in studying the plasmons supported by hybrid planar films comprising an atomically thin metal layer and MG. Due to translational symmetry along the plane of the structures, each plasmon oscillating at an optical frequency ω can be characterized by a well-defined in-plane wave vector Q, so that its associated electromagnetic field depends on time and inplane coordinates R = (x, y) just through an overall factor e iQ·R−iωt , which is implicitly understood in what follows. We express the response of the hybrid film in  terms of reflection and transmission coefficients for each of its constituting layers in a Fabry-Perot fashion (see Appendix). The calculation is simplified by the fact that the wavelengths of the plasmons under consideration are small compared with the light wavelength at the same frequency, therefore allowing us to work in the quasistatic limit, in which s-polarization components do not contribute.
We describe graphene through its RPA conductivity [51,52] (see Appendix), which, in virtue of the van der Waals nature of its binding to the surrounding materials, should not be affected by the electronic properties of the latter, apart from some possible doping associated with carrier transfer. This level of description has been shown to excellently describe the optical response of graphene down to atomic length scales [53]. In the main text we present results for graphene directly deposited on metal, while in the Appendix we provide additional calculations assuming a hBN intermediate layer, which is treated as a local, anisotropic dielectric film (see Appendix); this approach should capture the main modifications produced by the hBN layer on the response of the structure, which mainly consist of a featureless dielectric screening, accompanied by the signatures imprinted by its mid-IR phonon polaritons.
The response of the metal layer is however sensitive to the proximity of the graphene. As we show below, this demands an accurate description of its nonlocal re-sponse. To this end, we adopt a QM model for the metal layer consisting in calculating its transmission and reflection coefficients within the RPA. Full details of the formalism are provided in the Appendix. In brief, the response of valence electrons is incorporated by calculating the non-interacting RPA susceptibility χ 0 from the one-electron wave functions, while polarization of inner bands is accounted for through a local screened interaction. We further assume in-plane translational symmetry, so that the valence electron wave functions can be written as e ik ·R ϕ j (z), labeled by the 2D in-plane electron wave vector k and the out-of-plane band index j. The wave function component ϕ j (z) is obtained upon solution of the one-dimensional (1D) Schrödinger equation along the out-of-plane direction, specified for a confining potential V (z). We focus on gold and silver, for which a choice of effective mass m * = m e is appropriate, and more precisely, we apply this procedure to films consisting of a finite number N of either Au(111) or Ag(111) atomic layers (metal thickness d m = N a s , with interlayer spacing a s ≈ 0.236 nm). The results of this procedure depend on the specific choice of potential V (z), for which we use three different levels of approximation: • Atomic-layer potential (ALP). We adopt from Ref.
[54] a model potential that incorporates a harmonic corrugation in the bulk region and a smooth density profile at the surface, with parameters such that several important features of the electronic structure are correctly reproduced in the semiinfinite metal limit: the work function, the surfaceprojected bulk gap, and the position of the surface states relative to the Fermi level (see Fig. 6 in the Appendix, where the semi-infinite limit is approached for N ∼ 100). This potential therefore incorporates phenomenological information in a realistic fashion on (1) the out-of-plane quantization of electronic states, (2) the spill out of the electron wave functions beyond the metal film edges, and (3) the surface-projected electronic band gap produced by the bulk atomic-plane periodicity. Actually, this potential also describes surface and image states [54], and therefore should realistically account for electron spill out effects. In practice, we designate the positions of the thin film surfaces as z = 0 and z = d m , so that the first and last atomic planes are located at z = a s /2 and z = d m −a s /2, respectively (see panel (i) in Fig. 2a).
• Finite-barrier model (FBM). Neglecting the atomic periodicity of the metal, we assign the potential Fig. 2a), extracting the potential barrier depth V 0 from the ALP for each metal. This choice of V 0 should provide the correct offset for the obtained electron energies. The FBM describes (1) 1D quantization of electronic states and (2) electron spill out beyond the metal film edges.
• Infinite-barrier model (IBM). To further simplify the model, we consider an infinite potential well The IBM accounts only for the quantization of electronic states along the film confinement direction.
In Fig. 2b we show the SPP dispersion of self-standing gold films as predicted in the IBM, FBM, and ALP models, indicated by dotted, dashed, and solid curves, respectively, for gold films consisting of N = 1-10 atomic layers. In each case we plot the frequency corresponding to the maximum of the loss function Im{r m }, where r m is the reflection coefficient for p polarization (see Appendix) at a given in-plane optical momentum Q. We obtain excellent agreement among the plasmon dispersions described by each choice of binding potential V (z) within the QM model, even down to atomic-monolayer gold. We note that the dispersion curves predicted in the ALP are slightly redshifted with respect to those of the FBM, which in turn are redshifted relative to the dispersion of the IBM. Presumably, this redshift is related to electron spill out [28,55], although an interplay between spill out and interband polarization is known to control the actual sign of the plasmon shift for small wave vectors [56]. The plasmon lifetime in the gold film is characterized by its spectral full width at half maximum (FWHM) ∆ω plotted in Fig. 2c, as determined from the linewidth of the Im{r m } spectral curves; the results, which are rather independent of the choice of potential model, lie near the intrinsic phenomenological width γ = 71 meV introduced in the RPA formalism for gold, a value taken from a Drude model fit of optical data [47]. We find a qualitatively similar behavior in sil-ver films, but with smaller ∆ω than in gold (see Fig. 7a in the Appendix).
By depositing doped MG on the thin metal films, we introduce low-energy acoustic plasmon modes in the dispersion relation, for which metal screening confines light in the region between the graphene and the metal [57]. In Fig. 2d we plot the calculated acoustic-plasmon dispersion relations for a graphene Fermi energy E F = 1 eV obtained by using the three QM models considered in Fig.  2b to account for the metal. Interestingly the acoustic plasmons show a dispersion that is rather independent of metal thickness down to N = 1, apart from a minor shift towards the graphene intraband region (ω ≤ v F Q, where v F ≈ 10 6 m/s is the Fermi velocity in graphene) with increasing N . A single atomic metal layer is thus capable of pushing the graphene plasmon dispersion (dashed curve) toward this region and render it acoustic. The associated FWHM (Fig. 2e) of the acoustic plasmon is considerably reduced compared to its higher-energy bulk counterpart (Fig. 2c) at similar wave vectors although (unlike the bulk branch the acoustic plasmon linewidth) it strongly depends on the metal thickness and the model used for the potential V (z). It is remarkable that acoustic plasmons exist even by directly depositing the graphene on the metal without spacing, which confirms a recent prediction of this effect [24], defying the intuition that metal screening can quench the graphene response. At low mid-IR frequencies, the acoustic plasmons are predicted to exhibit smaller lifetimes as the film thickness increases; thicker metal films are thus providing more efficient screening and less relative weight of the field inside the metal, where the intrinsic inelastic rate is larger than in the graphene alone (cf. dashed curve, corresponding to the assumed intrinsic lifetime τ = 500 fs in graphene). At higher energies, the results become more involved, as metal screening is less effective, and the FWHM depends more strongly on the model used for the potential. We remark that the ALP model, which is expected to provide the most accurate description because it incorporates several phenomenological features of the electronic band structure, leads to generally higher damping than the other models.
In Fig. 3 we study the full plasmon dispersion of the MG-gold film heterostructure, which is schematically illustrated in Fig. 3a for a system comprised of only two atomic gold layers (see also Fig. 8 in the Appendix for contour plots of the loss function Im{R}). Fig. 3b indicates that, while the dispersion of the acoustic plasmon is only marginally influenced by the film thickness, the higher-energy plasmon retains a similar dependence as that of the isolated film considered in Fig. 2b. At high energies approaching the interband damping regime in graphene ( ω ≥ 2E F − v F Q), the high-energy plasmon dispersion undergoes a slight redshift, accompanied by an increase in the effective plasmon damping (see Fig.  3c).
It is remarkable that the presence of a single atomic metal layer (N = 1) is sufficient to support plasmons (see Figs. 2d and 3b), with a group velocity at low frequencies v g ≈ 1.7×10 6 m/s for E F = 1 eV that is nearly unchanged by the addition of further atomic metal layers. These modes can be modulated by varying the Fermi energy E F in MG, with relatively small damping in the ω E F region, while strong attenuation is produced by interband transitions at higher energies (Fig. 3d,e).
Examination of the effective damping rates in Figs. 2e and 3c indicates that an increase in film thickness can reduce the plasmon linewidth, owing to more effective screening by the metal layer. We are thus prompted to consider acoustic plasmons in the MG-gold film heterostructure in the limit of a semi-infinite metal layer, which is schematically illustrated in Fig. 4a. Our QM treatment of a thin film can also describe the response of a semi-infinite film if we approach the N → ∞ limit, which in practice is reached for N ∼ 100 atomic layers (see Fig. 6 in the Appendix).
We find it instructive to compare the above QM descriptions of the metal to classical approaches that are extensively used in recent literature. In particular, we consider the so-called specular-reflection model [58] (SRM), also known as semiclassical infinite barrier model [59], which allows us to express the response of a semi-infinite material in terms of the nonlocal dielectric function of the bulk, while a straightforward generalization can deal with arbitrary shapes [60]. Following the detailed prescription presented in the Appendix, we apply this model to three different choices for the bulk dielectric response function: the hydrodynamical model [26,27], the Lindhard dielectric function [61,62], and the Mermin prescription [63]. The bulk hydrodynamic model combined with the SRM coincides with the hydrodynamic model for finite geometries, which has been extensively used to discuss nonlocal effects in nanostructured metals [29,60,64]; it incorporates nonlocal effects though the hydrodynamic pressure, as well as a phenomenological damping rate γ. The Lindhard dielectric function is used here with a correction intended to account for d-band screening [60] (see Appendix), as well as a damping rate γ effectively introduced by replacing ω by ω+iγ/2 in the Lindhard formula. The Mermin prescription corrects the Lindhard formula in order to preserve local conservation of electron density during damping processes [63].
In Fig. 4b,c we compare the performance of the classical (i.e., SRM using hydrodynamic, Lindhard, and Mermin dielectric functions) and QM (i.e., RPA with ALP, FBM, and IBM potentials) approaches for MG deposited on an optically-thick gold surface. The dispersion relation is very similar in both cases. In contrast, we find significant discrepancies in the linewidth: the classical approach overestimates plasmon broadening when using the Lindhard and Mermin dielectric functions, but it produces a severe underestimate in the hydrodynamical model, presumably because the latter does not account for the efficient mechanism of Landau damping associated with decay into metal electron-hole pair (e-h) excitations; indeed, the bulk e-h region in gold (orange area  (e)). (f-i) Spatial distribution of the in-plane (Ex) and out-of-plane (−iEz) components of the plasmon electric field for various metal thicknesses (see legend in (f)) at the wave vector indicated by the black arrows in (b) (i.e., Q = 0.1 nm −1 ). Solid and dashed curves represent the fields for the high-energy and acoustic plasmons, respectively (see (b)). The graphene layer is approximated as a zero-thickness film (vertical dashed lines, separated dg/2 from the metal, where dg = 0.33 nm is the nominal graphene thickness); the near fields of the acoustic plasmons (h,i) are mainly concentrated in the graphene-metal region, whereas the higher-energy plasmons (f,g) exhibit more delocalized field profiles. We describe graphene in the RPA in all plots. in Fig. 4b, defined by ω ≤ ( /m * )(Q 2 /2 + Qv Au F ), where the gold Fermi velocity v Au F ≈ 1.4 × 10 6 m/s exceeds by 40% that of graphene) overlaps the plasmon dispersion for ω 0.8 eV, where ∆ω takes relatively large values. In addition to this, the finite damping introduced in the model and the lack of translational symmetry along the surface normal direction effectively extend coupling to the e-h region toward low in-plane wave vectors, thus increasing its overlap with the plasmon. Remarkably, the plasmon damping predicted by the QM approach using the IBM potential agrees well with the SRM approach using the Lindard and Mermin dielectric functions; the main difference between them lies in the neglect of quantum interference between outgoing and surface-reflected electron components in the SRM (see Appendix), which seems to play a minor role. Finally, in Fig. 4c,d we consider optically-thick silver surfaces, which exhibit qualitatively similar behavior as gold, although the reduced intrinsic inelastic damping of the argent metal ( γ = 21 meV) results in smaller plasmon linewidths.
We obtain qualitatively similar results when MG is separated from the metal by a thin layer of hBN (see Figs. 9 and 10 in the Appendix). The latter introduces sharp spectral features characterized by avoided crossing between the plasmons and the mid-IR photons of the insulator. Additionally, the MG-metal interaction is reduced relative to the structures without hBN because of their larger separation, therefore yielding low-energy plasmon bands with a less acoustic character. This effect is clearly observed when comparing contour plots of the loss function Im{R} in both types of structures (cf. Figs. 8 and 11 in the Appendix). Interestingly, although the acoustic plasmon dispersion depends strongly on the spacer thickness and level of graphene doping, it is not so sensitive to the permittivity of the dielectric spacer (Fig. 12 in the Appendix).
An interesting scenario arises when graphene is fully embedded in metal. We can easily simulate this type of structure through a trivial extension of the Fabry-Perot approach discussed in the Appendix. The results (Fig.  5, for graphene doped to EF = 1 eV and embedded in silver) reveal an interplay between three different plasmon modes provided by the graphene and each of the two metal regions. In particular, we find again an acoustic plasmon that is mainly localized in the graphene, now pushed even further toward the intraband transition region. Additionally, a higher-energy acoustic plamon emerges, and this mode becomes more acoustic and confined when one of the lower metal region is a single atomic layer (cf. Figs. 5b and 5c) because the plasmons in the two metal films are then closer in energy, so they undergo stronger interaction, similar to what we observed in Fig. 1. This second acoustic plasmon exhibits larger tunability with the number of metal atomic layers, although it is more lossy than the one dominated by graphene, as it has more weight in the metal, which is taken to have a lower intrinsic inelastic lifetime (31 fs for silver [47]) than graphene (500 fs). Similar conclusions are observed for lower graphene doping (see Fig. 13 in the Appendix, where we consider E F = 0.2 eV).

III. CONCLUSION
The acoustic plasmons supported by monolayer graphene when it is deposited on noble metal surfaces are strongly influenced by the metal optical response. We have identified quantum finite-size effects in the optical response of ultrathin films and semi-infinite metal layers that impact the acoustic plasmons formed by interaction with doped graphene. Remarkably, a single atomic metal layer is capable of rendering the graphene plasmons acoustic, displaying a group velocity ∼ 1.6 v F for a doping E F = 1 eV that is nearly insensitive to the addition of more metal layers. The level of metal screening, which influences the plasmon lifetimes, is however strongly dependent on metal thickness. We further reveal the important contribution of electron-hole pairs in the metal to the plasmon damping, a fair description of which can only be achieved by incorporating details of the metal electronic band structure, electron spill out, and surface-projected bulk gaps, such as we do in this work using a computationally efficient quantum approach. The present study could be improved by using a first principles estimate of the internuclear distance in the MG-metal interface, as well as by introducing the effect of in-plane atomic corrugation, although these are computationally demanding calculations that we expect to produce only minor corrections. We anticipate that the results presented here can elucidate the role of quantum and finite-size effects in acoustic plasmons, enabling a fundamental understanding of the damping and extreme spatial confinement associated with these excitations.
We consider a multilayer hybrid film consisting of MG deposited on a metal film, from which it is separated by a thin dielectric layer. The optical response is expressed in terms of the reflection (r i ) and transmission (t i ) coefficients of the graphene layer (i = g), the dielectric spacer (i = d), and the metal film (i = m) by describing the heterostructure as a Fabry-Perot resonator characterized by the total reflection and transmission coefficients [65] R = r gd + t gd t dg r m 1 − r dg r m , where the doubly-indexed coefficients combining graphene and the dielectric layer are defined as This prescription is convenient because we can use coefficients calculated for a surrounding vacuum by just taking the graphene-dielectric and dielectric-metal spacings as consisting of a zero-thickness vacuum layer. In the quasistatic limit here considered (see below), only p-polarization coefficients make a nonzero contribution.
Explicit values of these coefficients are given below for graphene and hBN, accompanied by a quantummechanical approach to deal with the metal film, which we compare with a nonlocal classical solution in the thickmetal limit. While in the main text we disregard the dielectric spacer by imposing r d = 0 and t d = 1, we provide in the Appendix complementary simulations including the effect of a thin hBN spacing layer between the graphene and the metal.

Appendix B: Quantum-mechanical description of the metallic film
We study the linear optical response of a self-sustained metal film of thickness d m contained in the 0 < z < d m region within the RPA [66], for which we construct the noninteracting susceptibility χ 0 in terms of one-electron wave functions, which are in turn calculated as the solutions of an effective potential V (z). For simplicity, we assume translational invariance along the in-plane directions (coordinates R = (x, y)) and adopt different models for the out-of-plane potential profile (see Fig. 2a). This allows us to express the optical response in terms of components of well-defined in-plane parallel wave vector Q and frequency ω, so we assume an implicit dependence on these variables, as well as a multiplicative factor e i(Q·R−ωt) . For example, the potential in real space-time has the form φ(r, t) = (2π) −3 d 2 Q dω e i(Q·R−ωt) φ(Q, z, ω). Additionally, we ignore retardation effects because the wavelengths of the plasmons under investigation are much smaller than the light wavelength, and consequently, we study the response using electrostatic potentials. Specifically, we consider an external (incident) potential φ ext (z) = (2π/Q)e −Q|z−z0| representing a source located below the film at z = z 0 (in practice we take z 0 = 0). The reflection coefficient is defined as the ratio of induced to external potentials r m = 1 − φ(z 0 )/φ ext (z 0 ), where an overall minus sign is introduced to make this definition coincide with the quasistatic limit of the Fresnel coefficient for p polarization. Likewise, the transmission coefficient is defined as t m = φ(z 1 )/φ ext (z 0 ), where the transmitted potential is evaluated at a position z 1 right above the film (we take z 1 = d m ).
The response of inner shells is accounted for by introducing a homogeneous film of local background permittivity b (ω) contained in the 0 < z < d m region, so that it extends half an atomic-layer spacing beyond each of the two outermost atomic planes (see ALP below). We note that b can take relatively large values in noble metals within the visible and near-IR spectral ranges (e.g., ∼ 9 in Au) due to d-band polarization. In practice, we obtain b (ω) by subtracting from the experimentally measured metal dielectric function [47] exp (ω) a Drude term, such that b (ω) = exp (ω) + ω 2 p /ω(ω + iγ), where ω p = 9.06 eV, γ = 0.071 eV for gold and ω p = 9.17 eV, γ = 0.021 eV for silver (see Fig. 14 in the Appendix for plots of exp (ω) and b (ω)).
We express the total potential as where the φ ext b term accounts for the potential created by φ ext in the presence of the background slab of permittivity b (i.e., including the background response, but not the response of the conduction electrons), while the integral gives the contribution due to the induced-charge density ρ ind (z ) associated with disturbances in the conduction electrons. The latter is mediated by the screened interaction v b (i.e., the Coulomb potential created at z by a point charge at z oscillating with frequency ω in Q space, including the effect of background polarization). In the absence of background polarization (i.e., for b = 1) we have v b (z, z ) = (2π/Q)e −Q|z−z | , while in the presence of an b = 1 film we still find a closed-form expression by direct solution of Poisson's equation [67]: The v dir b term captures the charge singularity in the interaction potential within each homogeneous region of space, whereas the addition of the v ref b term guarantees the continuity of the potential and the normal displacement at the interfaces. Using this expression, and noticing that the external potential originates in a source at z 0 ≤ 0, we can readily write the external potential including the interaction with the background film as Assuming linear response, we can express ρ ind in terms of the susceptibility χ according to where the external potential φ ext b driving the free electrons has been corrected by the direct background contribution as explained above. Using matrix notation with z acting as an index and a dot indicating integration over this coordinate, we can write where χ 0 is the non-interacting RPA susceptibility [66]. We use the well-known result [62] for the full spatial dependence of this quantity in terms of the one-electron metal wave functions ψ i (r), where the factor of 2 accounts for spin degeneracy, f i is the occupation of state i with energy ε i , and γ denotes a phenomenological inelastic damping rate [62,66]. In order to exploit translational invariance in the film, we multiplex the state index as i → {j, k }, where j labels eigenstates of the z-dependent out-of-plane 1D Schrödinger equation defined by the noted binding potential V (z) with an effective mass m * (see below), while k runs over 2D in-plane electron wave vectors. The electron wave functions have therefore the form ψ i (r) = A −1/2 e ik ·R ϕ j (z), where A is the film normalization area and ϕ j (z) are eigenstates of the 1D problem. Using these wave functions, we can recast Eq. (B2) as is the quantity actually used in Eq. (B1). Here, we define occupation numbers f j,k = Θ(E F − ε j − 2 k 2 /2m * ) that follow the Fermi-Dirac distribution at zero temperature for a Fermi energy E F , where the rightmost term inside the step function describes a parabolic dispersion along in-plane directions with effective mass m * as well.
We evaluate these expressions by expanding all quantities in sine-Fourier transform within an embedding infinitepotential box spanning a ∼ 1 nm vacuum region on each side of the film. The quantities χ 0 , v b , and χ then become square matrices in this representation, so we operate with them using linear algebra techniques. We have further corroborated the accuracy of our numerical results by comparing with a real-space discretization in the z coordinate, which further reduces spurious Gibbs' oscillations in the calculation of the near fields, although the sine basis generally converges with a smaller number of elements.
The average volumetric electron density in the film is n = [2/(Ad m )] j k f j,k , where the factor of 2 is again due to spin degeneracy. Using the customary transformation k → (2π) −2 A d 2 k , we obtain the self-consistent expression for the Fermi energy where M denotes the highest band index for which ε M < E F . To correctly define E F for arbitrary thickness, we fit the charge density n in the bulk limit (d m → ∞) by imposing experimentally-measured values of E F relative to vacuum in gold (−5.5 eV) and silver (−4.65 eV). The choice of potential parameters [54] in the ALP further ensures that the binding energies of surface states agree with available experimental measurements [68] (see thickness-dependent band structure in Fig. 6 in the Appendix). Under these conditions we obtain effective charge densities n = 70.5 nm −3 for gold and n = 59.6 nm −3 for silver. For the sake of consistency, we use these values of n in the calculations of χ 0 presented in this work for all models of V (z) (see Fig. 2a), and we also assume the same density per layer for finite-thickness films.
Appendix C: Nonlocal classical electromagnetic description of a semi-infinite metal film For the sake of comparison in the thick-film limit, we adopt the SRM [41,58,60] to compute a nonlocal reflection coefficient r m for a semi-infinite metal in the framework of classical electrodynamics. This model allows us to relate the surface response to the momentumand frequency-dependent dielectric function of the metal m (q, ω) under the assumption that conduction electrons undergo specular reflection at the surface without quantum interference between outgoing and reflected components. In practical terms, we calculate the reflection coefficient in the SRM as [58,60] plays the role of a surface response function. We further decompose the bulk metal permittivity as [60] m (q, ω) = b (ω) + free (q, ω) − 1, where b is the local background response defined above (see Fig. S9 in the Appendix) and free describes the nonlocal contribution of free conduction electrons. We adopt three different levels of approximation for the latter in the calculations presented in Fig.  4: where the subscripts refer to the hydrodynamic [26,27], Lindhard [61,62], and Mermin [63] models (see main text). Here, ω p is the bulk plasma frequency of the metal, E F = 2 k F 2 /2m * is the Fermi energy, k F = (9π/4) 1/3 (1/r s ) is the Fermi wave vector expressed in terms of the one-electron radius r s , and we use the function Incidentally, the reflection coefficient derived from the SRM combined with the choice free = hydro coincides with the solution of the classical hydrodynamic equations [69], where the β = 3/5 k F term accounts for the hydrodynamic pressure. We apply these models to gold and silver by taking m * = m e and r s = 3×Bohr-radius (i.e., r s = 0.16 nm), as well as values for ω p and γ as specified above.

Appendix F: Atomic layer potential (ALP)
In the ALP model we use the parametrized potential of Ref. [54] to obtain the one-electron states of a metal film including ad hoc band-structure information. Specifically, this potential consists of a harmonic bulk component, a differentiated region describing each outermost layer, and a long-range image tail, constructed in such a way that it reproduces several experimentally observed electronic structure features, namely: the work function, the surface-projected bulk gap, and the position of the surface states relative to the Fermi level, all of which depend on material and crystallographic orientation. For simplicity and to a good approximation, we take the effective electron mass as m * = m e in all directions. For a semi-infinite metal placed in the z > 0 region, the potential referred to the vacuum level can be written as [54] V surf (z) = where the normal coordinate z is given relative to the position of the outermost atomic plane (z = 0); a s is the inter-atomic layer spacing; the coefficients A 1 and A 10 are chosen to reproduce the width and position of the noted energy gap, respectively; the space between z = z i and z = 0 represents the transition from the solid bulk to free-space, where electron spill out takes place; and the parameters A 2 and η determine the positions of the Fermi level and the surface states relative to the vacuum level. We list the values of these parameters for Au(111) and Ag(111) in Table I. Five of the remaining six parameters are determined by imposing the continuity of the potential and its first derivative, so that A 20 = A 2 − A 1 − A 10 , A 3 = −A 20 + A 2 cos (ηz t ), α = ηA 2 sin (ηz t )/A 3 , λ = 2α, and z i = α −1 log (−λ/4A 3 ) + z t , while the intermediate point z t = −5π/4η is fixed with respect to the surface atomic layer in such a way that z i corresponds to the image plane position, which is important for describing image states [54]. The potential for a film with outermost atomic planes at z = a s /2 and z = d m − a s /2 can be expressed using Eq. (F1) as which is continuous at z = d m /2 by construction. Obviously, d m must be taken to be a multiple of the atomiclayer spacing a s .