Exploring charge-transfer effects at metal–molecule interfaces through modeling surface-enhanced Raman spectroscopy (SERS)

Abstract

Understanding charge transfer (CT) at metal–molecule interfaces is central to the design of new materials for catalysis, sensing, and energy applications. Surface-enhanced Raman spectroscopy (SERS) provides a powerful technique to probe metal–molecule interfaces since the SERS spectra reflect changes in geometry, orientation, and electronic structure which are influenced by CT. While CT can either be intrinsic or excitation-driven, we will focus on how the former contributes to SERS. However, it remains difficult to model these effects since it requires electronic structure methods capable of treating large systems. In this study, we present an efficient model to study these effects by combining a simplified time-dependent density functional theory (TDDFT) approach with a Raman bond model. The first-principles Raman bond model partitions Raman intensities into bond contributions such that the chemical mechanism in SERS can be interpreted as interatomic charge-flow modulations. Using this new model, we will examine how molecular orientation and intermolecular interactions affect the interfacial CT, employing N-heterocyclic carbenes (NHCs) as a model system. Using the Raman bond model, we will characterize the degree to which the interfacial charge flow is influenced by these effects. We then apply this model to characterize how charge flow affects molecular imaging using tip-enhanced Raman scattering (TERS). Using TERS, it is possible to image single molecules with sub-nanometer resolution. However, the role of CT in this process remains to be elucidated. We will show how the Raman bond model can characterize the importance of interfacial CT in TERS imaging. Overall, our results demonstrate that the Raman bond model, when combined with efficient first-principles calculations, offers a powerful approach for interpreting SERS spectra and providing new insights into interfacial CT.

---

1 Introduction

Understanding the binding and organization of molecules at interfaces is essential in many areas ranging from biology to materials science. Vibrational spectroscopy provides detailed information about the structure of molecules, how they organize at interfaces, and their interactions with the local environment. One approach to performing vibrational spectroscopy at metallic interfaces is based on surface-enhanced Raman spectroscopy (SERS).^1–4 In SERS, the Raman signal of molecules adsorbed on plasmonic metal substrates is enhanced by orders of magnitude compared to normal Raman scattering. The enhancement has been demonstrated to be large enough to enable single-molecule detection,^5–9 which illustrates the great potential of SERS for providing molecule-specific information even at very low adsorbate concentrations.

The SERS enhancement is typically attributed to two mechanisms, which in combination lead to the observed SERS spectra.^{1–4,10–14} The first mechanism is referred to as the electromagnetic mechanism (EM) and constitutes the majority of the signal enhancement.^12,15–18 This enhancement arises from the significantly increased and localized electric field at the metallic interface caused by the plasmon excitations in the substrate. In the simplest formulation of the EM, it can be shown that the enhancement scales as |E|⁴, where |E| is the enhancement in the electric field at the position of the molecule.^14,18,19 Since |E| can reach factors of 10–1000 (ref. 20) in junctions or crevasses of nanoparticle aggregates due to the plasmon excitation, it provides the simplest explanation for the large enhancement observed in SERS. The localization of the enhanced local field to small regions around the substrate is what gives SERS its surface selectivity and high sensitivity. Classical electrodynamics simulations are often sufficient for understanding EM,^21–24 although for very small separations of nanoparticles, QM effects become important because of short-range electronic tunneling, which leads to a reduced local field in a junction.^25–29 Another aspect of the large EM enhancement is that molecules or vibrational normal modes that are aligned with the local field are preferentially enhanced, which is referred to as the surface-selection rules.^18,30–32 This gives SERS selectivity towards molecules with a more upright adsorption geometry compared to molecules that lie flat, and this enables them to be detected even when they are the minority species on the surface.³³

The second mechanism arises from the specific interactions between the molecule and the substrate, and it is referred to as the chemical mechanism (CM). CM is often further divided based on whether it involves a resonance process with the incident light. The non-resonant mechanism is referred to as CHEM and describes the enhancement that can be attributed to the geometric and electronic changes in the molecule upon binding to the surface. These changes lead to shifts in the vibrational frequencies as well as the intensity of modes that are sensitive to the binding of the molecule. In contrast, the resonance mechanisms describe the enhancement arising from the incident light directly driving a charge-transfer transition between the molecule and the metal or a transition within the molecule. Since this is a resonance phenomenon, the enhancement is typically larger than that of CHEM. However, in contrast to CHEM, the resonance mechanisms also depend strongly on the incident light matching the electronic transition. Therefore, in this work, we will focus on CHEM, since it can provide insights that are general to molecules binding on metal surfaces.

Since CHEM depends strongly on the binding of the molecule to the substrate, it requires a quantum-mechanical treatment of the molecule–metal complex to correctly model it. For this reason, there has been a large focus on using electronic-structure methods to understand CHEM by modeling the Raman enhancement from molecules binding to small metal clusters.^33–37 Several models to interpret CHEM based on electronic-structure methods have been proposed that either consider a few important electronic transitions^38–41 or the charge flows near the molecule–metal interface by partitioning the induced electron density of the SERS system in real space.^42–44 An example of the real-space partitioning is the first-principles Raman bond model (RBM) that we have recently introduced to model the SERS enhancement.^42,43 In this model, the Raman intensity is partitioned into contributions from the induced dipole on each atom and the induced charge flow between atoms, so-called Raman bonds. This is done using a real-space partitioning of the induced charge density into atomic contributions based on a Hirshfeld analysis,⁴⁵ from which intrinsic and charge-fluctuation atomic polarizabilities can be defined.^46,47 The intrinsic atomic polarizability corresponds to the response of the atom to the external field, whereas the charge-fluctuation atomic polarizabilities reflect the charge redistributions due to the external field. To eliminate the origin dependence of the charge-fluctuation atomic polarizabilities, they are further re-expressed as charge flow between atoms using the LoProp method.⁴⁸ In this method, long-range charge flows are penalized, effectively assigning the charge flow to follow the bonds in the system. The RBM then partitions the Raman intensities based on the changes in the intrinsic polarizabilities and the induced charge flow between atoms due to a given molecular vibration. Thus, the RBM depends on both the Hirshfeld analysis and the reformulation in terms of charge flows, which are not uniquely determined. However, we have shown that stable charge-flow patterns can be achieved once long-range charge flows are eliminated.⁴² The RBM shows that CHEM correlates with the charge flow of electrons across the molecule–metal system modulated by the vibrations.^42,43 A major advantage of the first-principles Raman bond is that it is applicable whether the incident light is on- or off-resonance with transitions and thus provides insights into the interplay between EM and CM effects.^42,43 A common problem with using electronic-structure methods to model CHEM is the high computational cost of these methods. Therefore, it would be advantageous to combine methods to interpret SERS, such as the Raman bond model, with faster electronic-structure methods.

In this work, we combine the first-principles Raman bond model with a simplified time-dependent density functional theory (TDDFT) approach based on damped response theory to calculate and interpret SERS spectra. The simplified TDDFT approach was recently demonstrated to provide calculated SERS spectra in good agreement with results from full TDDFT but at a significantly reduced computational cost.⁴⁹ The combination of this computationally efficient method with the Raman bond model enables the understanding of CM in SERS by modeling large systems where both the EM and the CM are important. To illustrate this, we first demonstrate that the simplified Raman bond model (sRBM) produces an interpretation of SERS spectra consistent with RBM using full TDDFT. To this end, we will use the prototypical model system of pyridine interacting with an Ag₂₀ cluster. We will then use the sRBM to gain insights into the bonding of N-heterocyclic carbenes (NHCs) on gold nanoparticles. Specifically, we will demonstrate how CT effects play an important role in the SERS surface-selection rules. Finally, we will use the sRBM to gain insights into tip-enhanced Raman scattering (TERS) images of single molecules.

2 Methods and computational details

In this work, all calculations were performed using a local modified version of Amsterdam Modeling Suite (AMS).⁵⁰ In the TERS simulations, a triple-ζ polarized (TZP) basis set with a large frozen core from the ADF standard basis set library was used. For all ground-state properties, the Becke–Perdew (BP86) XC functional was used,^51,52 and relativistic effects are treated using the zeroth-order approximation (ZORA).⁵³ The geometry optimization was carried out for a free benzene molecule and a benzene molecule on a Ag₄₀ substrate. For the free benzene molecule, the analytical Hessian was used to obtain the vibrational normal modes.^54–56 For the benzene–Ag₄₀ complex, a mobile block Hessian was used,^57,58 which treats the entire Ag₄₀ substrate as a rigid body. In the SERS simulations of NHC, GFN1-xTB was used for both geometry optimization and frequency calculations.^59,60 For all NHC systems, the frequency calculations were also carried out using the mobile block Hessian, which treats the Au substrate, except for the adatoms, as a rigid body.

Based on the optimized geometry, the absolute Raman intensity is calculated as the differential cross-section under the short-time approximation, corresponding to an experimental setup with perpendicularly plane-polarized light and a 90° scattering angle for the detector:^61,62

where [small nu, Greek, tilde]_in and [small nu, Greek, tilde]_k are the incident frequency and the k-th vibrational frequency, and is the zz component of the polarizability derivative with respect to the k-th vibrational normal mode. Further is the differential cross-section, and h, c, and k_B are the Planck constant, speed of light, and Boltzmann constant, respectively. The polarizability derivatives were calculated based on numerical differentiation with respect to normal mode displacements,⁶² and only the zz component was used based on the surface selection rules.^18,30–32 For the TERS image simulations of the benzene molecule, the z-direction was chosen to be perpendicular to the molecule plane, while for the NHC simulations, the z-direction was chosen to be perpendicular to the Au substrate surface.

The static and complex dynamic polarizabilities used in eqn (1) were obtained using a simplified damped-response TDDFT implemented in the AOResponse module of a local modified version of the ADF engine.^49,63 In regular TDDFT, the central equation to solve for damped response is as follows:⁶³

where ω is the incident frequency and ω_ia = ε_i − ε_a is the difference in Kohn–Sham orbital energies, with V^ext being the external perturbation, and the δP_ia and δP_jb terms being elements of the perturbed density matrix. We use the convention that i, j, … are the indices for occupied orbitals and a, b, … are the indices for virtual orbitals. A phenomenological damping factor γ is used to represent the finite lifetime of the excited states. The most time-consuming part of solving eqn (2) is associated with the coupling matrix element K_ia,jb, which involves electron integrals in the formwhere f_XC is the exchange–correlation kernel, and the r and r′ terms are dummy variables for coordinates in real space. In this work, we took the DFTB approximation that was adopted in the so-called TDDFT+TB method for the coupling matrix such that^64,65where the transition charges q_ia,A and q_jb,B are evaluated using Löwdin population analysis and the γ_AB term describes the interactions between atomic transition charges q_ia,A and q_jb,B based on the interatomic distance |r_A − r_B|, as well as the chemical hardness of corresponding elements.⁶⁵ Although drastic simplifications have been made to the coupling matrix, previous work has shown that a calculated Raman cross-section in good agreement with regular TDDFT can be achieved.⁴⁹ Although TDDFT Raman intensities are typically semi-quantitative, the calculated SERS spectra are often of sufficient quality to allow reliable assignment when compared with experimental results.^66,67 However, we expect that, in addition to the shortcomings associated with using TDDFT, the neglect of anharmonic effects will influence the comparison with experimental results. We also expect that the neglect of the local environment and approximations in the adapted cluster model will also affect the comparison.

In this work, we apply a simplified Raman bond model to analyze SERS enhancement mechanisms. In the Raman bond model, the Hirshfeld charge analysis⁴⁵ is adopted to partition the total polarizability derivatives in eqn (1) into atomic and interatomic contributions:^42,43,68

where α,β ∈ {x,y,z} are indices of the Cartesian components, i and j are indices of the atoms, and Q_k is the k^th normal mode. The first term in eqn (5) is the atomic contribution where −r_β is the electric dipole moment operator in the β direction, R_i,β is the β component of the position of atom i and δρ_i,α is the perturbed electron density induced by an external field in the α direction, which is obtained using Hirshfeld partitioning.⁴⁵ The second term corresponds to the inter-atomic contribution where q_ij,α is the charge flow between atoms i and j induced by an external field in the α direction, which is calculated using the LoProp method⁴⁸ to ensure origin-independent results. The atomic contributions are termed as Raman atoms (R^atom) and the interatomic contributions are termed as Raman bonds (R^bond).⁴³

The overall Raman polarizability derivatives can be grouped into intra- and inter-fragment contributions utilizing the sRBM. These fragments can be either the molecule or the cluster. For fragments A and B, the intra-fragment contribution is

and the inter-fragment contribution is

To quantify the relative contributions of each term, a projection scheme was adopted such that

where R^total is the summation of all intra- and inter-fragment contributions, and p^intra and p^inter are intra- and inter-fragment Raman bonds, respectively. The original implementation of the Raman bond model was based on polarizability derivatives obtained using TDDFT.^42,43,68 In this work, we will instead use the polarizability derivatives obtained using the simplified damped-response approach in the analysis, and we will refer to this model as the simplified Raman bond model (sRBM).

3 Results and discussion

3.1 Raman bond model with simplified damped response

As shown in eqn (5), the Raman bond model is constructed based on partitioning of the response of the electron density and was originally based on results from regular TDDFT. However, it is unclear if the Hirshfeld partitioning scheme can generate qualitatively similar results when based on the simplified damped-response results. Specifically, one of the simplifications adopted in the simplified damped response involves approximations to the transition electron density. As shown in eqn (4), the transition charge q_ia,A can be regarded as a monopole approximation to the orbital product that corresponds to the transition density present in eqn (3). It has been shown that such an approximation leads to more localized electron density distribution; as a consequence, the weaker screening leads to an overestimation of the polarizability.⁴⁹ Although such overestimations can be largely reduced when the polarizability derivatives are explicitly calculated, it is unclear whether a similar error cancellation will occur in this sRBM.

Therefore, to verify that the Raman bond model is still applicable using the simplified damped response, we performed a direct comparison with the decomposition obtained using regular TDDFT. As a proof-of-principle, we constructed the Raman bond model based on the model system of py–Ag₂₀ used in the previous work on resonance with the plasmon-like excitation of the Ag₂₀ cluster.⁴³ Previous work showed that the simplified damped response can qualitatively reproduce the SERS spectra of this system using regular TDDFT.⁴⁹ The Raman bonds were decomposed into three components, namely cluster, molecule and inter-fragment contributions, and the comparison between the simplified damped-response theory and regular TDDFT is shown in Fig. 1.

Fig. 1 Comparison of Raman bond contributions for Ag₂₀–Pyr using different methods: (a) simplified damped response, and (b) regular TDDFT. The incident frequency is chosen as the plasmon-like resonance of Ag₂₀ (3.45 eV).

Comparing the top row of Fig. 1(a) and (b), we see that the summation of Raman bonds and Raman atoms calculated using simplified damped response agrees with regular TDDFT with only small differences in the relative magnitudes of the peaks. For example, the peaks around 1200 cm⁻¹ and 1450 cm⁻¹ are overestimated in the simplified damped-response results compared to those in the regular TDDFT.

More importantly, we see that the Raman bond model based on the simplified damped response theory and regular TDDFT produce qualitatively similar decompositions. Both methods show a dominant contribution from p^cluster, which indicates the importance of the EM at this frequency. This agrees with the assumption that when the incident frequency is on resonance with the plasmon, the EM should give the majority of the overall enhancement. Similar agreement is also achieved for p^inter, where the majority of the difference comes from the 1200 cm⁻¹ and 1450 cm⁻¹ regions. While for the p^mol contributions, both methods show only minimal magnitudes compared to the p^cluster contributions.

In summary, the projected Raman bond contributions obtained using the simplified damped response show good qualitative agreement with those obtained from regular TDDFT, although some discrepancies in the relative peak intensities remain. Specifically, the mechanisms of SERS enhancement can be captured using the sRBM. Thus, the Raman bond model can be extended to understand the interfacial CT of large systems when combined with the simplified damped response.

3.2 NHC on Au surfaces

As a first test of the sRBM to understand interfacial CT, we considered the binding of N-heterocyclic carbenes (NHCs) on gold nanoparticles. NHCs have emerged as alternatives to thiols, utilizing strong sigma bonds to form well-defined, highly ordered self-assembled monolayers on transition-metal surfaces and nanoparticles.^69–71 Depending on the specific NHC, they have been found to adapt orientations that are either flat or vertical on the surface.^72,73 Previous work using SERS has found that these orientations can co-exist on the surface of nanoparticles and that both configurations are necessary for explaining the observed spectra.^33,67 This was found to be true even when the flat orientation was the majority species, since the surface selection rules enables the detection of minor concentrations of the vertical orientation.^33,67 However, the role of the interactions between the molecules in determining the SERS spectra has not been examined. Therefore, we will use the sRBM to analyze the SERS spectra of the flat and vertical configurations of methyl-benzimidazolium N-heterocyclic carbene on a Au₅₈ nanocluster, as well as their interactions.

In Fig. 2, we plot the comparison between the sRBM decomposition for the SERS spectra of methyl-benzimidazolium N-heterocyclic carbene on a Au₅₈ nanocluster in the flat and vertical orientations. The SERS spectra have been calculated in the non-resonant limit and thus only includes contributions from the CHEM. We first begin with the Raman bond decomposition for the flat configuration, which can be seen in Fig. 2(a). The largest signals in the Raman spectrum are in the 1350–1450 cm⁻¹ and the 1000–1200 cm⁻¹ ranges, and at 1600 cm⁻¹ with additional strong bands around 650 and 740 cm⁻¹. From the decomposition, we see that most of the contributions come from p^cluster and p^inter, with only minimal contributions from p^mol. The fact that p^inter is large shows that the interfacial CT is providing a significant contributions to the Raman scattering. Interestingly, only two modes have significant contributions from p^mol. The mode around 740 cm⁻¹ corresponds to an out-of-plane vibration, which shows a large contribution from p^cluster. Therefore, for the flat orientation, the Raman signal originates from constructive interference between Raman bonds from the cluster and the inter-fragment contributions but with limited contributions from the Raman bonds in the molecule. This again is likely due to the flat orientation of the molecule where the vibrational modes that are predominantly aligned with the plane of the benzene ring cannot directly lead to large Raman bonds.

Fig. 2 sRBM decomposition for an NHC on an Au₅₈ surface in the (a) flat and (b) vertical configurations.

For comparison, in Fig. 2(b) we plot the sRBM decomposition for the NHC oriented vertically on the nanocluster. This spectrum is noticeably different from the flat configuration with the peaks in the region of 1350–1450 cm⁻¹ being the strongest. Also, we note that the magnitude of p^total is significantly greater than that of the flat configuration. This is in line with the expectations from the SERS surface selection rules which preferentially enhance molecules and vibrations that are aligned with the local field, which in this case is perpendicular to the metal surface. From the sRBM decomposition, we see that the Raman signal arises from constructive interference between p^mol, p^cluster, and p^inter. In contrast to the NHC in the flat orientation, we see that the strongest contribution for the vertical orientation is from the charge flow within the molecule. The other two contributions, p^cluster and p^inter, largely mirror the p^mol contributions leading to an overall constructive interference. Thus, the sRBM decomposition reveals the different character of the surface selection rules for molecules oriented parallel or perpendicular to the local fields. In the parallel case, the vibrations in the molecule lead to a coherent charge flow in the molecule and the cluster leading to a large Raman scattering. For the perpendicular case, the vibrations cannot effectively modulate the charge flow in the molecule and thus the largest contributions to the Raman scattering come from the interfacial charge flow.

As mentioned above, the flat and vertical configurations of NHCs can co-exist on nanoparticles. To understand how their interactions influence the SERS spectrum, we used the sRBM decomposition to study systems of two NHCs on an Au₅₈ cluster: two flat NHCs and a mixed configuration of one flat and one vertical NHC. The structure of these configurations can be found in the SI. A benefit of the sRBM decomposition is that it can distinguish the respective contributions of different molecules on the metal nanocluster. In Fig. 3, we show the comparison of the sRBM decomposition for the two different configurations of the two NHCs on the Au₅₈ cluster.

Fig. 3 sRBM decomposition for two NHCs on an Au₅₈ surface. (a) shows the decomposition for two NHCs in the flat configuration (F1 and F2), while (b) shows the decomposition for two NHCs in the flat and vertical configurations (F1 and V2), respectively.

Beginning with the two NHCs lying flat, we see that the p^total spectrum in Fig. 3(a) is similar to the p^total spectrum for the flat configuration in Fig. 2(a). As expected, we see that the intensities are increased by a factor of two. The sRBM decomposition still shows that the majority of the signal comes from the charge flow in the cluster and between the molecules and the cluster. Interestingly, the decompositions illustrate that the two flat NHCs, here referred to as F1 and F2, do not contribute equally to the signal. This is particularly noticeable for the peaks above 1400 cm⁻¹, where it can be seen that p^mol_F1 has the largest contributions. This is because F1 has a larger tilt angle relative to the surface than F2 does, which aligns the localized charge flow in F1 more closely with the local field, thus giving p^mol_F1 a stronger signal. It is also worth mentioning, that while p^cluster for one flat NHC compared to two NHCs is relatively similar, there is are noticeable differences for the peaks above 1400 cm⁻¹, which show the largest contributions from p^inter_F1. In contrast, the peaks around 600 cm⁻¹ are dominated by p^inter_F2. Overall, the sRBM shows that the two flat molecules behave fairly similarly due to the minimal interactions between the two molecules, although their individual contributions do depend sensitively on even small differences in tilt angles relative to the surface.

The sRBM decomposition for the mixed NHC configuration is shown in Fig. 3(b). In the following, we will refer to the flat configuration as F1 and the vertical configuration as V2. From the figure, we see that it is similar to the p^total spectrum for the vertical configuration in Fig. 2(b). This is expected as the surface selection rules dictate that the V2 configuration should have the largest Raman scattering and is further confirmed by the sRBM decomposition. The largest contributions come from p^mol, which is dominated by the contributions from the V2 configuration. The exception is the peak around 740 cm⁻¹, which is dominated by contributions from F1. This mode was also found to be strong in the single flat NHC on the surface, as well as for the configuration with two flat NHCs. Therefore, this mode serves as a good indicator of the flat configuration and is also seen in the experimental SERS spectrum for NHCs on gold surfaces.³³ Surprisingly, we see that p^inter is also dominated by contributions from V2. For the isolated flat configuration (or the two flat NHCs) the largest contribution is from p^inter. Therefore, the presence of the vertical configuration seems to reduce the interfacial charge flow for the flat configuration.

The sRBM decomposition of the NHCs on the metal surface illustrates how charge flow across the metal–molecule interface reflects the surface-selection rules. For molecules oriented perpendicular to the local field, the molecular vibrations are less effective in driving the charge flow, resulting in low enhancement. The larger enhancement for molecules oriented parallel with the local field occurs because the molecular vibrations effectively drive a coherent charge flow in the molecule and the surface. Finally, the sRBM decomposition shows that the interfacial charge-flow as reflected in p^inter is strongly influenced by interactions among molecules on the surface.

3.3 TERS images of benzene

To further demonstrate the usefulness of sRBM for understanding CT effects, we examine the modeling of single-molecule TERS images. TERS imaging relies on the high spatial resolution of the scanning tip combined with the strong local field enhancement to visualize individual vibrational modes of isolated molecules on a metal surface.^74–78 Since the current understanding of spatial resolution of the TERS images is mainly based on the confinement of the local field, most existing simulations focus on describing the EM only.^{79–81,81–84} However, several quantum mechanical models have been proposed to incorporate enhancement mechanisms beyond the EM, although often with simplified models for the tip due to the high computational cost.^85–89 Here, we will use the simplified damped-response method to conduct full quantum mechanical simulations of the TERS images such that both the EM and the CM are treated on an equal footing.

We will use a model system consisting of a Ag₂₀ tip interacting with a free-standing benzene molecule, in part due to its high symmetry and planar configuration, but also because it is a typical molecule used in the modeling of TERS images.^79,84,89 The sRBM decomposition of the TERS images for the 834 cm⁻¹ mode of benzene is shown in Fig. 4. Additional TERS images and their corresponding sRBM decomposition for the modes at 662 cm⁻¹ and 981 cm⁻¹ are shown in the SI. The incident frequency used in calculating the TERS images was chosen to be on resonance with the plasmon-like resonance of the tip at around 3.45 eV. Utilizing the symmetry of the vibrational modes, only half of the molecular plane ({(x,y)|0 ≤ x ≤ 4.8 Å,−4.8 ≤ y ≤ 4.8 Å}) was scanned. In this setup, a scanning height of 2.0 Å was used to be consistent with previous work.⁷⁹ A schematic of the tip–molecule structure is shown in Fig. 4(a).

Fig. 4 sRBM decomposition of the TERS images of the 834 cm⁻¹ mode of free-standing benzene. (a) A Ag₂₀ tip benzene model system. (b) TERS images generated from p^mol, decomposed into contributions from (c) the tip, (d) the molecule, and (e) the inter-fragment charge flow. The absolute values of sRBM decomposition are normalized with a scaling factor shown in the top right.

The sRBM decomposition of the TERS images is shown in Fig. 4(b)–(d). The TERS image based on p^total(Fig. 4(b)) is characterized by a dipolar-like pattern of the hotspot corresponding to the strong motion of the two H-atom vibrations that are out-of-plane but in opposite directions. From the decomposition we see that the large contributions to the signal are from p^tip, shown in Fig. 4(c), which is consistent with the EM mechanism due to the excitation of the plasmon-like resonance in the Ag₂₀ tip. Since the 834 cm⁻¹ mode itself is Raman inactive, the overall contribution from the molecule shown in Fig. 4(d) is very weak. The p^inter contribution from charge flow between the tip and the molecule, shown in Fig. 4(e), is also small but stronger than the contribution from the molecule. This is similar to what we found for the flat NHCs discussed above. Interestingly, the contributions from both the molecule and inter-fragment CT are destructive with respect to the charge flow in the tip. We also found this to be the case for the mode at 981 cm⁻¹ shown in the SI. The mode at 981 cm⁻¹ is also Raman inactive. However, for the Raman active mode at 662 cm⁻¹ we find constructive contributions from the tip, molecule and the inter-fragment charge flow leading to a high intensity. The modes at 834 cm⁻¹ and 981 cm⁻¹ are known to be enhanced by field-gradient effects.^79,90 The small contributions from the interfacial CT revealed by the sRBM decomposition support the idea that these modes likely obtain their enhancement from field-gradients effects. This is further supported by the observation that the hotspots in the TERS images are located away from the vibration atoms, whereas the contributions from p^mol and p^inter are localized closer to the molecule. Since the tip contributions are mainly attributed to the EM, while the molecule–tip contributions correspond to the CHEM, this difference indicates that the effects of CHEM are more localized than EM.

In the above simulations, the substrate was ignored since benzene only weakly interacts with the surface. However, recent TERS experiments have demonstrated a quenching effect of the TERS signal for molecules directly adsorbed on a metal surface.⁹¹ First-principles simulations have also identified an active role of the substrate in TERS images, although only for molecules that bond strongly to the underlying substrate.⁸⁹ Finally, tunneling between the tip and the substrate is also expected to be important for the determining the Raman scattering. Therefore, it is important to understand the role of the substrate in TERS imaging.

To do this we performed the sRBM decomposition for the Ag₂₀ tip interacting with benzene adsorbed on a small cluster and the results are shown in Fig. 5. To make the calculations feasible, a simple two-layer substrate consisting of 40 Ag atoms was used to support the scanning of the TERS image. See Fig. 5(a) for a schematic of the model system used. Utilizing the 2-fold symmetry, the overall TERS images were generated by mirroring the scanned results obtained from half the scanning area. The sRMB decomposition of the TERS images is shown in Fig. 5(b)–(h).

Fig. 5 sRBM decomposition of the TERS images of the 834 cm⁻¹ mode of benzene on a Ag substrate. (a) A Ag₂₀ tip benzene on a Ag₄₀ substrate model system. (b) TERS images generated from p^mol, decomposed into contributions from (c) the tip, (d) the molecule, and (e) the substrate. Inter-fragment Raman bond contributions of (f) molecule–tip, (g) molecule–substrate, and (h) tip–substrate. The absolute values of sRBM decomposition are normalized with a scaling factor shown in the top right.

As can be seen from Fig. 5(b), the TERS image based on p^total is different from that obtained from the free benzene molecule. Most noticeable, we see a distinct asymmetry between the two hotspots above and below the vibrating hydrogens. This can be attributed to the binding of the benzene molecule at a three-fold hollow site on coinage metals,⁹² which has also been reported in SERS studies of benzene on a silver substrate.⁹³ An illustration of the binding orientation can be seen in the TERS image in Fig. 5(b). This shows that the hydrogen atom near the strongest hotspot sits at a top-site, whereas the other hydrogen is closer to a three-fold hollow site. We also see from Fig. 5(b) that the intensity is significantly reduced compared to the TERS image of the free-standing benzene. This is consistent with the experimental observation of quenching of the TERS signal for molecules bound directly to a metal substrate.⁹¹ The absorption spectrum shown in the SI shows a significantly larger absorption peak for the tip–substrate configuration, likely due to the formation of the gap plasmon. Since the scanning height was chosen to be consistent with the free-molecule model, the distance between the tip and the substrate plane is close enough for quantum tunneling to occur, which is known to reduce the EM enhancement.^25–29 Therefore, the most likely explanation for the quenching is the chemical interactions between the molecule, substrate and tip.

To gain further insights into the TERS quenching and the asymmetric image, we will use the sRBM decomposition of the TERS image, which is shown in Fig. 5(c)–(e) for the tip, molecule, and substrate contributions, respectively. The inter-fragment contributions are shown in Fig. 5(f)–(h). We see from the TERS image based on p^tip that the asymmetry is less in the image based on p^tip but more pronounced in the images from the p^mol and the p^substrate contributions. Likewise, the inter-fragment contributions between the molecule and tip resemble that of the free-standing benzene, whereas the inter-fragment contributions between the molecule and substrate are asymmetric. The inter-fragment contributions between the tip and the substrate are also asymmetric. The contributions from p^mol are much stronger for benzene on the substrate than for free-standing benzene. This indicates that the interactions of benzene with the substrate break the symmetry of benzene and make the vibrational mode more Raman-allowed. Due to these combined results, the sRBM decomposition illustrates that the likely origin of the asymmetry in the TERS image is due to the CT interactions of the benzene with the substrate. Therefore, our results show that the TERS images are not only sensitive to the position of the tip but also reflect the electronic coupling to the substrate, even for molecules that are physisorbed.

From the decomposition, we see that the largest contributions still originate from p^tip; however, there are significantly larger contributions from both p^mol and p^substrate. We also see that inter-fragment contributions are much larger than for the free-standing benzene, indicating a more important role of interfacial CT. Interestingly, we find that there is also a significant contribution to the TERS images from the charge flow between the tip and the substate, indicating that quantum tunneling is important. The large quenching of the TERS image based on p^tip likely reflects the reduced EM field in the junctions due to the small separation between the tip and substate. This quenching is, to some degree, compensated for by the increase in the contributions from the charge flow between the tip, molecule, and substrate. We also note that most contributions are in-phase, especially around the strongest hotspot, which indicates that the molecular vibrations are more effective at driving the charge flow across the molecule–metal interface. However, it is not sufficient to overcome the large reduction in intensity from p^tip.

4 Conclusions

In this work, we demonstrated a combination of a simplified TDDFT method with a first-principles Raman bond model for the analysis and interpretation of SERS calculations. The simplified TDDFT method enables calculations of SERS spectra at a significantly reduced computational cost and has previously been shown to give calculated SERS spectra of similar quality to full TDDFT results. We denoted this new combination as the simplified Raman bond model (sRMB), which offers a computationally efficient approach for understanding SERS spectra when both the EM and the CM are important. Using the prototypical model system of pyridine interacting with an Ag₂₀ cluster, we first demonstrated that sRMB gives an interpretation of SERS spectra consistent with full TDDFT results. Specifically, the sRBM decomposition into contributions from the cluster, the molecule, and the inter-fragment charge-flow was shown to be in good agreement with the TDDFT results. The comparison demonstrated that sRBM provides reliable insights into the enhancement mechanisms and can be used to understand the importance of interfacial CT in SERS. The sRBM was then used to understand how interfacial CT influences the SERS spectra of NHCs with different orientations, as well as the modeling of single-molecule TERS images. For the NHCs, we showed that the surface-selection rules in SERS maps to the charge flow across the interface. Molecules that are aligned parallel with the local field leads to more coherent charge flow in the molecule and the surface, and thus a larger enhancement, whereas molecules that are oriented perpendicular cannot effectively drive this charge flow through their vibrational motion. Finally, the sRBM decomposition shows that the interfacial charge flow is strongly influenced by interactions among molecules on the surface. For the TERS images, the sRBM decomposition shows that the majority of the contributions come from the charge flow in the tip due to the EM mechanism. For benzene on a Ag substrate, we find that the TERS image is sensitive to the underlying lattice structure, which leads to an asymmetric image. Importantly, we see that the interaction with the substrate leads to a quenching of the TERS intensity, most likely due to quantum tunneling between the tip and the substrate. Overall, our results demonstrate that the sRBM offers a powerful approach for interpreting SERS spectra and providing new insights into interfacial CT.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data supporting the findings of this study are available within the article and its supplementary information (SI). The scripts to perform the Raman bond model analysis are available at https://github.com/jensengrouppsu/ramanbond. Additional raw data are available from the corresponding author upon reasonable request. Supplementary information: visual representations of the NHC systems, extended Raman bond decompositions, TERS images of Benzene, and sRBM of TERS images. See DOI: https://doi.org/10.1039/d5fd00117j.

Acknowledgements

The authors gratefully acknowledge financial support from National Science Foundation Grant CHE-2312222. Simulations in this work were in part performed on the Pennsylvania State University’s Institute for Computational and Data Sciences’ Roar supercomputer (https://icds.psu.edu/).

References

1. M. Moskovits, Rev. Mod. Phys., 1985, 57, 783–826.
2. A. Campion and P. Kambhampati, Chem. Soc. Rev., 1998, 27, 241–250.
3. K. Kneipp, H. Kneipp, I. Itzkan, R. R. Dasari and M. S. Feld, Chem. Rev., 1999, 99, 2957–2976.
4. S. Link and M. A. El-Sayed, Annu. Rev. Phys. Chem., 2003, 54, 331–366.
5. A. M. Michaels, M. Nirmal and L. E. Brus, J. Am. Chem. Soc., 1999, 121, 9932–9939.
6. H. Xu, E. J. Bjerneld, M. Käll and L. Börjesson, Phys. Rev. Lett., 1999, 83, 4357–4360.
7. J. A. Dieringer, R. B. Lettan, K. A. Scheidt and R. P. Van Duyne, J. Am. Chem. Soc., 2007, 129, 16249–16256.
8. K. Kneipp, Y. Wang, H. Kneipp, L. T. Perelman, I. Itzkan, R. R. Dasari and M. S. Feld, Phys. Rev. Lett., 1997, 78, 1667–1670.
9. S. M. Nie and S. R. Emory, Science, 1997, 275, 1102–1106.
10. A. Otto, Light Scattering in Solids IV: Electronics Scattering, Spin Effects, SERS, and Morphic Effects, Topics in Applied Physics, ed. M. Cardona and G. Güntherodt, Springer, Berlin/Heidelberg, 1984, vol. 54, pp. 289–418.
11. G. C. Schatz and R. P. Van Duyne, Handbook of Vibrational Spectroscopy, John Wiley and Sons, Ltd, New York, 2002, vol. 1, pp. 759–774.
12. L. Jensen, C. M. Aikens and G. C. Schatz, Chem. Soc. Rev., 2008, 37, 1061–1073.
13. S. M. Morton, D. W. Silverstein and L. Jensen, Chem. Rev., 2011, 111, 3962–3994.
14. D. V. Chulhai, Z. Hu, J. E. Moore, X. Chen and L. Jensen, Annu. Rev. Phys. Chem., 2016, 67, 541–564.
15. A. I. Pérez-Jiménez, D. Lyu, Z. Lu, G. Liu and B. Ren, Chem. Sci., 2020, 11, 4563–4577.
16. R. R. Jones, D. C. Hooper, L. Zhang, D. Wolverson and V. K. Valev, Nanoscale Res. Lett., 2019, 14, 231.
17. J. Langer, D. Jimenez de Aberasturi, J. Aizpurua, R. A. Alvarez-Puebla, B. Auguié, J. J. Baumberg, G. C. Bazan, S. E. J. Bell, A. Boisen, A. G. Brolo, J. Choo, D. Cialla-May, V. Deckert, L. Fabris, K. Faulds, F. J. García de Abajo, R. Goodacre, D. Graham, A. J. Haes, C. L. Haynes, C. Huck, T. Itoh, M. Käll, J. Kneipp, N. A. Kotov, H. Kuang, E. C. Le Ru, H. K. Lee, J.-F. Li, X. Y. Ling, S. A. Maier, T. Mayerhöfer, M. Moskovits, K. Murakoshi, J.-M. Nam, S. Nie, Y. Ozaki, I. Pastoriza-Santos, J. Perez-Juste, J. Popp, A. Pucci, S. Reich, B. Ren, G. C. Schatz, T. Shegai, S. Schlücker, L.-L. Tay, K. G. Thomas, Z.-Q. Tian, R. P. Van Duyne, T. Vo-Dinh, Y. Wang, K. A. Willets, C. Xu, H. Xu, Y. Xu, Y. S. Yamamoto, B. Zhao and L. M. Liz-Marzán, ACS Nano, 2020, 14, 28–117.
18. S.-Y. Ding, E.-M. You, Z.-Q. Tian and M. Moskovits, Chem. Soc. Rev., 2017, 46, 4042–4076.
19. E. C. Le Ru and P. G. Etchegoin, Principles of Surface-Enhanced Raman Spectroscopy, Elsevier, 2009.
20. E. Hao and G. C. Schatz, J. Chem. Phys., 2004, 120, 357–366.
21. W.-H. Yang, G. C. Schatz and R. P. Van Duyne, J. Chem. Phys., 1995, 103, 869–875.
22. B. T. Draine and P. J. Flatau, J. Opt. Soc. Am. A, 1994, 11, 1491–1499.
23. R. X. Bian, R. C. Dunn, X. S. Xie and P. T. Leung, Phys. Rev. Lett., 1995, 75, 4772–4775.
24. K. L. Kelly, E. Coronado, L. L. Zhao and G. C. Schatz, J. Phys. Chem. B, 2003, 107, 668–677.
25. J. A. Scholl, A. García-Etxarri, A. L. Koh and J. A. Dionne, Nano Lett., 2013, 13, 564–569.
26. K. J. Savage, M. M. Hawkeye, R. Esteban, A. G. Borisov, J. Aizpurua and J. J. Baumberg, Nature, 2012, 491, 574–577.
27. W. Zhu, R. Esteban, A. G. Borisov, J. J. Baumberg, P. Nordlander, H. J. Lezec, J. Aizpurua and K. B. Crozier, Nat. Commun., 2016, 7, 11495.
28. R. Esteban, A. G. Borisov, P. Nordlander and J. Aizpurua, Nat. Commun., 2012, 3, 825.
29. W. Zhu and K. B. Crozier, Nat. Commun., 2014, 5, 5228.
30. E. C. Le Ru, S. A. Meyer, C. Artur, P. G. Etchegoin, J. Grand, P. Lang and F. Maurel, Chem. Commun., 2011, 47, 3903–3905.
31. M. Moskovits and J. S. Suh, J. Phys. Chem., 1984, 88, 5526–5530.
32. M. Moskovits, J. Chem. Phys., 1982, 77, 4408–4416.
33. R. L. Thimes, A. V. B. Santos, R. Chen, G. Kaur, L. Jensen, D. M. Jenkins and J. P. Camden, J. Phys. Chem. Lett., 2023, 14, 4219–4224.
34. R. Li, I. Chaudhry, C. Tseng, S. Weng, Y. Wang, B. Zhao, I. Aravind, Z. Cai, J. Dawlaty, L. Jensen and S. B. Cronin, J. Phys. Chem. C, 2023, 127, 14263–14269.
35. J. R. Matthews, C. R. Shirazinejad, G. A. Isakson, S. M. E. Demers and J. H. Hafner, Nano Lett., 2017, 17, 2172–2177.
36. P. Z. El-Khoury and E. Aprà, J. Phys. Chem. C, 2020, 124, 17211–17217.
37. C. Van Dyck, B. Fu, R. P. Van Duyne, G. C. Schatz and M. A. Ratner, J. Phys. Chem. C, 2018, 122, 465–473.
38. S. M. Morton and L. Jensen, J. Am. Chem. Soc., 2009, 131, 4090–4098.
39. N. Valley, N. Greeneltch, R. P. Van Duyne and G. C. Schatz, J. Phys. Chem. Lett., 2013, 4, 2599–2604.
40. A. K. Kuhlman and A. T. Zayak, J. Phys. Chem. Lett., 2014, 5, 964–968.
41. A. T. Zayak, H. Choo, Y. S. Hu, D. J. Gargas, S. Cabrini, J. Bokor, P. J. Schuck and J. B. Neaton, J. Phys. Chem. Lett., 2012, 3, 1357–1362.
42. R. Chen and L. Jensen, J. Chem. Phys., 2020, 152, 024126.
43. R. Chen and L. Jensen, J. Chem. Phys., 2020, 153, 224704.
44. R. Chen and L. Jensen, Chem. Phys. Rev., 2023, 4, 021305.
45. F. L. Hirshfeld, Theor. Chim. Acta, 1977, 44, 129–138.
46. F. D. Proft, R. Vivas-Reyes, A. Peeters, C. V. Alsenoy and P. Geerlings, J. Comput. Chem., 2003, 24, 463–470.
47. A. Krishtal, P. Senet, M. Yang and C. Van Alsenoy, J. Chem. Phys., 2006, 125, 034312.
48. L. Gagliardi, R. Lindh and G. Karlström, J. Chem. Phys., 2004, 121, 4494–4500.
49. G. Hu and L. Jensen, J. Chem. Theory Comput., 2025, 21, 2546–2557.
50. E. J. Baerends, N. F. Aguirre, N. D. Austin, J. Autschbach, F. M. Bickelhaupt, R. Bulo, C. Cappelli, A. C. T. van Duin, F. Egidi, C. Fonseca Guerra, A. Förster, M. Franchini, T. P. M. Goumans, T. Heine, M. Hellström, C. R. Jacob, L. Jensen, M. Krykunov, E. van Lenthe, A. Michalak, M. M. Mitoraj, J. Neugebauer, V. P. Nicu, P. Philipsen, H. Ramanantoanina, R. Rüger, G. Schreckenbach, M. Stener, M. Swart, J. M. Thijssen, T. Trnka, L. Visscher, A. Yakovlev and S. van Gisbergen, J. Chem. Phys., 2025, 162, 162501.
51. A. D. Becke, Phys. Rev. A, 1988, 38, 3098–3100.
52. J. P. Perdew, Phys. Rev. B, 1986, 33, 8822–8824.
53. E. van Lenthe, E. J. Baerends and J. G. Snijders, J. Chem. Phys., 1993, 99, 4597–4610.
54. A. Bérces, R. M. Dickson, L. Fan, H. Jacobsen, D. Swerhone and T. Ziegler, Comput. Phys. Commun., 1997, 100, 247–262.
55. H. Jacobsen, A. Bérces, D. P. Swerhone and T. Ziegler, Comput. Phys. Commun., 1997, 100, 263–276.
56. S. K. Wolff, Int. J. Quantum Chem., 2005, 104, 645–659.
57. A. Ghysels, D. Van Neck, V. Van Speybroeck, T. Verstraelen and M. Waroquier, J. Chem. Phys., 2007, 126, 224102.
58. A. Ghysels, D. Van Neck and M. Waroquier, J. Chem. Phys., 2007, 127, 164108.
59. R. Rüger, A. Yakovlev, P. Philipsen, S. Borini, P. Melix, A. Oliveira, M. Franchini, T. van Vuren, T. Soini, M. de Reus, M. G. Asl, T. Teodoro, D. McCormack, S. Patchkovskii and T. Heine, AMS DFTB 2024.1, SCM, Theoretical Chemistry, Vrije Universiteit, Amsterdam, The Netherlands, 2024, https://www.scm.com.
60. S. Grimme, C. Bannwarth and P. Shushkov, J. Chem. Theory Comput., 2017, 13, 1989–2009.
61. L. Jensen, L. L. Zhao, J. Autschbach and G. C. Schatz, J. Chem. Phys., 2005, 123, 174110.
62. J. Neugebauer, M. Reiher, C. Kind and B. A. Hess, J. Comput. Chem., 2002, 23, 895–910.
63. L. Jensen, J. Autschbach and G. C. Schatz, J. Chem. Phys., 2005, 122, 224115.
64. R. Rüger, E. Van Lenthe, T. Heine and L. Visscher, J. Chem. Phys., 2016, 144, 184103.
65. N. Asadi-Aghbolaghi, R. Rüger, Z. Jamshidi and L. Visscher, J. Phys. Chem. C, 2020, 124, 7946–7955.
66. L. Zhao, L. Jensen and G. C. Schatz, J. Am. Chem. Soc., 2006, 128, 2911–2919.
67. S. Chowdhury, G. Hu, I. M. Jensen, A. V. B. Santos, D. M. Jenkins, L. Jensen and J. P. Camden, J. Phys. Chem. C, 2024, 128, 13550–13557.
68. R. Chen and L. Jensen, J. Chem. Phys., 2022, 157, 184705.
69. G. Kaur, R. L. Thimes, J. P. Camden and D. M. Jenkins, Chem. Commun., 2022, 58, 13188–13197.
70. T. Weidner, J. E. Baio, A. Mundstock, C. Große, S. Karthäuser, C. Bruhn and U. Siemeling, Aust. J. Chem., 2011, 64, 1177–1179.
71. P. Bellotti, M. Koy, M. N. Hopkinson and F. Glorius, Nat. Rev. Chem., 2021, 5, 711–725.
72. C. R. Larrea, C. J. Baddeley, M. R. Narouz, N. J. Mosey, J. H. Horton and C. M. Crudden, ChemPhysChem, 2017, 18, 3536–3539.
73. A. Inayeh, R. R. K. Groome, I. Singh, A. J. Veinot, F. C. de Lima, R. H. Miwa, C. M. Crudden and A. B. McLean, Nat. Commun., 2021, 12, 4034.
74. R. Zhang, Y. Zhang, Z. C. Dong, S. Jiang, C. Zhang, L. G. Chen, L. Zhang, Y. Liao, J. Aizpurua, Y. Luo, J. L. Yang and J. G. Hou, Nature, 2013, 498, 82–86.
75. J. Stadler, T. Schmid and R. Zenobi, Nano Lett., 2010, 10, 4514–4520.
76. N. Jiang, N. Chiang, L. R. Madison, E. A. Pozzi, M. R. Wasielewski, T. Seideman, M. A. Ratner, M. C. Hersam, G. C. Schatz and R. P. Van Duyne, Nano Lett., 2016, 16, 3898–3904.
77. N. Kumar, W. Su, M. Veselý, B. M. Weckhuysen, A. J. Pollard and A. J. Wain, Nanoscale, 2018, 10, 1815–1824.
78. Z. He, Z. Han, M. Kizer, R. J. Linhardt, X. Wang, A. M. Sinyukov, J. Wang, V. Deckert, A. V. Sokolov, J. Hu and M. O. Scully, J. Am. Chem. Soc., 2019, 141, 753–757.
79. P. Liu, D. V. Chulhai and L. Jensen, ACS Nano, 2017, 11, 5094–5102.
80. P. Liu, X. Chen, H. Ye and L. Jensen, ACS Nano, 2019, 13, 9342–9351.
81. Y. Zhang, B. Yang, A. Ghafoor, Y. Zhang, Y.-F. Zhang, R.-P. Wang, J.-L. Yang, Y. Luo, Z.-C. Dong and J. G. Hou, Natl. Sci. Rev., 2019, 6, 1169–1175.
82. X. Chen, P. Liu, Z. Hu and L. Jensen, Nat. Commun., 2019, 10, 2567.
83. S. Duan, G. Tian, Y. Ji, J. Shao, Z. Dong and Y. Luo, J. Am. Chem. Soc., 2015, 137, 9515–9518.
84. M. Takenaka, T. Taketsugu and T. Iwasa, J. Chem. Phys., 2021, 154, 024104.
85. K. Fiederling, M. Abasifard, M. Richter, V. Deckert, S. Kupfer and S. Gräfe, ACS Nano, 2023, 17, 13137–13146.
86. F. Latorre, S. Kupfer, T. Bocklitz, D. Kinzel, S. Trautmann, S. Grafe and V. Deckert, Nanoscale, 2016, 8, 10229–10239.
87. K. Fiederling, M. Abasifard, M. Richter, V. Deckert, S. Gräfe and S. Kupfer, Nanoscale, 2020, 12, 6346–6359.
88. R. L. M. Gieseking, J. Phys. Chem. C, 2022, 126, 11690–11700.
89. Y. Litman, F. P. Bonafé, A. Akkoush, H. Appel and M. Rossi, J. Phys. Chem. Lett., 2023, 14, 6850–6859.
90. D. V. Chulhai and L. Jensen, J. Phys. Chem. C, 2013, 117, 19622–19631.
91. B. Yang, G. Chen, A. Ghafoor, Y.-F. Zhang, X.-B. Zhang, H. Li, X.-R. Dong, R.-P. Wang, Y. Zhang, Y. Zhang and Z.-C. Dong, Angew. Chem., Int. Ed., 2023, 62, e202218799.
92. W. Liu, V. G. Ruiz, G.-X. Zhang, B. Santra, X. Ren, M. Scheffler and A. Tkatchenko, New J. Phys., 2013, 15, 053046.
93. M. Moskovits and D. P. DiLella, J. Chem. Phys., 1980, 73, 6068–6075.

---

Footnote

† These authors contributed equally.

---