Vibrational Energy Dissipation in Noncontact Single-Molecule Junctions Governed by Local Geometry and Electronic Structure

TL;DR

This paper studies how vibrations in single-molecule junctions dissipate energy through interactions with the substrate, using theory and simulations to reveal how geometry and electronic structure influence this process.

Contribution

A novel approach combining DFT, machine learning, and nonadiabatic molecular dynamics to analyze vibrational energy dissipation in noncontact single-molecule junctions.

Findings

Vibrational mode specificity is governed by electron–phonon and phonon–phonon coupling.

Electron–phonon relaxation rates vary by two orders of magnitude depending on tip–substrate geometry.

Electron–phonon coupling enhances phonon–phonon coupling in a weak nonadditive effect.

Abstract

The vibrational dynamics of adsorbate molecules in single-molecule
junctions depend critically on the geometric structure and electronic
interactions between the molecule and the substrate. Vibrations, excited
mechanochemically or by external stimuli, dissipate energy into substrate
electrons and phonons. Energy dissipation leads to the broadening
of spectral lines, vibrational lifetimes, and the coupling between
molecular and substrate phonons. It affects molecular manipulation,
giving rise to nanoscale friction, and contributes to scanning probe
and surface spectroscopy signals. We present an approach to disentangle
adsorbate vibrational dynamics in noncontact junctions by employing
Density Functional Theory, machine learning, and nonadiabatic molecular
dynamics. Focusing on the CO-functionalized Cu surfaces representing
a single-molecule junction, a widely studied system in scanning probe
and energy dissipation experiments, we reveal strong vibrational mode
specificity governed by the interplay of electron–phonon and
phonon–phonon coupling. Electron–phonon relaxation rates
vary by 2 orders of magnitude between modes and sensitively depend
on the tip–substrate geometry. We find evidence of a weak nonadditive
effect between both energy dissipation channels, where electron–phonon
coupling enhances phonon–phonon coupling. Our predicted vibrational
lifetimes agree with infrared spectroscopy and helium scattering experiments.
Finally, we outline how our findings can inform and enhance spectroscopy
and scanning probe experiments.

Figures

Schematic noncontact single-molecule junction experiment where a CO molecule on a surface is probed with a CO-decorated tip. Van-der-Waals (VdW) interactions and Pauli repulsion govern the interaction between sample and probe. Covalent bonds are present between the CO molecule and the metallic component. Energy exchange occurs primarily between chemically bonded components, in the form of EPC and PPC. a) Molecule is attached to the tip. b) Molecule is attached to the surface. c) Different adsorption environments considered in this work. All systems are calculated with periodic boundary conditions, with the depicted tip geometries extended by a further four layers of metal slab (not shown here).

Four eigenmodes of CO on Cu surface: frustrated rotation (FR); frustrated translation (FT); surface adsorbate stretch (SA); internal stretch of CO bond (IS).

VDOS of the four different adsorption environments. The total VDOS is shown in black. The contribution of the CO modes is shown in orange.

Coupling strength of the FR mode to other vibrational modes as a result of PPC.

DOS projected at the Cu atom to which the CO molecule bonds. Flatter surface geometry shows a more delocalized projected DOS.

a) Flowcharts showing the three different workflows used to calculate relaxation rates in this work. b) Relaxation rates of CO vibrational modes on different adsorption environments determined with the three methods shown in panel a. c) Contributions to relaxation rates by EPC and PPC when determined using the kinetic energy decay method: Light brown: electron–phonon relaxation rate extracted from MDEF; dark-brown: phonon–phonon relaxation rate determined with MD; slate blue: nonadditive difference of the sum of ECP and PPC with respect to total relaxation rate; black: total relaxation rate from MDEF.

Tables

1. Vibrational Energies of CO Molecule in Different Adsorption Geometries; Experimental Vibrational Energies of the CO Molecule on a Cu(111) Surface ,

	vibrational energies/meV
	FT	FR	SA	IS
sharp tip	2.74	31.38	49.06	271.48
blunt tip	2.89	32.28	50.35	270.86
slab adatom	3.17	31.74	49.13	272.26
slab top	5.29	35.92	47.97	264.51
slab top exp. ,	4.07	36.47	41.41	257.63

2. CO Adsorbate Vibrational Lifetimes due to PPC

	Lifetime/ps
	FT	FR	SA	IS
sharp tip	142.1	69.0	56.2	48.5
blunt tip	110.9	67.0	31.9	65.4
slab adatom	59.8	65.3	53.7	46.3
slab	19.3	24.6	26.8	53.7

3. Electron–Phonon-Induced Vibrational Lifetimes in the Harmonic Approximation, Including a Comparison to Theoretical Values Previously Reported in Literature

	Lifetime/ps
	FT	FR	SA	IS
sharp tip	185.9	21.2	80.0	2.8
blunt tip	99.1	6.0	39.0	3.3
slab adatom	110.8	10.8	52.9	4.9
slab (Cu(111))	47.9	3.6	12.5	1.9
CO on Cu(111)	84.4	2.1	7.8	2.2
CO on Cu(110)	76.9	2.7	10.8	2.2
CO on Cu(100)	125.0	2.6	20.0	3.2

4. Vibrational Lifetimes due to EPC and PPC Calculated by MDEF Equilibrium Correlation Analysis Compared to Experimental Lifetimes; Asterisks (*) Indicate that Experiments Correspond to Lower Bounds on the Vibrational Lifetime

	Lifetime/ps
	FT	FR	SA	IS
sharp tip	40.4	24.1	21.1	5.2
blunt tip	30.0	4.1	15.2	9.0
slab adatom	21.5	8.5	28.4	6.3
slab	19.8	3.6	20.7	3.6
CO on Cu(111)	-	2.5*	2.5*	2.5*
CO on Cu(100)	-	2.5*	1.9*	0.7*
CO on Cu(100)	8	-	-	-

Full text

Introduction

The vibrational dynamics of atoms and molecules adsorbed on a surface provide valuable insights into the interactions and energy exchange between adsorbate and substrate. These dynamics are governed by interactions between adsorbates and the surface and quantum-mechanical energy dissipation mechanisms. ?,? Energy dissipation mechanisms in noncontact single-molecule junctions significantly affect measured signals in atomic force microscopy (AFM), lateral force microscopy (LFM), ?,? scanning tunneling microscopy (STM),? inelastic electron tunneling spectroscopy (IETS), and infrared spectroscopy experiments. ?,? In such experiments, a key challenge is to disentangle geometric structure, dynamics, and energy dissipation mechanisms in the measured signal: For example, AFM-induced hopping of atoms and molecules leads to measurable energy dissipation, ?,?−? ? yet the underlying dissipation pathways remain unclear. Helium atom scattering (HAS) experiments of CO adsorbed on Cu surfaces have found a strong dependence of vibrational frequencies on both the surface facet and the presence of defects.?

It is well-established that the surface geometry plays a critical role in shaping molecule–surface interactions and energy transfer processes: Previous studies have shown the effect of defects on static interfacial charge transfer,? while adsorption site and facet affect the vibrational properties of adsorbates.? Establishing an understanding of how surface geometry alters different energy dissipation channels and affects the vibrational dynamics of molecules adsorbed on surfaces remains an important research area.?

Figure shows a schematic setup of an idealized single-molecule junction as it would be present in an AFM experiment. If vibrations of the molecule in the junction are excited, for example, through the motion of the tip, energy dissipation away from the molecule can occur by excitation of phonons or electron–hole-pair (EHP) excitations in the metal tip or substrate. The former is due to anharmonicity that induces phonon–phonon coupling (PPC) between adsorbate and substrate.? The latter is due to coupling between electronic and nuclear motions, commonly referred to as electron–phonon coupling (EPC).? Intramolecular vibrational relaxation and radiative processes, such as infrared emission, also contribute to energy dissipation, though their role is typically minor for small adsorbates. ?,? However, for larger adsorbates, intramolecular vibrational relaxation can become a significant dissipation channel.? Consequently, it is of great interest to explore how the geometric structure of a surface or junction affects the mechanical and electronic interactions with the adsorbate, which govern the relative contributions of EPC and PPC to the vibrational dynamics. Gaining control over these energy dissipation channels requires the ability to identify their distinct signatures and disentangle their interplayan essential step toward rational design of surface functionality at the atomic scale.

Schematic noncontact single-molecule junction experiment where a CO molecule on a surface is probed with a CO-decorated tip. Van-der-Waals (VdW) interactions and Pauli repulsion govern the interaction between sample and probe. Covalent bonds are present between the CO molecule and the metallic component. Energy exchange occurs primarily between chemically bonded components, in the form of EPC and PPC. a) Molecule is attached to the tip. b) Molecule is attached to the surface. c) Different adsorption environments considered in this work. All systems are calculated with periodic boundary conditions, with the depicted tip geometries extended by a further four layers of metal slab (not shown here).

In this work, we study competing energy dissipation mechanisms and their effect on the vibrational dynamics of adsorbate molecules for the exemplary case of CO-decorated tip–surface copper junctions. CO on Cu is a prototypical surface–adsorbate system. It serves as a benchmark in surface science, enabling the systematic investigation of fundamental processes such as adsorption, diffusion, and vibrational energy dissipation. ?,?,?−? ? Moreover, it is widely used in atomic-resolution scanning probe measurements.? Owing to the similar electronic structure of coinage metals, insights from CO on Cu are transferable to Ag and Au surfaces, making it relevant for experiments and applications in catalysis and energy conversion, where electron–phonon and phonon–phonon coupling govern surface reactivity.? Various previous theoretical studies of energy dissipation for CO on different, idealized Cu surfaces ?−? ? have been reported, but no clear analysis of the role of the local adsorption geometry exists to date. A specific challenge in this system for electronic structure theory is the correct description of the potential energy surface of CO on Cu(111). ?−? ? We address this by using the state-of-the-art screened hybrid functional HSE06? in combination with the MBD-NL van der Waals (vdW) correction.? From a first-principles perspective, we describe the details and main energy dissipation mechanisms, namely EPC and PPC, which govern the dynamics at the atomic scale.? We find that the surface geometry plays a significant role in these mechanisms, with flatter surfaces exhibiting stronger EPC and PPC. Additionally, we observe a strong mode dependence of EPC, with relaxation rates ranging over 2 orders of magnitude. The molecular dynamics with electronic friction (MDEF) approach ?,? allows us to study the interplay of the two dissipation channels and their concerted effects on the relaxation dynamics. We find signs of weak nonadditivity between the two effects, suggesting that subtle dynamical steering may arise from one or both mechanisms. Our results are in good agreement with existing experiments and support the interpretation of single-molecule dynamics in noncontact junctions.

Results

Models of the Tip Geometry in Noncontact Junctions

As shown in Figure, we can consider two interaction regimes for a noncontact CO tip-molecule–surface junction: In one case, the molecule is attached to the surface, in the other case, it is attached to the tip. The molecule can be excited by a variety of processes: The motion of the tip will displace the molecule with respect to the corrugated potential energy landscape of the underlying substrate.? The mechanochemical forces the molecule experiences in the process will excite its vibrations. The molecular vibrations can absorb light,? supported by field enhancement in the junction.? A potential bias and electron tunneling can lead to current-induced forces.? In all cases, molecular vibrations will be excited, leading to energy dissipation into EHPs and phonons in the metal. The case where the molecule is attached to the surface is not too different from the scenario in the absence of the tip. The tip and molecule are separated by a distance of multiple Angstroms, so mechanochemical forces on the molecule will be weak. Other perturbation scenarios have previously been studied. ?,? For simplicity, we consider only surfaces in both the slab and slab-adatom systems that correspond to the Cu(111) facet.

We turn our attention to the second case, where the molecule is attached to the tip. Regardless of the origin of the vibrational excitation, energy dissipation into the underlying substrate in this scenario will be almost negligible, as evidenced by calculations for CO on Cu in Section S1. Therefore, molecular vibrations will dominantly dissipate energy into EHPs and phonons in the tip. This motivates our model design: We consider four different possible adsorption environments: CO on a sharp tip (sharp tip), CO on a blunt tip (blunt tip), CO on an adatom on a surface slab (slab adatom), and CO on a slab (slab), which are depicted in Figuresc and S3. The underlying surface, in the case of lateral motion of the tip, only serves as a corrugated energy landscape that triggers molecular vibrations? and therefore does not have to be explicitly considered in the atomistic model. The “tip” models are therefore represented by two-dimensional periodic surface slab geometries. Further details on the models and the computational methodology are provided in the Methods and Section S1.

Vibrational Mode Analysis

CO adsorbed on a flat metal surface has four distinct types of vibrational modes: frustrated translation (FT), frustrated rotation (FR), surface–adsorbate stretch (SA), and the internal stretch of the CO bond (IS), as illustrated in Figure. The FT and FR modes consist of two degenerate modes each.

Four eigenmodes of CO on Cu surface: frustrated rotation (FR); frustrated translation (FT); surface adsorbate stretch (SA); internal stretch of CO bond (IS).

We determine vibrational energies from second-order force constants (i.e., the Hessian) for all atoms in the molecule and the top layers of the four surface models. Simulations include ten layers of freely movable Cu atoms, spanning hundreds of degrees of freedom. Energies, forces, and second-order force constants are evaluated with a MACE graph neural network model? trained on first-principles data (see Methods). The MACE potential accurately reproduces the underlying HSE06 data, correctly predicting the adsorption sites: the top site for CO on a (111) surface slab of copper and the hollow site for the adatom in the CO-adatom arrangement (see Section S2). The calculated vibrational frequencies for CO on the Cu (111) slab show good agreement with experimental data, as presented in Table. Interestingly, the stiffness of low-energy modes increases as the tip becomes blunter, as indicated in Table. Moreover, there is a jump in stiffness when transitioning from CO adsorbed on an adatom to CO adsorbed directly on the surface. HAS experiments estimated the energy of the FT modes on a flat Cu(111) surface to be 4.01 meV, while the energy dropped to 3.00 meV at step edges of the Cu(211) and Cu(511) surfaces.? This trend is in good agreement with our finding that more tapered (or corrugated) adsorption environments exhibit lower vibrational energies for the FT mode.

1: Vibrational Energies of CO Molecule in Different Adsorption Geometries; Experimental Vibrational Energies of the CO Molecule on a Cu(111) Surface ,

The vibrational density of states (VDOS) (see Figure) provides an overview of all vibrational modes. The VDOS is determined from the same molecular dynamics (MD) simulations that yield the results presented in Phonon–. We differentiate between modes localized dominantly at the molecule (orange in Figure), where the majority of the kinetic energy is contained in the CO and substrate modes (black in Figure), where most of the energy is in the motion of the Cu atoms. The FR, SA, and IS vibrational modes of CO exhibit higher vibrational energy compared to most copper substrate modes, which typically range between approximately 5 and 30 meV. We note that the VDOS shown here is normalized by the number of atoms. The substrate vibrational states are more localized in energy (or frequency) for slab and slab-adatom surfaces and more delocalized for the sharp and blunt tips. For sharp and blunt tips, a significant portion of the surrounding space is vacuum, whereas this space is occupied by atoms in the slab and slab-adatom surfaces. The presence of more atoms in these flatter adsorption environments results in a greater number of substrate vibrational states per volume. This facilitates stronger overall coupling between the CO and substrate vibrational modes, provided that the CO and Cu degrees of freedom are sufficiently coupled. As a rule of thumb, if the off-diagonal elements of the Hessian (for CO-Cu coupling) are of a similar magnitude, we can expect strong coupling. As shown in Figure S4, this is indeed the case. In the case of the sharp and blunt tip, we find a strong coupling for the CO to a small number of Cu degrees of freedom. The slab adatom and slab geometries exhibit weaker coupling but a larger number of states. When discussing PPC, we will return to this observation.

VDOS of the four different adsorption environments. The total VDOS is shown in black. The contribution of the CO modes is shown in orange.

Phonon–Phonon Coupling

The atomic motion of the CO molecule is coupled to the atomic motion of the Cu substrate through phonon–phonon interactions and to EHP through EPC.? Even at low temperatures, PPC plays a significant role in the dynamics of molecular vibrations. In our study, we use a temperature of 5 K, which is commonly employed in experiments such as AFM measurements with CO-decorated tips. PPC is determined by analyzing the dynamics of a CO molecule on adsorption environments where approximately 500 Cu atoms are free to move per unit cell, allowing energy transfer between CO modes and substrate modes. To determine coupling strengths and vibrational lifetimes, we use vibrational analysis of MD trajectories in thermodynamic equilibrium.? We refer to this approach as equilibrium correlation analysis. A conceptually related approach, which combines vibrational projection, statistical analysis, and local quantum vibrational embedding, was recently introduced to investigate ultrafast energy transfer pathways from lattice phonons to intramolecular modes in molecular crystals.? We first thermalize each system at a temperature of 5 K and then allow it to evolve in thermodynamic equilibrium over time. Over time, the velocity time series of the atoms becomes decorrelated, as the atoms transfer mechanical energy between each other. The rate of this decorrelation can be found in the cross-correlation function. For this analysis, we employ normal mode decomposition to project the velocity time series onto the basis defined by the vibrational eigenvectors. For the velocity time series of two modes v i(t) and v j(t) the velocity cross-correlation function C ij(τ) is as follows:

$$[eqn]$$

The Fourier transform of the velocity cross-correlation function yields the cross-spectrum:

$$[eqn]$$

The cross-spectrum provides valuable insights into vibrational lifetimes and relative coupling strengths. Let us first discuss the coupling strength, which can be determined by analyzing the off-diagonal elements C̃ ij(ω). These exhibit two Lorentzian-shaped peaks at the characteristic frequencies of modes i and j, with their heights reflecting the coupling strength between these modes (see Section S6.2). To further understand this interaction, we examine the coupling density ρ_j_(ω), defined as the sum of the cross spectra C̃ ij(ω) for all modes j correlated with the mode of interest i:

$$[eqn]$$

A broad peak in the coupling density at a given frequency means that modes at this frequency are correlated with the mode of interest, ω_i_, indicating strong vibrational coupling. Our analysis reveals that the coupling is strongest with adjacent substrate vibrational modes, as shown in Figures and S24–S27. This effect is less pronounced in case of the energetically decoupled SA and IS modes. The coupling between different CO modes is small in comparison, indicating that intramolecular vibrational redistribution plays only a minor role. Notably, the FR and IS modes exhibit the strongest coupling.

Coupling strength of the FR mode to other vibrational modes as a result of PPC.

The geometry of the surface has a substantial impact on the coupling density. Among the different geometries, the sharp tip exhibits the weakest coupling, while the slab geometry has the strongest coupling density. This observation can be attributed to the VDOS discussed earlier. As surface geometry becomes more atomically flat, the VDOS increases, offering more vibrational states to accommodate the energy of CO modes and thereby enhancing PPC.

The diagonal terms C̃ ii(ω) of the cross-correlation function contain the power spectrum. The power spectrum exhibits a single Lorentzian-shaped peak at the characteristic frequency of the mode i, where the full width at half-maximum Δf i is inversely proportional to the lifetime of mode τ_i_.

$$[eqn]$$

As shown in Table, the overall lifetimes tend to decrease as the surface becomes more atomically flat. Interestingly, this trend does not apply uniformly to all modes. The FT mode shows the strongest dependence on surface geometry, while this dependence is less pronounced for the FR and SA modes. The lifetime of the IS mode remains largely unaffected by changes in surface shape. Overall, the influence of the surface geometry is strongest in low-energy modes and becomes weaker as the energy of the mode increases. The largest phononic lifetime (FT on sharp tip model) is one order of magnitude larger than the smallest lifetime (FT on slab model).

2: CO Adsorbate Vibrational Lifetimes due to PPC

Electron–Phonon Coupling

EPC is a significant factor in vibrational energy dissipation on metal surfaces. ?,? Due to the absence of a band gap, the Born–Oppenheimer approximation does not strictly hold, and the motion of adsorbate atoms can induce EHP excitations that lead to nonadiabatic effects, sometimes also called electronic friction effects. In the weak perturbation regime, e.g., where adsorbate motion is triggered mechanochemically or through vibrational excitation, this effect leads to nonadiabatic energy dissipation. To describe the impact of these excitations on dynamics, different simulation methods have been developed. Time-dependent perturbation theory has emerged as a feasible approach for extended systems. This method enables the calculation of the electron–phonon-induced vibrational lifetimes (see Description of Electron–Phonon Coupling in the Methods),? which we discuss in the following.

We evaluate EPC at a temperature of 5 K, typical for noncontact single-molecule junction experiments. At such low temperatures, EPC dominates over PPC in metals ?,? and first-order EPC is the dominant contribution to vibrational lifetimes, allowing higher-order EPC effects to be neglected.? Its strength is sensitive to the electronic states that interact with nuclear motion, which we characterize by the electronic density of states (DOS) projected onto the copper atom bonded to CO (PDOS), which is obtained from Density Functional Theory (DFT). The PDOS is highly sensitive to the geometry of the system. Since the energy to excite electrons is supplied by vibrational states with frequency ω, only a small portion of the DOS near the Fermi level (E F ± ℏω)relative to the energy scale of electronic statesis accessible. Therefore, analyzing the DOS at the Fermi level is instructive for understanding the strength of EPC. As shown in Figure, the PDOS of the sharp and blunt tips exhibits states localized in energy. The magnitude of the PDOS close to the Fermi level is small for the sharp tip, while the blunt tip exhibits energetically localized spikes. In contrast, the slab and slab-adatom systems have a PDOS that is delocalized in energy. The slab PDOS exhibits a larger magnitude compared to that of the adatom on the slab. Since a larger magnitude of the PDOS at and close to the Fermi level allows for stronger EPC, we expect that EPC will be stronger at flatter surfaces.

DOS projected at the Cu atom to which the CO molecule bonds. Flatter surface geometry shows a more delocalized projected DOS.

The strength of EPC is reflected in the vibrational lifetimes, where shorter lifetimes indicate stronger coupling (see Table). The FT and SA modes experience relatively low EPC (resulting in longer lifetimes), whereas FR and IS modes exhibit EPC approximately an order of magnitude higher (with lifetimes an order of magnitude shorter) compared to FT and SA. Previous studies of CO on different Cu facets ?−? ? have reported similar trends of mode dependence, as we show in Table. Our results for the slab are in good agreement with literature values. In the case of the slab, there is an ongoing debate about whether the FR or the IS decays faster. Indeed, the lifetimes reported here and in other literature are very similar for these modes. However, as the tip becomes more corrugated, the IS mode clearly becomes the fastest dissipation mode, while the lifetime of the FR mode increases. As shown in Table, vibrational lifetimes are significantly longer on sharper tips, with CO on a flat surface exhibiting the shortest lifetimes. This trend aligns with the reduced magnitude of the PDOS close to the Fermi level on more corrugated surfaces.

3: Electron–Phonon-Induced Vibrational Lifetimes in the Harmonic Approximation, Including a Comparison to Theoretical Values Previously Reported in Literature

Electron–Phonon and Phonon–Phonon Coupling

The experimentally measured lifetime is governed by both EPC and PPC. To capture this combined effect, we employ Langevin dynamics, specifically MDEF.? We use a classical approximation, neglecting the zero-point energy of the phonons. This approach has been successfully applied to the dynamics of atoms ?,? and molecules ?,? at metal surfaces. In MDEF, EPC is treated as a friction force represented by the electronic friction tensor (see eq). The electronic friction tensor is calculated in the Markovian quasi-static limit (see Electron–Phonon Coupling). Here, the frictional force experienced by the nuclei depends only on their instantaneous velocities, and memory effects are ignored. This assumes that the electronic structure adjusts instantaneously to the nuclear configuration. In the quasi-static limit, the electronic friction tensor is computed at zero frequency in the harmonic approximation for a given nuclear geometry without explicitly propagating electronic dynamics. This approach is commonly employed for MDEF on metals. ?,? We determine vibrational relaxation rates and lifetimes using three different methods (see Figurea):

a) Flowcharts showing the three different workflows used to calculate relaxation rates in this work. b) Relaxation rates of CO vibrational modes on different adsorption environments determined with the three methods shown in panel a. c) Contributions to relaxation rates by EPC and PPC when determined using the kinetic energy decay method: Light brown: electron–phonon relaxation rate extracted from MDEF; dark-brown: phonon–phonon relaxation rate determined with MD; slate blue: nonadditive difference of the sum of ECP and PPC with respect to total relaxation rate; black: total relaxation rate from MDEF.

Equilibrium Correlation Analysis

Lifetimes are determined by analyzing the velocity cross-correlation function from an MD (without EPC) or MDEF (with EPC) trajectory in thermodynamic equilibrium. This method is introduced and applied in Phonon–Phonon Coupling to determine PPC relaxation rates and lifetimes from MD trajectories. When EPC is included via MDEF, the method yields total relaxation rates that capture both PPC and EPC.

Table shows the vibrational lifetimes from equilibrium correlation analysis of the MDEF simulation results. While the IS and FR lifetimes are dominated by EPC, the FT mode is dominated by PPC. The two effects are somewhat equal in contribution in the case of the SA mode lifetime. Lifetimes for the FR and IS modes are in good agreement with infrared spectroscopy experiment. However, the experimental lifetime of the SA mode is 2.5 ps, while the calculated lifetime is 20.7 ps for the Cu(111) surface. A similar discrepancy has been reported in literature and was attributed to inhomogeneous broadening:? Experiments conducted on rougher surfaces show a greater sensitivity of the mode amplitude and line width of the SA mode (and the amplitude of the FR mode) than for the IS mode. The literature lifetimes from infrared spectroscopy experiments should be considered as lower bounds (marked by asterisks (*) in Table). This is because the observed line widths are additionally broadened due to limitations of the measurement setup. Since lifetime is inversely related to line width, this broadening leads to a systematic underestimation of the lifetimes. Additionally, helium scattering experiments have measured the lifetime of the FT mode of CO on Cu(100) to 8 ps,? which is in the same order of magnitude as the lifetime for the FT mode of CO on a Cu(111) slab reported here.

4: Vibrational Lifetimes due to EPC and PPC Calculated by MDEF Equilibrium Correlation Analysis Compared to Experimental Lifetimes; Asterisks (*) Indicate that Experiments Correspond to Lower Bounds on the Vibrational Lifetime

Kinetic Energy Decay

We can compare lifetimes determined through equilibrium correlation analysis with those explicitly calculated from nonequilibrium kinetic energy decay simulations (see Methods). Hereby, lifetimes are predicted by monitoring how the kinetic energy of an initially excited vibration mode decreases over time during a nonequilibrium MDEF simulation. We excite vibrations with a specific amount of energy and monitor the decay of the kinetic energy (Figure S18). To determine the lifetimes, we assume that the vibrational modes behave like damped harmonic oscillators, where the envelope of the kinetic energy is assumed to exhibit an exponential decay. By including EPC using MDEF, total relaxation rates can be obtained. Overall, we observe that the kinetic energies of the CO vibrational modes indeed show an exponential decay (see Figure S18).

Additive Relaxation Rates

We sum two independent contributions to the vibrational relaxation rate: the PPC relaxation rates obtained from MD equilibrium correlation analysis, and the EPC relaxation rates calculated in the quasi-static limit using first-order time-dependent perturbation theory. The EPC contribution calculated in this way is described in detail in Electron–Phonon Coupling. The total relaxation rate γ_tot_ is the inverse of the lifetime γ_ x _ = 1/τ_ x _ and can be approximated by summing the contributions from PPC γ_ph‑ph_ (Table) and EPC γ_e‑ph_ (Table):

$$[eqn]$$

Figureb compares the three methods of determining total lifetimes (and relaxation rates). Additional technical details on the three methods are provided in Section S3. Within the estimated uncertainty of our simulations, the results are in good agreement. All three methods yield lifetimes that exhibit similar trends across the different systems and modes. CO adsorption on a terrace leads to overall larger relaxation rates and smaller lifetimes compared to adsorption on adatoms or single atom tips. Relaxation rates exhibit a strong mode dependence, with the largest rate being an order of magnitude greater than the smallest. The numerical values in Figureb can be found in Section S3.

Contributions to the Relaxation Rates

We decompose the total vibrational relaxation rate into EPC and EPP contributions to identify the dominant dissipation channel for each vibrational mode in each adsorption environment. This is done using the kinetic energy decay method. The phonon–phonon relaxation rates are obtained by monitoring the decay in kinetic energy during MD simulations (which do not include EPC). The electron–phonon relaxation rate is extracted from the MDEF trajectories. This is accomplished by integrating the electronic friction and random force components in eq along the trajectory.

Results for all adsorption environments are presented in Figure S22 and specifically for the slab adsorption environment in Figurec. The relative contributions of PPC and EPC to the total relaxation rate are relatively independent of the adsorption environment as shown in Figure S22 , but depend strongly on the vibrational mode. PPC primarily dominates the relaxation of FT and SA modes, whereas EPC predominantly influences the FR and IS modes.

Interestingly, the contributions from EPC and PPC do not fully add up to the total relaxation rate (see nonadditive component in Figurec) even when considering statistical uncertainties in the analysis. We observe a trend to underestimate the total rate, and the nonadditivity is more pronounced in modes where EPC dominates (FR and IS modes). EPC and PPC may not be fully independent, as shown in Figures S24–S27 . There, we find that PPC between the CO and the substrate vibrational modes increases if EPC is included in the equilibrium dynamics simulations. This result hints toward some level of nonadditivity and a deviation from eq when considering both effects simultaneously in the MDEF dynamics. Moreover, certain approximations taken in this work are expected to reduce non-additivity. For example, we have assumed a configuration-independent electronic friction tensor. In reality, the friction tensor is configuration dependent and changes as the geometry is displaced along different vibration modes, likely contributing further to nonadditivity.

Discussion

Vibrational Modes and Phonon–Phonon Coupling

Remarkably, all phononic lifetimes fall within the same order of magnitude. In contrast, electron–phonon lifetimes vary by 2 orders of magnitude (IS vs FT), highlighting a much stronger mode dependence compared to phonon–phonon interactions. Additionally, despite the significantly higher vibrational energy of the IS mode, its lifetime remains comparable to other modes, suggesting a surprising degree of energy independence in PPC. To understand this behavior, we analyze the VDOS and coupling density. In the slab and slab-adatom systems, energetically more localized substrate vibrational states lead to reduced PPC compared to the sharp and blunt tips. This is counteracted by a larger number of vibrational states per volume (approximately three times as many) in the slab and slab-adatom systems. The FT mode has an energy comparable to the lowest-energy substrate vibrational modes, while the FR mode is close in energy to the highest-energy substrate modes. Coupling of the molecular vibrations is most efficient to the substrate vibrations with similar energy, as shown in Figures and S24–S27. We can understand this in the following way: Energy conservation mandates that any newly excited vibrational mode arises from the creation or annihilation of N other modes, following $[eqn]$ . As the number of participating modes increases, the efficiency of the process decreases. Because the FT and FR modes overlap energetically with the substrate modes, they can access efficient N = 1 processes. However, individual substrate vibrations are strongly dependent on the surface shape. As a result, the phononic lifetimes of the FT and FR modes depend strongly on the surface shape, as shown in Table. The SA and IS modes are energetically relatively decoupled from the substrate modes compared to FT and FR. Therefore, their coupling of the substrate modes involves high-N processes. Figures S26 and S27 show that the SA and IS modes couple equally to all surface modes. Consequently, the influence of surface geometry on specific surface vibrations is effectively averaged out, as shown in Table.

Electron–Phonon Coupling

We find small electron–phonon lifetimes for the FR and IS modes and large lifetimes for the FT and SA modes, in good agreement with experimental and theoretical literature. ?,?−? ? For each adsorption environment, the electron–phonon lifetimes follow the trend of the vibrational frequencies of the respective modes, except for the FR mode. We can understand this behavior if we assume that the rate of EHP excitation is independent of the energies of the initial and final states. This corresponds to the wide-band, constant-coupling approximation, where one assumes a constant electronic DOS. In such a case, relaxation rates scale quadratically with the vibration energy. A similar explanation has been proposed by Persson.? As the vibrational energy increases, the available energy window for electronic transitions broadens, enabling more excitation pathways and thereby decreasing the lifetime. The vibrational energy of the IS mode is 2 orders of magnitude higher than that of the FT mode, consistent with the two-order-of-magnitude difference in their lifetimes. The FR mode deviates from this trend, with lifetimes comparable to the IS mode for the slab geometry. This can be explained by the large EPC reported by STM vibrational spectroscopy experiments for this mode.? Therefore, the EHP rate cannot be assumed to be independent of the involved energy states in this case.

Moving from the slab model to the sharp tip model, the lifetimes become larger. As Figure shows, the magnitude of the electronic PDOS becomes smaller with more corrugated adsorption environments. This leaves fewer electronic states available for electronic excitations. Consequently, EPC becomes less efficient, and the lifetime increases. However, the IS mode deviates from this trend. The available vibrational energy of the IS mode is approximately 270 meV, which is 1 order of magnitude larger than the energy of the SA and FR modes. The large vibrational energy of the IS mode increases the number of accessible electronic states for excitation. Even for the sharp tip, where the PDOS is highly localized, this energy range enables excitations between many more electronic states. Since the other modes carry significantly less energy, they are more strongly affected by the adsorption environment: The lifetime of the FT mode increases by 1 order of magnitude between the slab and the sharp tip, making it 2 orders of magnitude longer than that of the IS mode.

Combined Energy Dissipation Mechanisms

We have determined total relaxation rates using three methods: equilibrium correlation analysis, kinetic energy decay, and additive relaxation rates. Despite their significant conceptual difference, all three methods predict lifetimes that exhibit similar trends across the different systems and modes (see Figureb). This demonstrates the robustness of our analysis.

In Figurec, the sum of phonon–phonon relaxation rates from MD and electron–phonon rates from MDEF simulations tends to be smaller than the total relaxation rates obtained directly from MDEF. While the difference lies within the uncertainty and is therefore not statistically significant for most modes, it is nevertheless with the comparison of coupling densities with and without EPC: Coupling densities become larger if EPC is included in the simulations, indicating a weak nonadditive effect where electron–phonon-mediated PPC increases phononic relaxation. However, the approximations in our approach likely underestimate this nonadditivity. We assume a constant electronic friction tensor, whereas in reality, it is configuration-dependent. As the geometry is displaced along different vibration modes, the strength of EPC changes, likely enhancing the degree of nonadditivity.

The observation from Figurec that the relaxation rate of the FR (and IS) mode is dominated by EPC while the SA mode is not, is in good agreement with experiment. Asymmetric Fano line shapes have been reported for the FR mode, while symmetric Lorentzian line shapes have been measured for the SA mode.? Fano line shapes typically indicate that a discrete vibrational mode interacts with a continuum of electronic states, indicating significant EPC. Conversely, a Lorentzian line shape suggests that PPC dominates the relaxation of the vibration mode.?

However, it is important to note that the presence of EPC does not necessarily lead to a Fano line shape. Following Langreth,? the line shape of a vibrational mode in the presence of EHP excitations is governed by the dynamic dipole μ, with the asymmetry depending on the imaginary part of the equation:

$$[eqn]$$

Here μ_ions_ and μ_el_ are the ionic and electronic dipole, respectively. The electronic dipole consists of a real and an imaginary part. The imaginary part, $[eqn]$ , emerges due to the creation of EHPs as electrons and holes tunnel between the adsorbate and the substrate. Since this tunneling process is not instantaneous, the oscillations of the electronic dipole get out of phase with the driving electric field.

The imaginary part of eq becomes largegiving rise to an asymmetric line shapeif the imaginary electronic dipole $[eqn]$ is large compared to the ionic dipole μ_ions_. This occurs for modes with a dipole oriented parallel to the substratesuch as the FR modewhere screening by mirror charges significantly reduces the effective dipole strength. By contrast, modes with atomic motion perpendicular to the substrate, such as the SA and IS modes, exhibit a large ionic dipole. Subsequently, asymmetric line shapes would not be expected here, even if EPC is the dominant mechanism.

Our findings also inform the expected dynamic properties of CO on other coinage metals, as discussed in Section S4. The greater mass mismatch between CO and Ag or Au, compared to Cu, reduces PPC. Moreover, the phonon density of states shifts to lower energies from Cu to Ag to Au (Figure S23), leading to greater decoupling of the IS mode on Ag and Au than on Cu. Given the chemical similarity of the coinage metals, all of which are predominantly s-band conductors, this trend is expected to persist for different adsorption environments. Assuming comparable adsorption conditions, CO exhibits similar vibrational lifetimes on Cu, Ag, and Au. ?,? However, due to weaker interactions with Ag and Au, EPC is significantly reduced.? For adatoms or tip-like adsorption environments, trends observed on Cu likely extend to Ag and Au, although stronger interactions at low-coordination sites may facilitate shorter lifetimes relative to flat surfaces.

Perspective and Potential Experimental Validation

Finally, vibrational lifetimes calculated with equilibrium correlation analysis and kinetic energy decay are in good agreement, which demonstrates that both are suitable and equivalent approaches. Equilibrium correlation analysis can yield important information, such as coupling densities and phononic spectral functions. Our method for disentangling electron–phonon and phonon–phonon contributions to the vibrational dynamics of adsorbed molecules enables a better understanding of electronic and vibrational mass effects, as well as effects due to the local bonding environment. To connect theory with experiment, it is essential to identify techniques capable of validating mode-specific dissipation pathways. Most existing experiments for systems such as CO on Cu were conducted decades ago and are limited to flat surfaces. ?,? Data is lacking on low-energy modes, such as the FT. High resolution IR (infrared) and Raman spectroscopy allow for directly probing the dominant dissipation mechanisms via vibrational lineshapes: for vibrational modes with a small ionic dipole (such as those parallel to the surface where the dipole is screened), dominant EPC leads to characteristic Fano profiles (see eq), while PPC results in Lorentzian lineshapes.?

Tip-enhanced Raman spectroscopy (TERS) offers a powerful way to probe vibrational dynamics at the single-molecule level with atomic resolution, particularly when operated in the quantum tunneling regime of plasmonic enhancement.? By selectively tuning the excitation laser to match vibronic transitions, it becomes possible to achieve highly mode-specific resonance enhancement.? In STM-TERS, the application of a tunable bias voltage provides additional control: it can alter the tunneling gap and modulate the strength of EPC, which can be tracked via bias-dependent frequency shifts and intensity changes.? This could help to disentangle phonon–phonon and electron–phonon relaxation pathways and their contributions at specific adsorption sites.

To account for the adsorption site dependence (CO on different adsorption environments), experiments with spatial resolution are needed. STM-IETS (inelastic electron tunneling spectroscopy) can resolve low energy modes, is sensitive to IR and Raman inactive modes, and enables atomic spatial resolution.? Crucially, it is highly sensitive to EPC: Electrons lose energy by exciting vibrations during tunneling, and the resulting peak intensity scales with the electron–phonon coupling strength. ?,? Meanwhile, the peaks are broadened by both EPC and PPC. By combining IETS with a different method of vibrational excitation, such as IR, which excites modes via changes in the dipole moment, it becomes possible to disentangle the dissipation mechanisms. For instance, a mode showing a large IR peak but a small IETS signal would indicate that PPC dominates, while a strong IETS response would indicate that EPC is significant.

Finally, our work can improve the interpretability of vibrational line widths and energy dissipation observed in AFM, IETS, TERS, or IR spectroscopy. The spectral function calculated through equilibrium correlation analysis provides mode-specific line widths and lifetimes, while the separation of electron–phonon and phonon–phonon contributions enables clearer interpretation of lineshapes and mode-specific enhancements, particularly in IR and TERS. By disentangling dissipation pathways, our method can facilitate cross-correlation between IR and IETS signals and help assign adsorption sites based on the adsorption environment. Understanding which relaxation pathway dominates can enable experimental tuning. If the relaxation of a vibrational mode is primarily governed by EPC, adjusting electronic properties, such as the tip material, functionalization, or the bias voltage, can enhance signal detection and reveal details about the coupling. If PPC dominates, tuning temperature or using a different isotope (^14^C or ^13^C instead of ^12^C in CO, for instance) can shift vibrational frequencies and affect the coupling. In scanning probe experiments and single-molecule manipulation, where vibrational excitation can induce molecular motion or reactions, understanding electron–phonon and phonon–phonon relaxation supports design conditions that favor specific molecular transformations. Similarly, in catalysis, targeting specific vibrational modes can selectively trigger reactions, enhancing the control and efficiency of chemical processes.

Conclusion

The vibrational dynamics of CO-decorated copper tips that range from atomically flat to maximally sharp exhibit pronounced vibrational mode specificity, driven by the interplay of EPC and PPC, as well as significant geometric effects. Lifetimes due to EPC differ by two orders of magnitude between modes due to vibrational energies, the symmetry of the mode, and the DOS. Overall, CO vibrations of flatter tip geometries exhibit stronger electron–phonon relaxation, due to the increased energetic delocalization of the local electronic DOS. However, the high-energy IS mode does not show this trend due to the large energy window to excite electrons. Phonon–phonon lifetimes also decrease as the surface becomes flatter, although the IS mode does not follow this trend. The FT, FR, and SA modes are coupled to the substrate vibrations, leading to significant geometry dependence, while the IS mode is decoupled from the substrate vibrations. The sum of individual electron–phonon and phonon–phonon relaxation rates tends to underestimate the total relaxation rate computed by the same method. This weak nonadditivity arises from electron–phonon-mediated PPC, which enhances phonon–phonon interactions. Finally, the lifetimes obtained from equilibrium correlation analysis and kinetic energy decay are in good agreement between the two methods and also agree well with experiment.

Future work will explore the role of nondiagonal electronic friction on the mode selectivity and surface dependence. Incorporating a fully configuration-dependent electronic friction tensor will enable a more detailed investigation of nonadditivity.? Finally, extending the current framework to account for zero-point energy contributions will further improve the fidelity of atomic-scale friction models.

Methods

Modeling CO on Different Adsorption Environments

Electronic structure calculations are performed using DFT with the FHI-aims code,? employing periodic boundary conditions and a repeated-slab approach including a dipole correction. The sharp and blunt tip geometries are placed on a four-layer slab for better convergence. Convergence tests and further details are shown in Section S1.

Commonly used methods for first-principles surface simulations, such as Generalized Gradient Approximation (GGA) functionals, fail to accurately describe CO adsorption on Cu and other metal surfaces. ?−? ? Due to the many-electron self-interaction error, the energies of the states responsible for bonding the CO molecule to the Cu surface are misrepresented. This results in an incorrect preference for the fcc-hollow site over the top site, which experimental studies have shown to be the correct adsorption site.? The Random Phase Approximation (RPA) has been proposed as a method to yield qualitatively correct adsorption sites, while the PBE0 functional shows a slight preference for the correct adsorption site.? However, RPA calculations are unfeasible due to their high computational cost, while PBE0 is not ideal for modeling bulk metals since it tends to overestimate the delocalization of electrons, which can lead to an inaccurate description of metallic systems, particularly their electronic structure and density of states near the Fermi level. ?,?

In our work, we use the HSE06 screened hybrid functional,? which treats short-range interactions with exact exchange and reduces the amount of exact exchange at long-range. This helps mitigate some of the issues, such as overdelocalization of electrons and poor convergence in calculations.? We combine this function with the MBD-NL long-range dispersion correction?a combination that was previously shown to provide accurate adsorption structure and energetics.? With this combination, we predict an accurate PES (see Table S1). Geometry optimizations were performed for all systems with the four bottom layers frozen, prior to determining other properties (PDOS, EPC, PPC, etc.).

Interatomic Potential

We use the MACE equivariant neural network potential? and train it at the HSE06+MBD-NL level of theory, which accurately describes the potential energy surface of CO on Cu(111) as well as CO adsorbed on an adatom. Using a training set of 621 data points, the potential achieves an energy error of 1.6 meV/atom and a force error of 15.9 meV Å^–1^, which was determined using k-fold cross-validation, where we split our data set into five equally sized segments. The training data contains structures of all four adsorption environments (sharp tip, blunt tip, slab adatom, slab). Details regarding hyperparameters, training/test splits, and validation during dynamics simulations are shown in Section S2.

Vibration Calculations

Second-order force constants are calculated using finite displacements of 0.0025 Å. Energies and forces are determined using the MACE potential. During the vibration calculation, the bottom four layers of our systems are frozen. The second-order force constants are calculated for the top 10 mobile Cu layers as well as the CO molecule (approximately 500 atoms).

Description of Electron–Phonon Coupling

In first-order perturbation theory, we can express the (3 × 3) atomwise block of the electronic friction tensor Λ ij associated with atomic degrees of freedom R i and R j as follows: ?,?

$$[eqn]$$

The electronic friction tensor associated with energy ε = ℏω of a given vibration mode is calculated by summing over all possible excitations between effective single particle Kohn–Sham states ε _ kν _ and ε _ kν′ _. The indices ν and ν′ indicate the initial and final electronic state, respectively, and k is the associated wave vector. The terms in the sum are given by the nonadiabatic coupling matrix elements (first row of eq), the occupation factors f(ε _ kν _) and f(ε _ kν′ ) of the initial and final states (second row), and the conservation of energy (third row). Since we are dealing with electrons, the occupation factors follow the Fermi–Dirac distribution, which introduces the temperature dependence. The term ℏω accounts for the perturbing frequency of the vibration mode of interest. We calculate electronic friction with the electronic structure code FHI-aims.? Convergence tests for electronic friction calculations can be found in Section S1.6. We explicitly calculate electronic friction for C and O atoms. For the Cu atoms, we use isotropic electronic friction and a diagonal friction tensor. This friction value is calculated for a single atom in bulk copper (see Section S1.7). EPC calculations are performed using the PBE? exchange-correlation functional, as has been done in previous simulations of CO on Cu substrates. ?,? Details can be found in Section S1.7. We note that employing the HSE06 functional for electronic friction calculations in systems of this size is currently infeasible given available high-performance computing resources. Vibrational relaxation rates are determined from the friction tensor and the mass-weighted vibrational eigenmodes, ũα, at the relevant vibrational frequency ω_α:

$$[eqn]$$

Molecular Dynamics

We perform MD simulations using the atomic simulation environment (ASE).? A time step of 0.1 fs is employed. Since we are interested in the dynamics of a single CO molecule, we conduct MD simulations in an extended unit cell containing a single adsorbed CO molecule, shown in Figure. The extended unit cell facilitates the back-folding of substrate phonon bands, enabling effective sampling of the phonon band structure.

To account for EPC, we use MDEF, where the dynamics are governed by the Langevin equation:?

$$[eqn]$$

Here, V(R) represents the potential energy surface, which depends on the atomic positions R. The term Λ ij denotes a (3 × 3) atomwise block of the electronic friction tensor, while $[eqn]$ is the random force vector responsible for maintaining detailed balance at a given electronic temperature T, as dictated by the second fluctuation–dissipation theorem.? For the dynamics, the electronic friction tensor is evaluated in the harmonic approximation at zero perturbation frequency, based on a fixed nuclear geometry. A comparison of lifetimes with and without a perturbing frequency can be found in Section S3.1. During the calculations, the constant friction approximation is used for the electronic friction tensor, where the friction tensor values have been determined at the relaxed equilibrium geometry of the molecule on the surface. To determine relaxation rates and lifetimes from the MD and MDEF trajectories, we use three different methods:

Equilibrium Correlation Analysis

For each system, we first perform equilibration at 5 K, followed by computing a single, long MD or MDEF trajectory. To analyze the velocity time series, we use normal mode decomposition and determine the cross-spectrum using the Welch method. Thereby, the velocity time series is cut into multiple, partially overlapping segments. A window functiona Hann window in our caseis then applied to each segment before the cross-spectrum is calculated. The cross-spectra of the different segments are then averaged. This improves the smoothness of the cross-spectrum. The height and width of the peak(s) in the cross-spectrum are then determined by fitting a Lorentzian function. A detailed workflow is provided in Section S3.2.

Kinetic Energy Decay

For each system and vibrational mode, relaxation rates and lifetimes are extracted from ensemble averages over 64 MD or MDEF trajectories. We first equilibrate the system at 5 K. Then, we excite the vibrational mode of interest by a displacement amplitude equivalent to 50 meV along the respective eigenvector. Using normal mode decomposition, we project the velocity time series onto the basis of the vibration modes. We subsequently monitor the temporal decay of the envelope of the ensemble-averaged kinetic energy. Assuming exponential decayas would be the case for a damped harmonic oscillatorwe determine relaxation rates via an exponential fit. See Section S3.3 for full details.

Additive Relaxation Rates

Total relaxation rates are calculated by summing the individual relaxation rates associated with phonon–phonon and EPC. Phonon–phonon relaxation rates are derived from the inverse values listed in Table, which were determined through equilibrium correlation analysis. Electron–phonon relaxation rates are obtained from the inverse values in Table, which were directly calculated from the electronic friction tensor. See Section S3.4 for the computational workflow.

Supplementary Material