Optical response of noble metal nanostructures: Quantum surface effects in crystallographic facets

Noble metal nanostructures are ubiquitous elements in nano-optics, supporting plasmon modes that can focus light down to length scales commensurate with nonlocal effects associated with quantum confinement and spatial dispersion in the underlying electron gas. Nonlocal effects are naturally more prominent for crystalline noble metals, which potentially offer lower intrinsic loss than their amorphous counterparts, and with particular crystal facets giving rise to distinct electronic surface states. Here, we employ a quantum-mechanical model to describe nonclassical effects impacting the optical response of crystalline noble-metal films and demonstrate that these can be well-captured using a set of surface-response functions known as Feibelman $d$-parameters. In particular, we characterize the $d$-parameters associated with the (111) and (100) crystal facets of gold, silver, and copper, emphasizing the importance of surface effects arising due to electron wave function spill-out and the surface-projected band gap emerging from atomic-layer corrugation. We then show that the extracted $d$-parameters can be straightforwardly applied to describe the optical response of various nanoscale metal morphologies of interest, including metallic ultra-thin films, graphene-metal heterostructures hosting extremely confined acoustic graphene plasmons, and crystallographic faceted metallic nanoparticles supporting localized surface plasmons. The tabulated $d$-parameters reported here can circumvent computationally expensive first-principles atomistic simulations to describe microscopic nonlocal effects in the optical response of mesoscopic crystalline metal surfaces, which are becoming widely available with increasing control over morphology down to atomic length scales for state-of-the-art experiments in nano-optics.


Introduction
Metals support collective oscillations of their conduction electrons, known as plasmons, with light-trapping and manipulation capabilities at nanometer length scales, i.e., well below the diffraction limit imposed by traditional op-tics. 1,2 The wealth of fundamental explorations in plasmonics over the last couple of decades has contributed to shape the field of nano-optics, 3 holding great promises for light-based technologies including theranostics, 4,5 photocatalysis, 6,7 plasmonic colors, 8 solar energy harvesting, 9,10 and quantum information. [11][12][13] Advances in modern nanofabrication techniques have driven plasmonics research further by enabling the realization of nano-optical devices that operate on deep-subwavelength scales. 14,15 As current capabilities can routinely pattern metallic nanostructures down to the few-nanometer regime, [16][17][18] where the frontiers of quantum and classical physics coalesce, new routes towards next-generation plasmon-based technologies begin to emerge, while also posing new challenges in understanding and modeling their optical response at nanometric scales. 13,19 The realization of thin crystalline noble-metal films [20][21][22][23] is key to cutting-edge explorations of novel plasmonic devices: metallic nanostructures with a high degree of crystallinity are anticipated to suffer reduced Ohmic losses when compared to their amorphous counterparts, 24 with the recent observation of plasmons in fewatom-thick crystalline silver films partially confirming this intuitive concept. 18 Furthermore, it is well-established in surface science 25 that (111) noble metal surfaces possess Shockley surface states, 26 whose features resemble those of a two-dimensional electron gas (2DEG). Shockley surface states can support 2D-like plasmon modes 27 that can be characterized, e.g., by angle-resolved spectroscopy, [28][29][30] while they can potentially play a role in near-field light-matter interactions at such surfaces.
First-principles simulation methods capture nonclassical effects in the optical response of ultra-thin metal films or few-atom metal clusters, 31,32 but necessitate intensive computational effort that rapidly becomes unfeasible for structures with characteristic sizes 10 nm; unfortunately, precision within ∼ 10 is what is currently afforded by state-of-the-art top-down nanofabrication techniques. One of the overarching challenges in theoretical nano-optics is thus to describe the optical properties of nanostructured metals by solving Maxwell's equations while accounting-in the response functions entering the constitutive relationsfor quantum-mechanical effects that emerge when electrons are confined in low-dimensional systems, ideally without resorting to overlydemanding numerical approaches that obscure the underlying physics. Fortunately, the situation is somewhat simplified in metals by their ability to effectively screen electromagnetic fields, which leads to an optical response dominated by surface effects. In this context, the concept of microscopic surface-response functions, like the Feibelman d-parameters, [33][34][35][36][37][38] offers a practical and scalable recipe to simultaneously incorporate quantum mechanical phenomena, such as electronic spill-out, nonlocal effects, and surface-enabled Landau damping, into the optical response of metal nanostructures, [33][34][35][36]39 as has been recently demonstrated experimentally. 40 Here, we apply a quantum-mechanical framework that describes the metal as a vertical stack of (homogeneous) atomic layers which are modeled as a realistic one-dimensional (1D) potential. The wave functions obtained by solving the corresponding Schrödinger equation are then used to compute the nonlocal optical response of selected noble metal (gold, silver, and copper) films with specific crystallographic orientations-namely, the (100) and (111) facets-from which we extract the associated Feibelman d-parameters. We demonstrate that the d-parameters obtained for a thick film (tantamount to a semi-infinite metal), once tabulated (as we do here), can be straightforwardly incorporated into a wide range of electromagnetic problems, ranging from analytical solutions for simple geometries to full-wave numerical electromagnetic solvers of realistic particles, to accurately describe intrinsic quantum mechanical effects affecting the system's optical response. 35,40 We anticipate that the results presented herein can be widely deployed to describe ongoing experiments and engineer future nanoscale plasmonic devices.

Results and discussion
Classically, the metal response to a monochromatic electromagnetic field of frequency ω can be described by the Drude model through its local-response dielectric function 41 where b (ω) is an ad hoc correction that accounts for the screening associated with the core electrons in the metal, while the second term describes the response of free electrons characterized by a plasma frequency ω p and a phenomenological scattering rate γ exp typically extracted from experiments. To maintain fidelity with experimental data, we follow ref 41 and construct b (ω) by subtracting the free-electron component from the experimentally tabulated dielectric function, 42 exp m , so that b (ω) = exp m (ω)+(ω exp p ) 2 /(ω 2 +iωγ exp ); the parameters used to characterize the noble metals that will be considered in this work, namely silver (Ag), gold (Au), and copper (Cu), are specified in Table 1 below. In classical electrodynamics, a metal surface is commonly described by a dielectric function that changes abruptly from the bulk, local dielectric function of the metal, m ≡ m (ω), to that of the adjacent dielectric, d ≡ d (ω). However, this naive prescription can be augmented with d-parameter-corrected boundary conditions that incorporate quantum effects associated with the optical response of metal surfaces. Specifically, for a p-polarized electromagnetic field impinging on a metal surface from the dielectric side, the nonretarded reflection and transmission coefficients read [34][35][36][37] respectively. Here, Q is the in-plane wave vector (i.e., parallel to the interface), while d ⊥ and d denote the frequency-dependent, complexvalued quantum surface-response functions introduced by Feibelman [33][34][35][36] (see Methods). On the other hand, for p-polarized light impinging on the interface from the metal side, the corresponding nonretarded reflection and transmission coefficients read Note that, in general, r md = −r dm when d α = 0, where α =⊥, . Equipped with eqs 2 and 3, the overall Fabry-Pérot (FP)-like reflection coefficient of the composite dielectric-metaldielectric heterostructure can be determined via where L denotes the metal film thickness. The above quantum-surface-corrected description only requires as inputs the Feibelman response functions, d ⊥ and d , which can be obtained for a particular dielectric-metal interface either experimentally 40 or from sophisticated quantum mechanical methods [e.g., timedependent density-functional theory (TDDFT) or empirical quantum-corrected models]. Here we obtain the d -parameters associated with crystalline noble metal surfaces from nonclassical optical response simulations of their reflection coefficients based on the random-phase approximation (RPA), employing the QM model reported in ref 43 and detailed in the Methods section. More specifically, we construct the RPA susceptibility from single-particle wave functions Ψ i (r) satisfying the Schrödinger equation for a one-dimensional (1D) phenomeno- logical potential V (z) characterizing each material; 44 this so-called atomic layer potential (ALP) is fitted to reproduce salient features of the bulk and semi-infinite surface electronic structures, so that the ALP-RPA description captures the effects of electronic band structure, electron spill-out, binding energies, and transverse atomic corrugation in the optical response of layered Ag, Au, and Cu films with (100) or (111) crystallographic orientation (see Table 3). See the Methods section for a comprehensive description of the calculation, in which the metal background polarizability b incorporating core electron screening is also included when describing the electron-electron interaction.
Using the ALP-RPA model, we follow the prescription outlined in the Methods section to extract the Feibelman d-parameters presented in Figure 1. The d ⊥ associated with crystalline noble metal surfaces are contrasted with those obtained within the specular reflection model (SRM), first proposed by Ritchie and Marusak to study nonlocal effects on surface plasmon dispersion, 45 and later generalized by others to deal with more complex structures; 36,41,[45][46][47] the SRM-also known as the semiclassical infinite barrier model (SCIB)-incorporates bulk spatial dispersion (i.e., nonlocality) but assumes a homogeneous electron gas and thus neglects atomic corrugations. In the case of the (100) orientation and the SRM, the absence of surface currents for charge-neutral materials fixes d = 0, 33 while d is introduced heuristically in order to incorporate the response of the Shockley surface states in the (111) facet; as explained in Methods, when extracting d ⊥ for the (111) facets we explicitly omit intraband transitions involving surface states to avoid double-counting the effect of the 2DEG. The Feibelman d ⊥ -parameters presented in Figure 1 clearly distinguish the surface response functions of a metal's crystallographic facets. In particular, there are two main regions at low energies ( ω < 1 eV) and at energies around ω cl sp = ω p / b (ω cl sp ) + d where the classical nonretarded surface plasmon is spectrally centered.
Nonretarded surface plasmon dispersion. Equipped with the d-parameters of various noble metals and their crystal facets, along with the analytical expression of eq 2, we illustrate their ability to reproduce the nonretarded surface plasmon dispersion (given by the poles of r dm ) obtained directly from the ALP. While the ALP method actually describes a crystalline metal film of finite thickness, we obtain well-converged results for N 50 monolayers,. Figure 2a shows the loss function Im{r} for silver with a (111) crystallographic orientation in the ALP model as a function of optical in-plane wave vector and energy; peaks in Im{r} indicate the surface plasmon dispersion, which tends toward zero frequency for small Q, in agreement with eq 4 for a film of finite thickness and non-vanishing e −2QL . Figures 2b-d display the nonretarded surface plasmon dispersion for Ag, Au, and Cu, obtained from the ALP (colored dots) along with the Feibelman d ⊥ -parameter (solid curves), together with available experimental data. [48][49][50] Here, the d ⊥ -parameter-results have been extracted by comparing eq 2 with the ALP reflection coefficient in the thick-film limit; this procedure, however, needs to be carried out judiciously as the conditions QL 1 and Q|d ⊥ | 1 must be simultaneously fulfilled. Chiefly, our results show that the optical response of a metal surface is determined by the surface's specific crystallographic orientation and that its surface response can be well-described in terms of the Feibelman d ⊥ -parameter, as evidenced by the overlapping dispersion relation calculated via the ALP model.
From eq 2, it follows that the nonclassical surface plasmon dispersion in the nonretarded regime-keeping only terms up to first-order in qd α -exhibits an approximately linear behavior with in-plane wave vector Q, namely [35][36][37] Re where ω cl sp ≡ Re ω p / b (ω cl sp ) + d is the classical nonretarded surface plasmon frequency. 47 While the surface-corrected Ag(111) response obtained from both ALP and SRM models is in good agreement with experiment, determination of the gold surface plasmon dispersion is complicated by an immediate onset of broadening in the loss function at low Q; the situation is further compounded for copper, where no clear maximum emerges in neither the response described by eq 2 with d -parameters nor in the direct ALP calculation, and only the SRM exhibiting well-defined maxima.
Acoustic surface plasmons due to Shockley surface states. At low frequencies, a feature exhibiting a nearly linear dispersion emerges in the loss function associated with (111)-faceted metal surfaces, indicating the existence of acoustic surface plasmons formed by Shockley surface states. 29,47,[51][52][53][54] Figure 3 shows the loss function of a Au(111) surface, which in the low-frequency regime is marked by the presence of a well-defined but relatively broad feature associated with its acoustic surface plasmon (c.f. scales of Figure 3a and Figure 2a). Next, we present the dispersion relation of acoustic surface plasmons akin to the Au(111), Ag(111), and Cu(111) surfaces obtained within the ALP framework (solid curves), and whose slope-corresponding to the acoustic surface plasmon velocity-is then determined through a linear fit (dashed curves); see Table 2). Note that the intrinsic acoustic surface plasmons supported by noble metal surfaces of well-defined crystallographic orientation have been characterized experimentally under different conditions: For Au(111), a phase velocity v φ /v 2D = 1.7 was observed at room temperature, 55 while a value v φ /v 2D ≈ 0.8 was reported at 78 K; 56 the extracted value in the present work (see Table 2) is close to unity. Ag(111) layers (i.e., in the semi-infinite limit) simulated with the ALP model. In panels (b-d) we present the dispersion relations of the indicated metals and facets as determined from maxima in their computed loss functions; we extract the maxima directly from Im{r} for ALP calculations (color-coded points), while the dispersion of the d-parameter corrected model is given by the poles of the denominator of r dm in eq 2 (solid curves). Triangular symbols represent experimental data (for Ag, 48 Au, 49 and Cu 50 ).  The Shockley 2DEG supported by the (111)facets can be accounted for through the Feibelman d -parameter (since it can mimic a surface conductivity). We exploit this by introducing, in an ad hoc fashion, a heuristic expression for d as described in the Methods section. We emphasize here that, because we account for the 2DEG heuristically, intraband transitions involving surface states are omitted in the ALPbased d ⊥ calculations. Then, in possession of both d ⊥ and d , we use eq 2 to reproduce the optical response calculations performed with the ALP model. As observed in Figure 5c-d, where the real and imaginary parts of the reflection coefficient for Au(111) are compared for different values of Q, the reconstruction is satisfactory for small Q, although the amplitudes of these already weak features are not well-reproduced.
Nonclassical optical response of ultrathin metal films. The practical utility of the d-parameter framework for mesoscale electromagnetism becomes apparent by recognizing that, once obtained for a specific dielectricmetal interface, they can be readily incorporated in a broad range of optical response calculations, either via d-parameter-corrected scattering coefficients 35 or through d-parametermodified boundary conditions. 39,40 As a concrete example, we now investigate the nonclassical optical response of ultrathin silver films comprising N (111) atomic monolayers. In Figure 4a-b we present the loss function for ultrathin silver films with thicknesses N = 5 and N = 20, which is dominated by the surface plasmon supported by the films; incidentally, the plasmon dispersion relation obtained from the ALP model closely resembles that obtained in a simple Fabry-Pérot description, even without including d-parameters. vectors predicted in the ALP model with FP models that include or neglect the d-parameter correction. Here, we observe an unremarkable effect resulting from d ⊥ for small values of the parallel wave vector (Q = 0.1 nm −1 ), while at larger in-plane momenta, for which the plasmon resonance approaches the plasma frequency (ω → ω p ), noticeable differences emerge (e.g., for Q = 0.8 nm −1 ). The excellent agreement between the calculated curves based on the d-parameter and ALP frameworks underscores how the optical response obtained analytically using eq (4) together with the Feibelman d-parameters (c.f. eqs 2-3) accurately accounts for quantum effects impacting the film's electromagnetic response. In particular, the optical response for extremely-thin Ag films, down to N = 5 atomic planes, appears to be wellreproduced by the surface-corrected thin film reflection coefficient, although the application of the d-parameters to describe such films is questionable.
Graphene next to crystallographically faceted metal films: Acoustic graphene plasmons We consider the extrinsic acoustic plasmons produced by the hybridization of a closely-spaced graphene layer with a crystalline metal film. Unlike the intrinsic acoustic plasmons supported by the (111)-facets, the introduced graphene layer and its opto-electronic tunability provides an additional knob to actively modulate the optical response of the emerging low-energy acoustic plasmon modes with linearized dispersion. 58 The experimental capability to position graphene within ∼ 1 nm of a noble metal layer is launching explorations of extreme light concentration within the gap region, 59-61 which could be further improved by employing crystalline noble metals. 43 In what follows, we summarize the semi-analytical FP description of the optical response based on the extracted d-parameters.
For a zero-thickness 2D graphene monolayer, the reflection and transmission coefficients in the quasistatic limit 62 read where σ(Q, ω) is the nonlocal conductivity of graphene, which we treat here at the level of the nonlocal RPA 36,63-65 (using Mermin's prescription for the relaxation-time approximation which conserves local particle number; 66 we take τ = 500 fs).
Similarly to eq (4), we compute the reflection coefficient of an extended graphene sheet on top of a semi-infinite metal via the FP resonance model as where s is the spacing between the metal surface and graphene, and r g and t g are the reflection and transmission coefficients of graphene. In our calculations we follow the prescription of ref 43 to account for spatial dependence of the carbon 2p orbitals ϕ 2p (r) extending outwards from the graphene monolayer plane, leading to the corrected graphene reflection and transmission coefficients r g = r 2D g C 2 g e −Qdg and t g = t 2D g C 2 g e −Qdg , where d g = 0.33 nm is the interlayer spacing of graphite and C g is a coupling factor defined in ref 43. Taking into account the aforementioned effective graphene thickness, the separation distance s actually corresponds to the distance between the edge of the graphene and the metal surface, i.e., s = 0 corresponds to a distance d g /2 between the graphene center and the metal surface.
Acoustic plasmons are anticipated to emerge even by depositing graphene directly on metallic films, 58 heralded by the prominent lowenergy linear dispersion feature in the reflection coefficient of Figure 5a. For the considered hy-brid graphene-Au(111) surface, the excited extrinsic acoustic plasmons are characterized as before by ω = v gr Q with v gr denoting the associated group velocity; in such a heterostructure, v g is determined by graphene Fermi energy E F and the graphene-metal spacing s, as illustrated in Figure 5b. Accordingly, a given dispersion velocity can be obtained by different combinations of E F and s. Figure 5c reveals that among Ag, Au, and Cu (not shown), neither the choice of metal nor the considered crystalline facet strongly influences the acoustic plasmon dispersion characteristics; at low energies, these metals are all good conductors that effectively screen the graphene plasmon and render its dispersion acoustic. We remark, however, that, for the same heterostructure, the dispersion relation of the higher energy plasmon mode is indeed dominated by the metal properties. 43 However, inspection of the linewidths of the acoustic plasmon in Figure. 5d reveals a substantial dependence on the quantum-mechanical effects aris- ing from the various crystalline facets, which is underlined by the underestimation of spectral widths in the SRM. 43 In particular, the obtained results suggest that crystalline Au(111) gives rise to additional surface-enhanced damping when compared to Au(100), presumably due to the presence of a surface state, and warranting further study of the acoustic plasmons supported by heterostructured crystalline metal films.
Crystallographically faceted nanoparticles. Going beyond planar, layered media, we explore the role of crystallographic orientation in faceted noble metal nanoparticles (NPs). The optical response of metallic NPs is dominated by the localized surface plasmon (LSP) resonances supported by it; the most prominent of which are typically those of dipolar character, as they can couple to far-field radiation. As the NP size is reduced towards nanometric dimensions, the ensuing NP's surface-to-volume ratio grows and leads to more pronounced nonclassical corrections associated with the NP's quantum surface-response. To illustrate the importance of using the appropriate Feibelman d-parameters for determining the quantum surface-response associated with specific facets, we consider in Figure 6 a realistic faceted silver NP. Any natural NP-especially those with characteristic dimensions 10 -20 nmno matter how carefully synthesized, will always deviate from a perfect sphere, as a consequence of its growth in a sequence of specific crystallographic planes. 67 The shape closest to a sphere is that of a truncated octahedron, as depicted schematically in Figure 6; it is characterized by large hexagonal (111) surfaces and smaller (100) facets.
To compute the NP's nonclassical extinction cross-section σ ext , we implement in a finiteelement method (FEM) solver the mesoscopic boundary conditions within the Feibelman formalism; 35,36,39,40 in practice, this is tantamount to the introduction of surface electric and magnetic currents (see Methods) characterized by the d-parameters presented in Figure 1. Comparing with the classical spectra (red curves), it is clear that the effect of the d-parameters is to capture the nonlocal optical response of such an NP. As expected for silver, the spectra are shifted to higher energies as a result of an inward shift (Re{d ⊥ } < 0) of the screening charges. Studying different NP sizes, from 3 to 7 nm in "radius" (meaning here the circumscribed radius R circ ), a consistent trend is observed, with the resonance broadening (as a result of increased Landau damping) and under- Figure 6: Extinction spectra for a Ag truncated octahedra including quantum surfacecorrections. Optical extinction cross section (normalized to the radius of the circumscribed sphere R circ for Ag truncated octahedra with R circ = 3, 5, and 7 nm (a, b, and c respectively), as shown in the schematics on top. Red curves correspond to the response of bulk Ag (classical with no d parameters), and black curves to facets (100) and (111) described by their corresponding d parameters shown in Figure 1a and b. The blue and yellow curves in the inset of (a) show the corresponding spectra if the entire particle is described entirely by the Feibelman parameters of (100) and (111) facets, respectively. going stronger blueshifts, as the size decreases, which is compatible with the predictions of nonlocal hydrodynamics, 41,68,69 and also in accordance with electron energy-loss spectroscopy (EELS) experiments. [70][71][72][73] The observed behavior is mainly attributed to the (111) facets, as shown in the inset of Figure 6a, where blue and yellow curves show the corresponding spectra assuming that the entire NP is described solely by the (100) or (111) d-parameters, respectively; since the associated surface area of the (100) facets is smaller, and their corresponding Feibelman parameters are significantly lower in magnitude than those of the (111) surface, they only induce a small frequency shift in the spec-tra. Then, while the main effect is due to the (111) facets, the corresponding spectrum almost coincides with the "mixed" one, where each facet is described by its own parameters. It is also worth noting that, because the truncated octahedron constitutes a highly symmetric shape, the optical response of such NPs resembles that of spheres, and thus changing the angle of incidence is not expected to lead to significant differences.

Conclusions
The inherently large losses sustained by noble metals is often regarded as the Achilles heel of nano-optical functionalities based on plasmonics, and has motivated intensive efforts to identify new material platforms that can support long-lived polaritons. Crystalline noble metal films constitute one appealing possibility that is now becoming increasingly available. We have introduced surface response corrections in the form of Feibelman d-parameters extracted from quantum-mechanical optical response calculations of crystalline noble metals. As demonstrated here, the tabulated dparameters for gold, silver, and copper surfaces with specific crystallographic orientations can be straightforwardly incorporated in analytical models as well as in computational electromagnetic solvers for the calculation the nonclassical optical response of various nanoplasmonic systems of current interest that contain these facets. We envision that the d-parameters reported here can be widely deployed to describe quantum surface effects in crystalline noble metal surfaces which are actively explored for novel nanophotonic functionalities and applications.

Microscopic surface-response functions: Feibelman d -parameters
The surface-response functions introduced by Feibelman, 33 d ⊥ and d , are formally given, respectively, by the first-moment of the quantum mechanical induced charge (ρ ind ) and of the parallel component of the corresponding current density (J ind x ) 33-36 in the long-wavelength limit. The Feibelman parameters can be rigorously incorporated in electrodynamic problems by appropriately adjusting the boundary conditions at the surface, 33 both in analytical treatments [34][35][36]39 and in numerical implementations. 39,40 It is implicit that in the presence of more mechanisms, d ⊥ j d (j) ⊥ and d j d (j) . In the following we separately consider contributions from bulk spatial dispersion-relevant for any metal surface with a compressible electron gas-and Shockley surface states-relevant to the (111) noble-metal surfaces.

Contributions from bulk spatial dispersion
The surface-response functions associated with the spatial dispersion (nonlocal response) of the bulk response functions of the metal, i.e., the wave vector dependence of the dielectric function, can be expressed through the specular reflection model 36,46,47 Note how the vanishing of d is a consequence of the charge-neutral interface. 37 Below, we discuss how this is changed in the presence of a Shockley surface state.

Contributions from Shockley surface states
We consider the effect of a Shockley surface state, while for simplicity leaving out the response associated with the spatial dispersion of the bulk states. We note that, with our sign convention of the surface-normal, the surface conductivity σ 2D is related to the Feibelman d parameter as σ 2D (ω) = −iω 0 ( m − d )d (ω). 36 Furthermore, in the nonretarded limit, the 2D plasmon dispersion relation associated with σ 2D is generally given by q(ω) = iω 0 ( m + d )/σ 2D (ω). 36 Anticipating that we have a Shockley surface state that supports an acoustic plasmon, 27 corresponding to a dispersion relation ω(ω + iγ 2D ) = v φ q, we thus find the connection between this phase velocity v φ and the Feibelman parameters to be: For low frequencies (ω ω p ) this simplifies to d (ω) v φ / ω(ω + iγ 2D ), while it naturally vanishes at high frequencies. When implemented in the mesoscopic boundary conditions for the electrodynamics, these Feibelman parameters-by construction-support an acoustic plasmon with the desired phase velocity. Considering the poles of the scattering coefficients (see eqs 2 and 3) and substituting in eqs 9, we obtain and we indeed find two decoupled solutions: the "classical" surface plasmon resonance (defined by m + d = 0) and the (added ad-hoc) acoustic one with ω(ω + iγ 2D ) = v φ q. Table 2 gathers the obtained phase velocities by fitting the acoustic surface plasmon featured in the optical response (see Figure 5a), and the damping γ 2D by associating it to the width of a Lorentzian; specifically, the widths have been computed numerically from the second derivative of the imaginary part of the reflection coefficient in order to remove the background contribution (see Figure 5d).

Extraction of the d-parameters for crystalline metal surfaces
We describe crystalline metal films quantum mechanically, simulating their optical response at the level of the random-phase approximation (RPA) following the procedure of ref 43: As explained therein, metal films are considered to have translational symmetry in the R = (x, y) plane, so that their electronic wave functions are amenable to expansion in a plane wave basis according to Ψ j,k (r) = A −1/2 e ik ·R ϕ j (z), with A denoting the normalization area, k the 2D electron momentum, and ϕ j (z) the spatial dependence of state j in the quantization direction z; the latter quantity is obtained by solving the eigenvalue problem Hϕ j (z) = ε ⊥ j ϕ j (z) to obtain the associated energy eigenvalues ε ⊥ j of the 1D Hamiltonian H = − 2 ∂ 2 z /2m e + V (z) determining the band dispersion ε j,k = 2 k 2 /2m e + ε ⊥ j . The 1D potential V (z), here referred to as the atomic layer potential (ALP), is selected from those reported in ref 44 that characterize faceted metals of thickness t composed of N atomic planes stacked in the z-direction with inter-layer spacing a s (naturally, t is an integer multiple of a s ).
Electronic bands are populated by successively filling the lowest bands until the effective bulk electronic density n eff is reached, thereby determining the Fermi energy E F as where the sums over j terminate when ε ⊥ M < E F / < ε ⊥ M +1 , i.e., j = M is the highest partiallyoccupied band. The electronic densities n eff are determined by imposing the experimentallyestablished value of E F for a given noble metal in the bulk limit (i.e., for a sufficiently thick film). For consistency with experimental observations, 74-76 we impose a linear variation in the effective mass of the parabolic bands as a function of their quantized energies ε ⊥ j according to m * j /m e = a ε ⊥ j + b, thereby avoiding artifacts due to an unrealistic number of excitation channels for vertical transitions introduced by perfectly-aligned parabolic bands; specific parameters used in our calculations are reported in Table 3-note that the surface states for (111) noble metal facets are assigned specific experimentally-determined effective masses. We characterize the optical response of noble metal films by the reflection coefficient R(Q, ω), expressed as a function of the optical in-plane wave vector Q and frequency ω. Considering that the relevant length scales are far smaller than the involved optical wavelengths, we invoke the quasistatic approximation to compute the reflection coefficient in terms of electrostatic potentials as R = 1 − φ(z)/φ ext (z), where φ = φ ext + φ ind is the sum of external and induced potentials, the former exciting the system and the latter obtained as φ ind (z) = dz χ(z, z )φ ext (z ) in terms of the film susceptibility χ, which we compute in the ALP-RPA formalism reported in ref 43. In principle, the RPA response function is constructed by summing over all possible transitions between electronic states; however, because the Shockley surface state of the (111) facet is incorporated in d following the ad-hoc prescription of the previous section, the reflection coefficient that is used to extract d ⊥ is computed by excluding intraband transitions involving surface states, thereby avoiding double-counting of the associated 2DEG.
Following the RPA description of ref 43, we correct the Coulomb interaction to incorporate screening from core electrons using the experimentally-extracted polycrystalline dielectric functions b (ω) parametrized in Table 1 and plotted below in Fig. 7. The d ⊥ associated with a given noble metal facet is extracted from ALP-RPA simulations of a sufficiently thick film, so that the optical response is converged with the number of atomic planes. More specifically, we obtain d ⊥ by fitting eq 2 to the ALP-RPA reflection coefficient of the thick film, employing the corresponding bulk dielectric function m of eq 1; importantly, the bulk plasma frequency ω p for each metal facet that enters eq 1 in the fitting is obtained from the ALP-RPA response of a finite film for sufficiently small in-plane wave vector, e.g., Q ∼ 0.005, so that nonlocal effects are safely neglected and the surface plasmon resonance (∼ 1 eV) is captured in an uncorrected Fabry-Pérot description; this procedure enables a stable parametrization of a crystallographic facet's bulk properties when the contribution from the d-parameters is negligible. In practice, the surface plasmon for a finite film with N = 10 − 40 atomic planes appears at lower energies than the surface plasmon for the semi-infinite film, and the associated resonance is undamped by interband transitions, thereby giving rise to a well-defined peak (see, for instance, Figure 4) from which ω p is obtained by fitting m = d (1 + r dm )/(1 − r dm ), c.f. eq 2 in the Q → 0 limit.   Table 3 for the characteristic parameters of each metal.
Once the bulk properties for each facet are set, we construct d ⊥ (ω) at a given ω by fitting eq 4 for a given value of Q to the ALP-RPA response. We maintain a large number of layers to avoid quantum finite size effects emerging in thin films (< 20 atomic layers). We then confirm that convergence is maintained with the calculated d ⊥ (ω) as the number of layers and/or parallel wave vector is varied, where the latter condition is typically satisfied for 0.3 < Q < 1.5 nm −1 . After the parameters are obtained, they are applied to thick ( Figure 2) and thin films (Figure 4) to confirm their applicability.
In the case of the (100)-facet metals we set d = 0, whereas for the (111)-facet the prescription of eq 9 is employed to describe the intrinsic low-energy acoustic plasmon; as anticipated, d ⊥ (ω) approaches zero when ω → 0, and hence such a condition was imposed after solving for d ⊥ (ω). Table 3: Characterization of quantum well states in noble metals. The parameters defining the electronic bands of noble metals entering our optical response calculations are presented. The quantities a and b define the linear variation in effective mass for band j as a function of its associated energy ε ⊥ j ; values for the Fermi energy E F , the effective mass associated with surface states m * (SS), and the effective mass for the bottom of the conduction band m 0 , are extracted from experimental reports, while the effective electron density n eff is fitted to match E F , thus fixing the plasma frequency ω p .

Finite-element implementation
To calculate the extinction spectra of the truncated octahedra, we used the commercial FEM solver Comsol Multiphysics 5.4. As it has been shown elsewhere, 35,39,40 Feibelman parameters can be incorporated in any computational method by adjusting the boundary conditions. More specifically, the Feibelman parameters introduce discontinuities in the parallel components of the electric and magnetic fields, which can be expressed through 35,39,40 n × (E 2 − E 1 ) = −d ⊥n × ∇ n · (E 2 − E 1 ) , n × (H 2 − H 1 ) = iωd n × ∇ n · (D 2 − D 1 ) ×n , where E i , H i , and D i are the electric, magnetic, and displacement fields on side i of an interface between two media 1 and 2, andn is the unit vector normal to the interface. These conditions can be readily implemented in version 5.4 of Comsol Multiphysics, through surface current and surface magnetic current densities, expressed through the down and up functions in Comsol Multiphysics for the fields at sides 1 and 2. Since the expressions for these currents (right-hand sides of eqs 12a and 12b) contain the fields themselves, the problem needs to be solved iteratively, starting with the currents due to the incident plane wave. To calculate scattering and absorption cross section, we need to integrate the Poynting flux of the scattered and total field over a surface (a sphere) enclosing the NP, with radius large enough (2 − 3 nm more than R circ ) to ensure that numerical noise due to the currents close to the surface will be minimum. For sharp-edged NPs like the octahedra studied here, it is also necessary to introduce some rounding, to ensure that any spurious edge/corner modes will be absent. This is needed in the classical calculation, in the absence of surface currents; when these are present, the additional damping they introduce smooths things nicely. However, for a direct comparison between the two cases, it is necessary to include the same rounding in both cases. This, however, causes an additional numerical problem, because the iterative method diverges when surface currents are added in such small rounded elements. For this reason, surface currents are only used to describe the square (100) and hexagonal (111) facets. This is in practice not a bad approximation, as one needs to somehow introduce a smooth transition between the two different current densities. In terms of set-up parameters, a cubic physical domain of side 300 nm was used, surrounded by 300 nm-thick perfectly-matched layers. For the finite-element discretization, a mesh of 30000 domain elements with maximum element size 20 nm and minimum element size 0.5 nm provided converged spectra.