Vibrational energy transfer of excited CO molecules on NaCl(100): a non-equilibrium ab initio molecular dynamics study

Abstract

CO/NaCl(100) at monolayer coverage is a weakly bound adsorbate system for which unusual phenomena after infrared laser excitation of adsorbate vibrations have been observed. Among them are very long vibrational lifetimes in the CO stretch mode (∼4 ms), efficient vibrational energy pooling, and orientational isomerization from a “C-bound” to an “O-bound” configuration. Here we use non-equilibrium ab initio molecular dynamics (neq-AIMD), coupled with time correlation function (TCF) techniques, to study vibrational relaxation dynamics and related processes after selectively exciting specific normal modes (CO internal stretch, frustrated rotation) across submonolayer and monolayer coverages of this system. We are mostly concerned with picosecond dynamics, and in one case we use recently proposed machine-learned, AIMD-based potentials for CO/NaCl, to explore significantly longer timescales (∼ns). To monitor the dynamics, we present atom- and mode-resolved transient vibrational spectra, kinetic energies and qualitative time-dependent local mode couplings. Time-dependent vibrational sum frequency (VSF) spectra are reported for a flipping CO in the CO/NaCl(100) monolayer, following frustrated bending mode excitation. Important highlights of this work are as follows: (i) localization of vibrational energy in the CO stretch modes on ps timescales and signatures of vibrational energy pooling on ns timescales; (ii) sub- to few-picosecond lifetimes of frustrated bending modes of both “C-bound” and “O-bound” CO adsorbates resulting from two dissipative channels, namely, the surface phonons and neighbouring adsorbates causing successive flipping, desorption and diffusion events; and (iii) identification of vibrational non-adiabatic channels in metastable “O-bound” geometries following CO stretch excitation.

---

1 Introduction

Vibrational energy transfer is an ubiquitous phenomenon which mediates chemical reactions at gas/solid and gas/liquid interfaces. The study of vibrational energy pathways can shed light on adsorbate–adsorbate and adsorbate–surface inter- and intra-molecular couplings. This has implications for interface science in general. More specifically, the CO/NaCl(100) system to be studied here has also been mentioned as a model system for astrochemistry.¹ In the context of astrochemistry (see, e.g., ref. 2–4), vibrational energy relaxation is widely recognized as a necessary condition for the stabilization of newly formed molecules under low-temperature conditions, such as those encountered in dense molecular clouds and on interstellar dust grains, where excess formation energy must be efficiently dissipated into the surrounding condensed phase to prevent prompt desorption or dissociation. Microscopic mechanisms of vibrational energy transfer involving mode-specific vibrational excitation, desorption, diffusion and dissipation in astrochemically relevant systems have been thoroughly investigated, for instance, in ref. 5–7.

Monolayer CO on NaCl(100), in its most stable, low-temperature p(2 × 1) “tilted/antiparallel” (T/A) phase,⁸ is a condensed-phase system with weak surface–adsorbate interactions and long vibrational lifetimes of the CO internal stretch (IS) mode⁹ (∼4 ms), compared to CO on other surfaces, like CO/Pt(111)¹⁰ or CO/Si(100)¹¹ with vibrational lifetimes in the order of picoseconds and nanoseconds, respectively. This is due to a much slower relaxation to the predominant dissipation channel, i.e., the phonon bath, in the former case, due to a ten-fold frequency mismatch of the IS mode (∼2125 cm⁻¹, see below) with the Debye frequency of NaCl (∼220 cm⁻¹). The long lifetime leads not only to narrow linewidths in vibrational absorption spectra,¹² which allow one to resolve Davydov splittings (of about 6 cm⁻¹ for the T/A phase¹³), but it is also responsible for the observability of IR fluorescence signals,¹³ and a plethora of other interesting processes found for this system. One example is vibrational energy pooling (VEP) responsible for forming “base camps” of highly vibrationally excited CO molecules throughout the monolayer.^13,14 What is the origin of such a pooling process? Which molecules are prone to forming these base camps? Is it surface-mediated? To address such questions, Corcelli and Tully reported relevant vibrational energy flow pathways in this system and calculated rate constants for processes such as resonant and non-resonant excitation transfer within the CO monolayer, vibrational relaxation via phonons and vibrational dephasing.¹⁵ They demonstrated, using perturbation theory, that vibrational energy transfer between the CO molecules occurs much faster compared to vibrational relaxation to the surface phonons. The resonant V–V energy transfer (i.e., energy transfer between high-frequency stretch modes) was found to be the fastest, occurring on a ∼15 ps timescale, followed by non resonant V–V energy transfer on the sub-microsecond timescale, and finally vibrational population relaxation rates as well as single quantum vibrational energy loss via fluorescence on sub-second timescales. Localization of vibrational quanta was also found to occur via dephasing on a 10 ps timescale.

Furthermore, rotational isomerization (from C-bound to O-bound¹⁶ and vice versa¹) dynamics in CO/NaCl(100), succeeding VEP, makes the vibrational dynamics in this system even more interesting. One of the theoretical works in this direction is reported by Guo and coworkers,¹⁷ where the isomerizing CO molecules in highly excited vibrational states undergo “roaming” behaviour. That is, the CO molecule surpasses the isomerization saddle point, going farther away from the NaCl surface before transitioning to the O-bound configuration, thereby highlighting the importance of mode couplings in mediating such processes. The “O-bound” configuration, once created, is metastable, undergoing a back reaction even at low temperatures (∼20 K) to “C-down”.^1,18 Mechanistically, for this process, a thermally activated, phonon-assisted tunnelling-dominated mechanism involving the low-frequency frustrated rotational mode has been proposed.¹⁸ The low-frequency rotational mode has an energy of ∼150 cm⁻¹ for the “C-down” and ∼80 cm⁻¹ for the “O-down” configurations¹⁹ (see below) and, alongside the high-frequency CO stretch mode, will be a major focus of the present work.

In order to analyze vibrational dynamics, we shall resort to vibrational signatures and spectra, notably the vibrational density of states (VDOS) and vibrational sum frequency (VSF) spectra, and their time dependence. Similar information for both “C-bound” and “O-bound” configurations was recently already reported by us,¹⁹ which provided some insight into mode couplings, for example, of this system, albeit under equilibrium conditions. In order to identify vibrational energy distribution pathways and transient spectra and mode couplings within the system, we extend our study here to non-equilibrium situations. Specifically, time-resolved pump probe spectra will be calculated by probing the evolution of a vibrationally excited state by a pair of pulses with a delay time, over which the transient response of the system can be followed. Several examples showcase the power of this tool, both experimentally and theoretically, to study energy transfer in adsorbate–surface systems (see, e.g., ref. 20–23, and references therein). In this work, we utilize non-equilibrium ab initio molecular dynamics (neq-AIMD) simulations and a velocity-swapping method (in the absence of an explicit IR electric field) to pump energy into the desired vibrational mode and probe its fate over fixed spectral windows as a function of a delay time, Δt, signifying the action of a probe pulse. This approach allows us to disentangle the effects of neighbouring CO molecules and the NaCl phonons responsible for the vibrational energy distribution from the pre-excited CO molecule on a picosecond timescale. Adopting a recently developed machine-learned potential for CO/NaCl,²⁴ we will also go to nanosecond timescales in one example.

Our paper is organized as follows. Section 2 presents the structural models of CO/NaCl(100) as well the AIMD and TCF approaches used to obtain time-resolved vibrational spectra and other analysis tools used for post-processing the vibrational relaxation dynamics of the system. A brief description of the velocity-swapping approach in the context of the present system is also given. In Section 3, we present results pertaining to (i) relaxation dynamics of vibrationally excited (ν_stretch ≥ 1) CO stretch modes of CO(s) in pure “C-bound”, “O-bound” and mixed “C-bound”/“O-bound” CO(s) in monolayer coverages; (ii) relaxation dynamics of the same mode in pure “C-bound” submonolayer coverages; and (iii) relaxation dynamics of highly excited frustrated bending modes for an “O-bound” and “C-bound” CO in monolayer coverages. Time- and subsystem-resolved vibrational spectra and kinetic energy traces are reported for the above cases; additionally presenting an extension to the work of Melani et al.²⁵ from time-dependent VDOS to time-dependent VSF spectra. Other quantities, such as time-dependent distance correlation matrices used to highlight internal mode couplings in the system with reference to some of the above cases of pre-excitation, are also discussed. Finally, Section 4 concludes this work.

2 Models and methods

2.1 Adopted models and the electronic structure method

We mainly study four configurations, as shown in Fig. 1. The models are all based on supercells, with dimensions of 8.017 × 8.017 × 30 Å (∼20 Å vacuum along the z-axis, the direction perpendicular to the surface), if not stated otherwise. Specifically, we consider a submonolayer (Θ = 0.5) coverage comprising two CO molecules per cell (model A) in a tilted/antiparallel (T/A) phase; the others being monolayer coverages (Θ = 1), the T/A “all C-bound” form (model B), the most stable monolayer configuration at temperature below 35 K,^8,26 as well as metastable T/A forms with one (model C) or all four (model D) CO molecules within the unit cell inverted to the “O-bound” orientations. There, configurations were optimized using periodic DFT, employing the PBE exchange correlation functional²⁷ and a projector wave augmented (PAW) basis,²⁸ as implemented in the Vienna Ab initio Simulation Package (VASP).^29–31 Additionally, Grimme's D2 correction³² was applied to account for van der Waals interactions, with substitution of the van der Waals parameters for Na with those of Ne (to mimic Na⁺ (ref. 33)). Further details on optimization, energetics and frequencies obtained from normal mode analyses (NMA) can be found in ref. 19.

Fig. 1 Optimized (PBE+D2) structures for various adsorbate models for CO/NaCl(100) – see text. A is a half-coverage phase model (Θ = 0.5), and B–D are high-coverage phase models (Θ = 1). See also ref. 19.

2.2 Machine-learned potential for CO/NaCl(100)

In ref. 24, we have constructed various machine-learned potentials for CO/NaCl(100), adopting an equivariant, graph neural network (NN) approach, NequIP,³⁴ trained on AIMD trajectories (canonical and non-equilibrium ones – the latter are analyzed in detail and form a central part of the discussion in the present work). While introducing only small additional errors, the machine-learned NN potentials allow us in principle to treat arbitrary coverages and significantly longer propagation times in MD simulations – the latter being particular useful for long-lived CO stretch modes. Below, we shall use the so-called Generation 6 (Gen 6) version of the graph neural network potential for CO/NaCl proposed in ref. 24 to study the relaxation of model B described above (monolayer T/A, “all C-down”), in which the CO stretch modes of carbon monoxide were pre-excited over ns timescales (Section 3.1.3). (In the SI, Section S11, model C is also considered.)

2.3 Non-equilibrium ab initio molecular dynamics (AIMD)

In general, we run AIMD simulations starting from each of these models in three phases: (1) equilibrium, canonical phase (NVT ensemble, typically at T = 30 K), (2) equilibrium microcanonical phase (called eq-NVE) in what follows, and (3) a non-equilibrium microcanonical phase (neq-NVE), with the aim to compute non-equilibrium transient vibrational spectra and kinetic energies to obtain vibrational energy distribution pathways. Further details on the three phases can be found in SI Fig. S1 and in ref. 19. Monolayer configurations were found to be unstable at 300 K, displaying desorption, flipping and roaming of the CO adsorbates. However, no flipping of CO was observed for monolayer “O-bound” phases at 30 K in the NVT trajectories since the available energy is much lower (∼21 cm⁻¹) compared to the classical isomerization barrier (∼565 cm⁻¹ or 0.15 eV^19,24) in the monolayer system. From NMA performed in ref. 19, the characteristic frequencies of interest in this work are as follows:

• CO stretch in 1 ML “C-bound” CO/NaCl(100) ∼ 2125 cm⁻¹ (0.263 eV);

• CO stretch in 1 ML “O-bound” CO/NaCl(100) ∼ 2116 cm⁻¹ (0.262 eV);

• CO frustrated bending in 0.25 ML “C-bound” CO/NaCl(100) ∼ 150 cm⁻¹ (0.019 eV);

• CO frustrated bending in 0.25 ML “O-bound” CO/NaCl(100) ∼ 80 cm⁻¹ (0.010 eV).

These numbers reflect a red-shift of stretch and rotational modes upon “C-down” → “O-down” inversion, the former in agreement with experimental findings.¹³

Within our model, IR excitation energies corresponding to 0.26 eV are expected to excite both the “C-bound” and “O-bound” CO stretch modes by one vibrational quantum each. Additionally, energies in the same order exclusively pumped into the low-frequency bending/rotational modes can lead to highly excited frustrated vibrations, which can thereby induce flipping of the CO molecules (see below). Note that the 0.26 eV-excited bending mode(s) may also be thought of as the result of a (slow) transfer of a CO-stretch quantum into bending. Here, we aim to excite mainly the CO stretch mode and one of the frustrated bending modes of CO on NaCl(100), in one or multiple CO adsorbates per unit cell, by one or more vibrational quanta as specified below, neglecting any specific polarization axis with respect to the NaCl surface.

2.4 Vibrational excitation by velocity swapping

In order to mimic IR laser excitation of a particular vibrational mode such that the other vibrational modes are not directly coupled to the laser and remain at equilibrium, we sampled configurations from the equilibrium NVE AIMD trajectories at a given temperature (T = 30 K) and exclusively pumped excess energy into two different modes of the CO adsorbate(s) as mentioned. Specifically, we “pump” the CO stretch mode with energy equivalent to vibrational quanta, ν_stretch = 1, 2, 12, for one or more CO adsorbate molecules, or pump multiple quanta of energy equivalent to that of a CO-stretch quantum, 0.26 eV, into the bending mode of one of the CO adsorbates in the monolayer. In the latter case, for “O-bound”/“C-bound” CO, since one vibrational quantum of the frustrated bending mode is ∼80/150 cm⁻¹, respectively, 0.26 eV corresponds to pumping energy equivalent to about 26/14 quanta into the respective bending modes. Details of the velocity-swapping approaches utilized to pre-excite the CO stretch and frustrated bending modes can be found in the SI, Section S2.

2.5 Analysis of neq-AIMD trajectories

2.5.1 Instantaneous kinetic energy analysis. One quantity useful in probing energy redistribution within the system is the instantaneous kinetic energy, decomposed into different subsystem contributions, evaluated as , where N is the total number of movable atoms in the model system. In this case, the subsystems comprise individual CO adsorbates and the phonon bath provided by the NaCl surface. For energy flow from the pre-excited CO molecule (*CO), the respective (envelope of the rapidly oscillating) kinetic energy trace is fitted, for t ≥ t₀, as:where t₀ is the instant of time when pre-excitation occurs, τ₁ is the relaxation time of the respective pre-excited normal mode and is the limit of within the time frame of our non-equilibrium simulation phase. Relaxation rates, k, can be evaluated from the lifetimes via k = 1/τ₁. Similarly, kinetic energy traces for acceptance of energy by neighbouring subsystems are fitted via:

E^subsystem_kin(t − t₀) = −E^subsystem_a,kin(t₀)e^{−(t−t₀)/τ₂} + ΔE^subsystem_kin, (2)

where τ₂ is the corresponding vibrational excitation lifetime of the subsystem. Subsystems to be considered below are neighbouring CO molecules or the NaCl surface, respectively.

2.5.2 Time-dependent VDOS. In pump–probe spectroscopy, energy transfer to other degrees of freedom (DoF) or subsystems can manifest as bleaching of the IR intensity of the excited vibrational mode, and consequential spreading of the spectral progressions (so called “spectral diffusion”).³⁵ Analogous to such pump–probe techniques, we evaluate time-dependent vibrational density of states, spanning fixed time windows over the timescale of neq-AIMD trajectories, for different pre-excitation scenarios. This corresponds to the so-called “non mass-weighted” counterpart of instantaneous kinetic energy spectral densities (KESDs), given as:where T is the time window over which the velocity–velocity autocorrelation function (VVAF) is computed at time t. N_t is the total number of timesteps sampled within [t − T, t] and N is again the total number of moving atoms. Furthermore, t_j = (j − 1)δt, with the timestep δt = 0.5 fs. 〈⋯〉 indicates an ensemble average. The VVAF in eqn (3) for the total system is normalized for each time window T, chosen in the range between 1 and 2 ps. Additionally, a Hann window function w(t) = cos²(πt/2T) is multiplied with the correlation function before the Fourier cosine transform.

Both Cartesian and internal VVAFs are evaluated, as in ref. 19. While the former is decomposed into molecular subsystems, the latter are decomposed into internal DoFs, which then yield molecule-/mode-resolved transient VDOS curves. Note that although the VVCF (velocity–velocity correlation function) approach to vibrational spectra is rigorously justified only for equilibrium scenarios, its reliability for appropriate time windows has been demonstrated earlier for calculation of physical properties in non-equilibrium scenarios.^25,35 However, we prioritize statistical exponential fitting of subsystem energies for the calculation of vibrational lifetimes throughout the manuscript; a comparison in some cases to lifetimes obtained via fitting VDOS decay has been discussed for completeness.

A more distinct observation of the modulated intensities and shifts in the spectral profiles as a result of vibrational relaxation is provided by the integrated VDOS difference spectrum, as reported earlier in ref. 25, given as:

where the delay time (Δt) represents the time delay between the pump and probe. The integration limits are over the entire frequency range (0,∞), or a frequency window [ω₁,ω₂] for partial ΔVDOS curves evaluated for specific fingerprint regions, such as the ν_stretch = 1 region. A rise/dip in integrated VDOS difference curves would indicate flow into/from the corresponding subsystem/internal mode, allowing for a clearer representation of the changes in vibrational dynamics resulting from the redistribution of excess energy. Details regarding the normalization of the time-dependent VDOS can be found in the SI, Section S3. (Note that Δt = 0 corresponds to the point of pre-excitation; therefore our ΔVDOS (Δt = 0) curves always start from zero.)

2.5.3 Time-dependent VSF response. We also extend the neq-AIMD/TCF approach for computing transient VDOS towards modelling time-resolved VSF spectra based on the formalism of Ohto et al.³⁶ for one of the interesting scenarios: flipping of “O-bound” CO in monolayer CO/NaCl(100) (model C), following frustrated bending mode excitation (Section 3.3.2). We are particularly interested in the xxz-polarized VSF response, governed by the susceptibility tensor element χ_xxz⁽²⁾, where z is the direction perpendicular to the surface, and x is parallel to the NaCl surface. The latter is computed within an ssVACF (surface-specific velocity autocorrelation function) approach³⁶ as

Here, c′ = −0.889e, a′ = 2.99a₀², b′ = 1.55 rad⁻¹, e′ = −2.146 rad and d′ = 5.62a₀² are parameters obtained from fitting calculated dipole and polarizability derivatives, respectively, as detailed in ref. 19. (e and a₀ are the elementary charge and the Bohr radius, respectively.) v^CO_z,i and v̲^CO_i are the z-projected bond velocity and bond-projected velocity vectors of the ith CO molecule, r̲^CO_i is the bond distance vector and r̂^CO_z,i(t_j) is the instantaneous bond unit vector component projected along the z direction. The ensemble-averaged terms 〈⋯〉 (except the orientational terms) are the aforementioned ssVACFs, which yield zero signal for centrosymmetric molecules and, in this case, for molecules aligned parallel to the surface. Additionally, T′ is the time window over which the ssVACF is computed at time t. Every reported VSF is normalized for the given time interval T′ = 2 ps. Again, t_j = (j − 1)δt, where δt = 0.5 fs in this work. The quantum correction factor Q(ω,T′;T_avg) in the above expression, which is normally given as Q(ω) = βħω/(1 − exp(−βħω)), with β = 1/k_BT and T being the temperature, now also depends on the chosen time window T′. Therefore, an average temperature, T_avg, is obtained over each time window, T′, which is also ensemble averaged, and then used for Q(ω). Below, we shall plot the square of the total response (|χ_xxz⁽²⁾(ω,t;T_avg)|² = Re(χ_xxz⁽²⁾)² + Im(χ_xxz⁽²⁾)², as well as the real (Re(χ_xxz⁽²⁾)) and imaginary responses (Im(χ_xxz⁽²⁾)). |χ_xxz⁽²⁾|² is a measure for the total VSF spectrum and Im(χ_xxz⁽²⁾) is expected to yield orientational information about the adsorbates.

2.5.4 Distance correlation matrices and distribution probabilities. Another quantity that can aid in identifying transient mode couplings between internal DoFs of CO adsorbates during neq-AIMD simulations is the time-dependent distance correlation matrix between the internal DoFs, a and b, of two individual CO adsorbates (inter) or for the same molecule (intra), evaluated over fixed-time windows analogous to the transient VDOS and VSF spectra. Following ref. 19, the difference distance correlations are computed as:

ΔR(a,b;Δt) = R(a,b;Δt) − R(a,b;Δt = 0), ∀a,b ∈ {r,θ,ϕ,X,Y,Z}, (6)

where R(a,b;Δt) is the distance correlation over a time window Δt, and R(a,b;Δt = 0) is the distance correlation just before the pre-excitation of the adsorbate(s) over a time window of 2 ps. The definition of the distance correlation matrix can be found in ref. 19. Note that the difference distance correlation matrices are of dimension 6 × 6 (owing to the six internal DoFs of each CO) and are evaluated over fixed time windows of 2 ps. This quantity is also ensemble-averaged, similar to the VVAFs in the evaluation of spectra.

Quantities like the angular distribution probabilities, P(θ), P(ϕ), and the joint angular distribution probabilities, P(θ,ϕ), which are defined in ref. 19 and in the SI, Section S4, are also reported occasionally for some of the studied scenarios.

3 Results and discussion

Except for the long-time runs using NN potentials in Section 3.1.3, all results discussed below follow from the analysis of the non-equilibrium AIMD trajectories, 5 ps before and 10–20 ps after excitation of the corresponding vibrational modes. Results shown are averaged over three trajectories if not stated otherwise. We mainly discuss the following scenarios: (i) vibrational relaxation of the CO internal stretch of one or more “C-bound”/“O-bound” CO adsorbates in monolayer and submonolayer coverages over the NaCl(100) surface, (ii) vibrational relaxation from highly excited vibrational levels (ν_stretch > 10) of the CO stretch mode of a single “C-bound”/“O-bound” CO in the monolayer, and (iii) vibrational relaxation of a frustrated bending mode of a “C-bound”/“O-bound” CO in a monolayer coverage.

3.1 Vibrational relaxation of the CO stretch mode in monolayer CO/NaCl(100)

3.1.1 Vibrational relaxation from low vibrational states.
3.1.1.1 Monolayers with “all C-down” orientation. We first consider the post-excitation picosecond dynamics of monolayer T/A “all C-down” (model B in Fig. 1) for two cases where we pre-excite the internal stretch mode(s) of (i) one CO by one vibrational quantum (0.26 eV), and (ii) all COs in the monolayer by one quantum each (1.04 eV). For both, we find an energy gain equivalent to the amount of vibrational quanta pumped into the system.

Fig. 2 displays the instantaneous kinetic energies of the total system as well of the individual subsystems (CO molecules and NaCl surface) for both eq-NVE and neq-NVE phases, for both cases (i) and (ii). We find for both, after pre-excitation, there is a periodic to and fro transfer of the excess kinetic energy into the potential energy of the system and a neglible amount of energy transfer from the adsorbates to the NaCl surface.

Fig. 2 Kinetic energy analysis of the non-equilibrium NVE-AIMD simulations starting from model B, averaged over three trajectories, in the cases where (a) one CO molecule is pre-excited by one vibrational quantum (0.26 eV) and (b) all four CO molecules are pre-excited by one vibrational quantum each (1.04 eV). The right panels in both cases showcase the averaged kinetic energies of the subsystems: the individual CO molecules and the NaCl surface. Note that the CO(s) were perturbed after the initial 5 ps of phase (3) of the neq-AIMD approach, as shown in SI Fig. S1.

For the case of excitation of a single CO adsorbate by one quantum (case (i), Fig. 2(a)), the energy stays more or less localized on the pre-excited molecule, ¹CO, with small amounts of energy exchange with the CO adsorbates, ²CO–⁴CO. This is also manifested as a rise and dip in the integrated ΔVDOS curves, as shown in Fig. 3(a). Energy inflow results in a redistribution of the VDOS of the total system, causing an enhanced intensity in the pre-excited molecule VDOS and a mitigated intensity for all other subsystems initially. The molecular VDOS of the pre-excited CO (not shown) shows a red-shifted peak (Δω ∼ 40 cm⁻¹) with respect to the equilibrium CO-stretch frequency due to vibrational heating of the molecule leading to mode softening, which then blue-shifts over time, indicating energy outflow from this molecule. Consequently, at certain Δt values, we see spectral broadening of the CO stretch peaks and the appearance of similar red-shifted shoulder peaks for the neighbouring adsorbates as well.

Fig. 3 Integrated change in the decomposed Cartesian VDOS curves, ΔVDOS, cf. eqn (4) with integration limits (0,∞), for model B when (a) one “C-bound” CO in a ML CO/NaCl(100) system is pre-excited by one vibrational quantum (0.26 eV); (b) all four “C-bound” CO molecules are pre-excited by one quantum (1.04 eV) along the corresponding CO stretch modes.

Now coming to the case where we excite all four CO adsorbates per cell in the monolayer by one vibrational stretch quantum each (case (ii), Fig. 2(b) and 3(b)), we still find a negligible transfer of energy to the NaCl surface, with significant energy exchange only among the CO adsorbates. The extent of excitation in each CO adsorbate can be determined either from the change in magnitude of their integrated ΔVDOS curves (cf. Fig. 3(b)), or the intensity of the respective red-shifted peak observed in the CO stretch region (cf. SI Section S5 and Fig. S3). SI Section S5 also provides further details for case (ii), among them, a slight population of the ν_stretch = 2 in the VDOS curves, indicating an onset of VEP.

Summarizing, on the ps timescales shown, we have no significant energy transfer from CO stretch to surface phonon modes. Even on nanosecond timescales, there is no such energy transfer (see Section 3.1.3). This is in agreement with the experimental long (∼ms) lifetimes mentioned above. However, there is rapid exchange of energy between the adsorbed CO molecules. In fact, by computing the ΔVDOS curves according to eqn (4) in a characteristic energy window [ω₁, ω₂] = [2050, 2200] cm⁻¹ centered around the CO stretch frequency only, we find a near-exponential post-excitation decay of ΔVDOS curves for individual molecules, giving a lifetime of a localized CO stretch vibration of 3.0 ps due to intermolecular vibrational energy transfer, for the case when all four molecules were excited by one quantum, cf. SI Section S5 and Fig. S4 there. Note that from the Davydov splitting Δ_Dav ∼ 6 cm⁻¹ for T/A 1 ML CO/NaCl(100) mentioned earlier, one expects a Rabi oscillation time of τ = 2πħ/Δ_Dav ∼ 5.5 ps, which is of the same order of magnitude. This timescale is also consistent with the resonant V–V energy transfer (within ∼15 ps) between two CO molecules adsorbed on NaCl(100) computed by Corcelli and Tully,³⁷ albeit for the ¹³C¹⁶O/NaCl(100) system.

Similarly, using an energy window characteristic for librational modes (ϕ, θ), one finds that these modes are energy-accepting, and a characteristic “lifetime” of 1.23 ps is found for the latter; see also Fig. S4. Note that the above lifetimes are extracted from VDOS decay curves and should therefore be regarded as preliminary estimates, as they may depend on the finite time window chosen for the transient VDOS analysis.

It is interesting to compare the above observations with NEMD simulations of Ferrari et al. in ref. 38 to study vibrational energy transfer processes in crystalline and amorphous CO clusters. In their study, the CO stretch mode was found to couple predominantly to the CO stretch vibrations of neighboring CO molecules, with essentially no direct energy transfer to low-frequency rotational or translational modes, highlighting the inefficiency of direct multiphonon dissipation pathways, in concordance with our results. The V–V transfer (resonant and non-resonant) was found to occur, albeit on much longer timescales from sub-nanoseconds to a few nanoseconds, respectively.

3.1.1.2 Monolayers with “O-down” CO orientations. Next, we examine how the orientation of the CO adsorbates (“O-down” or “C-down”) influences VET dynamics in the monolayer CO/NaCl(100) system. A pre-excited ν_stretch = 1 “O-down” CO in a one-“O-bound” monolayer (model C, Fig. 4(a)) follows similar dynamics as the “C-down” CO in a monolayer discussed above, at least when the negligible energy transfer to the surface is concerned. However, closer inspection shows there is much less VET to the neighbouring CO(s) than before due to the red-shifted stretch mode of the inverted CO.

Fig. 4 Kinetic energy analysis of the neq-NVE-AIMD simulations, averaged over three trajectories in the case where one “O-bound” CO molecule is pre-excited by one vibrational quantum (0.26 V) in the presence of “C-bound” neighbouring adsorbates ((a), corresponding to model C), and (b) in the presence of neighbouring “O-bound” adsorbates (corresponding to model D). The left and right panels showcase the averaged kinetic energies of the two subsystems—CO molecules and the NaCl surface—and the averaged kinetic energies of the individual CO adsorbate molecules, respectively. Note: ¹CO was perturbed after the initial 5 ps of the NVE run.

Next, energy equivalent to one vibrational quantum (0.26 eV) is pumped into the CO stretch mode of an “O-bound” CO in an all-“O-bound” monolayer configuration (model D), Fig. 4(b). We note that this configuration is a metastable minimum in itself and even at equilibrium conditions, during the last phase of the eq-NVE simulations, adopts a configuration where one of the “O-bound” CO molecules becomes nearly parallel to the NaCl surface. Inspecting non-equilibrium conditions after pre-exciting one of the tilted “O-bound” CO adsorbates, the instantaneous kinetic energies (KE) in Fig. 4(b) illustrate a rapid outflow of energy from the pre-excited molecule during the initial 5 ps, beyond which the fluctuations in the KE stabilize. Approximately 30% of the excess energy undergoes decay, primarily channeling to the NaCl surface and neighbouring CO adsorbates. After about 3 ps of the CO excitation, the surface phonons remain vibrationally “hot”. A complementary perspective is provided by the integrated ΔVDOS curves reported in Fig. 5(b), indicating a stronger dynamic involvement of the NaCl surface in this case. A further glance into the molecule-resolved, integrated, but internal-mode-resolved ΔVDOS reported in SI Section S6 and Fig. S6 illustrates a vibrational energy flux between the CO IS mode and the frustrated librational modes (ϕ, θ).

Fig. 5 Integrated change in Cartesian VDOS curves, ΔVDOS, cf. eqn (4) with integration limits (0, ∞), in the case where one “O-bound” CO molecule is pre-excited by one vibrational quantum (0.26 V) (a) in the presence of “C-bound” neighbouring adsorbates (model C), and (b) in the presence of neighbouring “O-bound” adsorbates (model D).

Since, in the case of model D, excitation of a single CO stretch mode leads to a clear time evolution in the kinetic energies E_kin(t) of various subsystems, we can perform a single exponential fit of E_kin(t) for the excited molecule, ¹CO, according to eqn (1); cf. Fig. S7(a) in the SI. There it is demonstrated that the fit gives a de-excitation lifetime τ₁ = 0.89 ps with residual excess energy ΔE_kin of 0.08 eV for the relaxing ¹CO. (Fitting the ΔVDOS for the pre-excited molecule given in Fig. 5(b) yields a similar effective lifetime of ∼0.95 ps.) Using the fit function (2) for E_kin(t) of the accepting subsystems reveals that the surface, NaCl, accepts energy with an excitation lifetime of τ₂(surf) ∼ 1.3 ps, and the neighbouring (unexcited) CO molecules do so with a lifetime of τ₂(rest CO) ∼ 3.7 ps; cf. SI Fig. S7. Note also that 1/τ₁ ∼ 1/τ₂(surf) + 1/τ₂(rest CO) (one finds 1/(1/1.3 + 1/3.7) ps = 0.96 ps). Given our limited sampling, and following the methodology of ref. 38, we performed an individual analysis of the nonequilibrium trajectories by fitting the time-dependent energies, E_kin(t), over different time windows. This procedure yields a range of de-excitation and excitation lifetimes for the subsystems discussed above. We obtain de-excitation lifetimes of [0.2, 1.5] ps for the initially pre-excited CO, with flipping observed in 33% of the trajectories, and excitation lifetimes of [1.4, 1.8] ps and [0.7, 3.3] ps for neighboring CO molecules and the surface, respectively.

In summary, for low vibrational excitation energies, “O-down” excitation leads to similar dynamics as for “C-down” if the neighbours are still “C-down” (albeit with slower CO–CO VET), whereas if the neighbours are also “O-down”, another relaxation channel (to the surface) may open on the ps timescale, in addition to VET to CO neighbours. At higher vibrational excitation energies (>3 eV), as discussed in the next section, “O-down” configurations may be transiently trapped in metastable minima geometries similar to those found in gas phase CO within the initial few picoseconds, unlike all “C-down” configurations where only localization of energy on the pre-excited molecule is observed. Orientation of the CO molecules can therefore have a strong impact on isomerization rates as well as on relaxation pathways depending on the energetic landscape around the isomerizing CO. We are aware, however, that our predictions are based on methods and sampling that need to be improved in the future to make fully quantitative statements for comparison with experiments to come.

In the SI, Section S10, we look at neq-NVE-AIMD trajectories more closely. For the above-studied case of an internal CO stretch mode of an excited “O-bound” adsorbate in an all “O-bound” monolayer, we observe, as a result of the VET, subsequent flipping of the pre-excited molecule, followed by flipping of neighbouring adsorbates. Flipping events are found to occur predominantly via roaming and lateral displacements of the CO molecule, with instances of vibrationally induced desorption (cf. SI Fig. S14 in one of the initially unexcited molecules [which still weakly interacts with the surface and could adsorb back on it]).

Further analysis of time-dependent distance correlation matrix elements ΔR(a,b) as defined in eqn (6) sheds light into the transient mode couplings and illustrates possible intra- and inter-molecular vibrational non-adiabatic channels of energy transport (see the SI, Section S9 and Fig. S12), with energy flux between the high-frequency (HF) CO IS mode and low-frequency (LF) librational modes.

3.1.2 Vibrational relaxation from higher vibrational states. Next, we studied a scenario where we pump energy equivalent to ∼12 CO stretch quanta (3.12 eV) into an “O-bound” CO in a one-“O-bound” monolayer (model C) configuration; cf. Fig. 6. In the figure, the excitation time is again at 5 ps.

Fig. 6 Model C: (a) averaged instantaneous Z values (distance of the CO molecules from the surface); (b) instantaneous bond lengths; (c and e) joint angular probabilities, P(θ,ϕ), averaged over three trajectories, (c) before and (e) after pre-excitation of the “O-bound” CO on NaCl with twelve vibrational quanta. (d) A metastable configuration in which the system is trapped after pre-excitation.

From Fig. 6(b), one first of all notes that, in this high-energy scenario, significant lengthening and large fluctuations of the CO bond length r of the excited (inverted, “O-down”) molecule occur, while the neighbouring (“C-down”) molecules are not as strongly affected on the ps timescale shown. One thus observes a localization of energy in the pre-excited molecule and large oscillations in the time-dependent kinetic energy curves, also as a consequence of large-amplitude motion. Inspecting the joint angular probability distributions P(θ,ϕ), as mentioned in Section 2.5.4 (Fig. 6(c) and (e)), and CO–surface distances Z (cf. Fig. 6(a)) reveals another interesting finding: the pre-excited CO molecule gets somewhat “trapped” in a metastable minimum. The mean value of the dynamic tilt angle of the excited CO shifts from 140° to 90° with a decrease in azimuthal disorder. The metastable configuration looks similar to what could be a saddle point during the transition from the “O-bound” to “C-bound” isomer (cf. Fig. 6(d)). This scenario is reminiscent of the energy transfer dynamics observed in highly vibrationally excited CO molecules in the gas phase studied by Guo and coworkers,³⁹ who found new minima structure(s) of CO molecules emerging when one of the CO molecule(s) was in a highly excited vibrational state (ν_stretch = 10) and the other was in the vibrational ground state; they attributed it to the formation of collisional complexes owing to enhanced dipole moments at higher vibrational excitation energies. Such metastable structures observed here could, therefore, be important in the isomerization dynamics for such highly excited CO molecules on the monolayer of NaCl(100).

3.1.3 Vibrational relaxation of the CO stretch modes: nanosecond dynamics. As mentioned, we also considered, for the 1 ML T/A CO/NaCl(100) model B, long(er)-time dynamics using the machine-learned potential Gen 6 of ref. 24, constructed from the short-timescale AIMD simulations of this work. In Fig. 7, we show kinetic energies of the surface and individual CO molecules on a ns timescale (left panels), and on shorter, ps timescales (right), in the case where initially all four CO molecules were excited by one ν_stretch quantum of energy each. This corresponds to the situation shown in Fig. 2(b).

Fig. 7 Model B, single-trajectory long-time (ns) propagation after exciting every (“C-down”) molecule with a vibrational quantum ν_stretch each (1.04 eV). The left panels show surface- and individual-CO-resolved kinetic energies on the timescale of 2.4 ns, while the right panels show the corresponding short-time (30 ps) dynamics. On the left, in the fourth panel, the energies of one (0.26 eV) up to four (1.04 eV) vibrational quanta are indicated as dashed horizontal lines.

The striking result we observe here is that, particularly on longer timescales of hundreds of ps to several ns, individual CO molecules can pool vibrational energies (left panels of Fig. 7). For example, molecule ²CO has acquired a vibrational energy of about three quanta after around 1 ns, whereas the other three CO molecules lost a significant portion of their initial vibrational energy. As the dashed lines in the ²CO panel (left) indicate, one may even anticipate seeing some vibrational energy quantization in the accepting molecule. No significant vibrational energy transfer to the surface is observed, even on the ns timescale, as expected. While the conclusions that can be drawn from single trajectories are limited, Fig. 7 indicates the occurrence of VEP from classical MD on a ns timescale. Note that the extent of pooling is limited by the finite simulation cell size, as only four vibrational quanta are available from the four CO molecules adsorbed per unit cell in our simulations. Making the simulation cell larger would therefore relax the finite-size bottleneck on VEP by increasing the number of available vibrational quanta and enabling intermolecular coupling pathways beyond nearest-neighbour COs, thereby potentially altering both the extent and timescale of pooling.

The NN potentials used here offer a way to a more systematic study, on extended spatial domains and even longer propagation times, to consider this and related phenomena in a time- and atom-resolved fashion, while still based solely on ab initio information.

3.2 Vibrational relaxation of the CO stretch mode in submonolayer CO/NaCl(100)

We now study the relaxation dynamics of a pre-excited “C-bound” CO in a submonolayer coverage (0.5 ML, model A in Fig. 1). In this case, VET occurs from the pre-excited CO molecule (termed as ²CO here) to the neighbouring CO molecule predominantly after 8 ps of the simulations, on average. Fig. 8(a) illustrates the molecule-resolved Cartesian ΔVDOS. In addition, mode-resolved ΔVDOS for the different CO adsorbates are reported in Fig. 8(b)–(d). The initial transfer of energy from the CO stretch mode of ²CO to ¹CO is slow, followed by a gradual increase in the VET. Some excitation also transfers to the frustrated azimuthal mode (ϕ) of both adsorbates. Inspecting the integrated ΔVDOS curves (Fig. 8(c) and (d)), we find that the vibrational energy flows into the angular DoF (ϕ) of the neighbouring CO molecule between 8–10 ps. This energy flows minorly from the CO IS mode of the same molecule and majorly from the CO stretch mode of the pre-excited molecule. The IS mode of the neighbouring molecule also simultaneously accepts energy. There is also a very small amount of energy transfer between the other internal DoFs such as the tilt angle (θ) and perpendicular CO–surface mode (Z) of the unexcited molecule. Again, the surface phonons do not exchange energy with the adsorbates, similar to the all-“C-bound” monolayer case.

Fig. 8 Integrated change in (a) Cartesian VDOS over a time window of 1 ps and (b) internal VDOS over a time window of 2 ps for 0.5 ML CO/NaCl(100), where the CO stretch of one CO molecule (²CO) is pre-excited with one vibrational quantum. Decomposed integrated change in internal VDOS for (c) ¹CO and (d) ²CO. Δt = 0 indicates the time of pre-excitation.

Time-dependent difference distance correlation matrices evaluated for this case (Fig. S13 in the SI) also hint at vibrational non-adiabatic energy redistribution channels similar to the earlier reported scenario of a relaxing “O-bound” CO in an all-“O-bound” monolayer, albeit with no CO flipping observed for the present case within the few-picoseconds timescale.

3.3 Vibrational relaxation of a frustrated bending mode in monolayer CO/NaCl

3.3.1 Vibrational relaxation dynamics. So far, additional energy was pumped exclusively into the CO IS mode for different “C-down” and “O-down” CO molecules adsorbed on a NaCl(100) surface. Despite the excitation energy always exceeding the barrier for “C-down” ↔ “O-down” inversion, only in one such case (all-“O-bound”-CO monolayer on NaCl(100), model D) was flipping of the CO adsorbates observed. Now we study direct pre-excitation of the frustrated rotational mode, which induces flipping dynamics as we shall see. For this purpose, we performed a direct excitation of a frustrated rotational mode of one of the CO adsorbates in the monolayer system, and follow the energy flux in the system. High-lying vibrational states of a frustrated rotational mode of a “C-bound” or “O-bound” CO surrounded by “C-bound” adsorbates are excited (models B and C, respectively), such that the injected energy overlaps with one vibrational quantum of the CO internal stretch mode (0.26 eV). This situation can also be seen as a model when a decaying CO stretch vibration has already transferred all its energy into frustrated rotations. Results pertaining to both scenarios, exciting “C-bound” or “O-bound” molecules, are discussed below. Note that the results are ensemble-averaged over six trajectories due to two different pre-excitation approaches being utilized (see SI Section S2.2 for further details).

When a frustrated bending mode of an “O-bound” CO in a one-“O-bound” ML configuration (model C) is pre-excited with 0.26 eV, corresponding to ∼26 quanta, the kinetic energy curves in Fig. 9(a) show significant energy transfer to the NaCl surface within the first 5 ps, beyond which the phonons are in a “vibrationally hot” state. This results from the efficient coupling of the vibrationally excited frustrated bending modes with the NaCl phonons. Consecutive flipping of the pre-excited CO molecule in all trajectories is observed in the initial ∼3 ps, accompanied by roaming dynamics: the CO molecule undergoes precessing motion, flipping to the other isomer while moving away from the surface with lateral translation, and eventually returning nearly on top of the same Na⁺ site; cf. Fig. S15(a) in the SI, Section S10.

Fig. 9 Kinetic energy analysis of the neq-NVE-AIMD simulations, averaged over six trajectories for the cases where (a) a frustrated rotational mode of an “O-bound” CO molecule (¹CO) is pre-excited by 26 vibrational quanta (0.26 eV, model C), and (b) a frustrated rotational mode of a “C-bound” molecule (¹CO) is pre-excited by 14 vibrational quanta (0.26 eV) in the presence of “C-bound” neighbouring adsorbates (model B). The top and bottom panels showcase averaged kinetic energies of two subsystems—CO molecules and the NaCl surface—and the averaged kinetic energies of the individual CO adsorbate molecules, respectively. Note that pre-excitation occurs at 5 ps.

The VET dynamics can be further followed from the time-dependent VDOS curves reported in Fig. 10. One can see a modulation in the phonon VDOS; for instance, the emergence and shift of the shoulder peak around 150–180 cm⁻¹. In Fig. 10(a), a low-intensity broad peak is observed for the pre-excited molecule in the CO stretch region around 2100 cm⁻¹ at Δt = 2 ps, spanning over 200 cm⁻¹, reflecting spectral diffusion. One can discern here the excitation of the higher vibrational states of the librational modes, which overlap and subsequently transfer excitation to the CO stretch mode, further affirmed by the integrated (internal) ΔVDOS curves reported in SI Fig. S8. At later delay times, the reduction in the spectral width of the CO stretch peak can be attributed to the loss of excitation of the frustrated bending modes, leading to a narrower and well-defined spectral profile, now arising from the disentangled r mode. Further microscopic insights obtained via the internal ΔVDOS can be found in SI Section S7.

Fig. 10 Time-dependent decomposed Cartesian VDOS, averaged over six trajectories of 20 ps each, for (a) the pre-excited CO molecule, ¹CO, as well as (b) the NaCl surface after pre-excitation of a bending mode of one “O-bound” CO in a one-“O-bound” monolayer CO/NaCl(100) by ∼26 vibrational quanta (0.26 eV, model C). Here Δt = 0 refers to the time (t = 5 ps) at which the pre-excitation is done.

We also performed an exponential fit to the instantaneous kinetic energy (cf. SI Fig. S10) of the pre-excited molecule. Considering the statistical spread of the fit, this yields an effective relaxation lifetime in the range of [0.5, 4] ps for the frustrated bending mode. A single exponential fit to the integrated ΔVDOS of the pre-excited CO gave an effective relaxation lifetime of 0.4 ps, which is smaller compared with that obtained from the instantaneous energy fitting. Again, for the energy-accepting subsystems, neighboring CO and the NaCl surface, energies can be fitted in individual trajectories over different time segments, yielding [0.5, 7] ps and [0.5, 3] ps, respectively. See also SI Fig. S10, where average lifetimes are reported.

In addition, note that we also utilized the recently developed NN Gen 6 potential²⁴ to simulate the present scenario for an extended time period of 200 ps and found similar dynamics of the relaxing “O-down” CO, albeit with a de-excitation lifetime of 2.6 ps. More details can be found in SI Section S11.

Similar observations are found when pre-exciting a frustrated rotational mode of a “C-bound” CO in an all-“C-bound” monolayer (model B) by (about) 14 vibrational quanta, corresponding also to 0.26 eV excitation energy. The instantaneous kinetic energy curves reported in Fig. 9(b) illustrate a loss of nearly 80% of the excess kinetic energy pumped into the CO molecule within the first ps of the simulation. In this scenario, the excess kinetic energy is primarily accepted by the NaCl surface, with some energy transfer to neighbouring CO molecules, similar to the previous scenario involving pre-excitation of the bending mode of an “O-bound” CO. Also, a portion of the excess kinetic energy is sourced into the potential energy at the time of pre-excitation, and therefore effects arising from reconstruction of the geometry cannot be neglected.

Flipping of the pre-excited “C-bound” CO is observed in almost all trajectories (recall the “C-bound” to “O-bound” isomerization barrier of ∼0.15 eV^24,40). Some trajectories involve single or multiple flipping events, where sometimes the flipped “O-bound” CO configuration is retained, and other times it flips back to the “C-bound” configuration. Roaming is again prevalent in all the flipping trajectories (see Fig. S15(b) in the SI), with CO–surface distances reaching as high as 7 Å in some of them.

The relaxation time of the pre-excited “C-down” CO, obtained from the instantaneous kinetic energy curves (see Fig. S11), spans the range [0.1, 3.2] ps, compared to the “O-bound” CO bending mode relaxation, with effective lifetimes in the range of ∼[0.5, 4] ps. Comparatively, the excitation time constant of the NaCl surface, [0.5, 5.3] ps, as well as that of the neighbouring CO molecules ([0.3, 5.2] ps) is in a similar range. (Note that from fitting the ΔVDOS of ¹CO, an effective lifetime of 3.6 ps is obtained.)

In summary, when exciting CO molecules in the frustrated rotation/bending by 0.26 eV, we find relaxation times on the ps timescale, caused by energy transfer both to the NaCl surface and to neighbouring CO molecules. The relaxation time of “C-bound” molecules is in a similar range compared to that of “O-bound” molecules, and is dominated by energy transfer into the NaCl surface. Further, excitation of both “C-bound” and “O-bound” molecules in the rotational mode(s) leads to roaming dynamics and inversion, i.e., “C-bound” → “O-bound” and “O-bound” → “C-bound”, respectively. The picture emerges that if the bottleneck of energy transfer from an excited CO stretch mode to frustrated rotations has been passed, subsequent relaxation processes (also to the surface) and reactions (e.g., inversion) can proceed fast, on a ps timescale. The situation is similar to other systems, e.g., H/Si(100)⁴¹ and CO/Si(100),⁴² where the slow relaxation of high-frequency modes (the H–Si stretch or CO stretch mode) proceeds indirectly via low-frequency adsorbate modes (e.g., bending), which then couple efficiently to the low-frequency phonon bath of the surface. Again, for CO/NaCl, this is still a preliminary prediction from theory, which needs to be confirmed by experiment. One possible way to do so could be time-resolved VSF spectroscopy.

3.3.2 Time-dependent VSF xxz spectra. To make progress in this direction, we finally report on the time-dependent VSF response of the CO adsorbate(s) for one of the above scenarios, namely vibrational relaxation of a frustrated bending mode in an “O-bound” CO in a one-“O-bound” monolayer (model C). The transient xxz VSF response is computed following the parameterized ssVACF approach, as outlined in Section 2.5.3 and eqn (5), after the “O-down” molecule had been excited with ∼26 bending modes (0.26 eV). Note that only non-parallel and less orientationally disordered molecules are expected to yield a non-zero VSF intensity, and therefore, unlike VDOS curves, capture phase/orientation effects. Being a directly measurable quantity, theoretically computed transient VSF spectra, therefore, offer an advantage over the transient VDOS curves reported above. We emphasize again that our VSF spectra primarily model the CO stretch response, without intermolecular coupling effects. In Fig. 11, we report the transient VSF intensity (|χ_xxz⁽²⁾|²), as well as real and imaginary responses (Re{χ_xxz⁽²⁾}, Im{χ_xxz⁽²⁾}), the latter capturing phase effects, for individual CO molecules after the pre-excitation. In addition, in SI Fig. S16, trajectories are analyzed with the help of distribution functions.

Fig. 11 Molecularly decomposed transient vibrational sum frequency response spectra, obtained from χ_xxz⁽²⁾, of the CO stretch in a one-“O-bound” CO/NaCl(100) monolayer system (model C), where a frustrated rotational mode of the ¹CO (“O-bound”) is pumped with energy equivalent to 26 vibrational quanta (0.26 eV) at Δt = 0 ps. VSF intensity (|χ_xxz⁽²⁾|²), and the imaginary and real parts of χ_xxz⁽²⁾ are plotted in the upper, middle and lower panels, respectively. The spectra are averaged over six trajectories.

From Fig. 11 it is observed that, within a few picoseconds of pre-excitation, a broadened peak in the VSF intensity plot appears for the pre-excited molecule, which can be attributed to the spectrally diffused CO stretch peak. Note that insignificant VSF intensity in the initial few ps results from the orientational and translational disorder (roaming) induced in the respective CO adsorbate (¹CO). After Δt = 12 ps, a faded and narrowed CO stretch peak, with a blue-shifting tail, is found to appear due to dissipation of excitation from this molecule. Consequently, red-shifted trails along the time axis signify the gain of excitation by neighbouring molecules. The loss/gain in intensity therefore depends on the following factors: (i) the overall angular orientation adopted by the molecules, (ii) their perpendicular distance from the surface, (iii) their lateral motion along the surface, and (iv) finite temperature effects. Additionally, all the neighbouring molecules display a negative imaginary response, corresponding to “C-bound” orientations, while the pre-excited “O-bound” molecule shows a positive imaginary response. One can even find some negative low-intensity peaks on either side of the positive peak of ¹CO in the first 8 ps after pre-excitation. Both real and imaginary responses contribute to the resulting lineshape in the intensity plots.

4 Conclusions and outlook

In this work, we studied picosecond-timescale (and, in one example, nanoscond-timescale) relaxation/vibrational energy transfer dynamics exhibited by monolayer/sub-monolayer coverages of CO/NaCl(100) when non-equilibrium vibrationally excited states of the system at particular energies (targeting internal stretch and frustrated rotational modes) are generated. We obtain a somewhat detailed picture of the CO vibrational dynamics on an insulating substrate, disentangling effects from different inter/intra-molecularly coupled degrees of freedom, as well as the moving surface atoms. We do so via quantities like averaged time-dependent internal DoFs, such as CO–surface distances, lateral displacements, angular probabilities, time-dependent vibrational density of states and VSF responses of the CO stretch, as well as the transient difference distance-correlation maps. Specifically, the transient (internal) VDOS curves aid in pinpointing energy flow into various internal modes and track the path of this excess energy pumped into the vibrational mode in the frequency domain. The main findings are as follows:

• There is no significant energy flow to the surface phonons from the CO stretch modes of a “C-down” and “O-down” CO, particularly not within the few-picoseconds timescale, in the presence of only “C-down” neighbours. VET occurs predominantly between the CO stretch modes of neighbouring CO in the monolayer in this case, on a few-ps timescale.

• We see signatures of vibrational energy pooling when exciting CO stretch modes, notably on longer (ns) timescales.

• We predict that exciting the frustrated bending modes to higher levels (with an energy corresponding to a CO stretch quantum) leads to significant, direct energy transfer to the surface phonons, resulting in flipping and desorption events, most often via roaming trajectories. Relaxation times for the corresponding, highly excited frustrated rotations are on the order of a few ps, with the “C-down” CO exhibiting similar relaxation dynamics to those of an “O-down” CO.

• Vibrationally non-adiabatic energy channels are more likely to be prevalent in metastable configurations like the “all-O-bound” monolayer, when energy leaks from the CO stretch mode to low-frequency frustrated librational modes, and also in sub-monolayer coverages, where the flipping threshold is in general known to be lower.

• Metastable geometries found upon pre-excitation of a CO adsorbate in the monolayer to high CO stretch vibrational levels (>10) may be important in the isomerization and pooling dynamics of CO/NaCl(100).

• Time-dependent difference non-linear correlation matrices can offer valuable insights into inter- and intra-mode couplings and the resulting dynamics.

As an outlook, we propose to study in greater detail the dynamics of the CO adsorbates on the NaCl surface on nanosecond timescales and beyond and in larger unit cells/domains, using the NN potential,²⁴ in order to derive further microscopic insights into the isomerization and vibrational pooling dynamics of the CO/NaCl(100) system. Using larger supercells will provide more diversity of oriented configurations and the possibility to investigate VET between non-neighbouring CO(s). While CO ices relevant for astrochemistry are not the focus of the present study, we note that benchmark studies⁴³ have shown that CO–CO intermolecular interactions (which are dominated by dispersion) and anharmonic couplings relevant for V–V energy transfer are sensitive to the choice of exchange–correlation functionals. In particular, range-separated hybrid meta-GGA functionals including explicit dispersion corrections provide a more reliable description than PBE+D2, which in turn neglects higher-order dispersion effects beyond two-body dipole–dipole interactions and may lead to systematically too-weak CO–CO interactions. While the PBE+D2 approach employed here is computationally efficient and benefits from partial error cancellation, allowing semi-quantitative insight and the evaluation of probable trends in V–V transfer mechanisms, more quantitative predictions will require more accurate electronic-structure methods. Therefore, for future studies at higher coverages, further refinement of our NN PES via Δ-machine learning approaches (currently in progress), as well as the potential energy surfaces developed by Guo et al. in ref. 17 and 39, combined with a many-body expansion framework, offers a promising route toward substantially improved quantitative accuracy for vibrational energy pooling and intermolecular energy transfer. Another direction would be to also extend this study to include buried monolayer configurations with distinct overlayers (CO/N₂) as done in experiments,^1,16,18 to investigate isotope-specific VET between overlayer and monolayer molecules. Finally, at low temperatures/energies, reactions such as “O-down” → “C-down” have been suggested to be tunnel-dominated.¹ Quantum effects, which can also cause large isotope effects,¹ thus appear to be important at low temperatures, perhaps also for vibrational relaxation and VET. Some work in these directions is currently underway in our group.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

Supplementary information (SI): technical details on AIMD overflow, velocity swapping approach, VDOS computation, angular distribution probabilities; additional details on vibrational relaxation from CO stretch and frustrated rotational modes in monolayer configurations, time-dependent non-linear correlations, analysis of non-equilibrium trajectories and dynamics with neural network potentials. See DOI: https://doi.org/10.1039/d5cp03781f.

The trajectories run in the non-equilibrium NVE phase starting from different pre-excited configurations can be found in Zenodo under https://doi.org/10.5281/zenodo.17791367. Any other relevant data can be obtained upon reasonable request from the authors.

Acknowledgements

SS and PS thank Giacomo Melani (Pisa) for fruitful discussions on non-equilibrium spectra calculations. This work was funded by the European Union, under ERC Synergy Grant 101163161, “IRASTRO”, “Molecular quantum dynamics in low temperature condensed phase astrochemistry”. Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them. We also acknowledge the North-German Supercomputing Alliance (HLRN) for providing HPC resources. SS also thanks the International Max Planck Research School for Elementary Processes in Physical Chemistry (Berlin) for support.

References

1. A. Choudhury, J. A. DeVine, S. A. Lau, A. Kandratsenka, D. Schwarzer, P. Saalfrank and A. M. Wodtke, Condensed-phase isomerization through tunnelling gateways, Nature, 2022, 612, 691–695.
2. A. Das, Astrochemistry: The study of chemical processes in space, Life Sci. Space Res., 2024, 43, 43–53 , Life Sciences in Space Research Tenth Anniversary Special Issue.
3. A. Ishibashi, G. Molpeceres, H. Hidaka, Y. Oba, T. Lamberts and N. Watanabe, Proposed Importance of HOCO Chemistry: Inefficient Formation of CO2 from CO and OH Reactions on Ice Dust, Astrophys. J., 2024, 976, 162.
4. Y. Oba, T. Tomaru, T. Lamberts, A. Kouchi and N. Watanabe, An infrared measurement of chemical desorption from interstellar ice analogues, Nat. Astron., 2018, 2, 228–232.
5. G. Molpeceres, V. Zaverkin, K. Furuya, Y. Aikawa and J. Kästner, Reaction dynamics on amorphous solid water surfaces using interatomic machine-learned potentials -Microscopic energy partition revealed from the P + H → PH reaction, Astron. Astrophys., 2023, 673, A51.
6. J. C. Santos, K.-J. Chuang, J. G. M. Schrauwen, A. Traspas Muiña, J. Zhang, H. M. Cuppen, B. Redlich, H. Linnartz and S. Ioppolo, Resonant infrared irradiation of CO and CH3OH interstellar ices, Astron. Astrophys., 2023, 672, A112.
7. L. Slumstrup, J. D. Thrower, A. T. Hopkinson, G. Wenzel, R. Jaganathan, J. G. M. Schrauwen, B. Redlich, S. Ioppolo and L. Hornekær, CO desorption from interstellar icy grains induced by IR excitation of superhydrogenated PAHs, arXiv, 2025, preprint, arxiv:2507.07896 DOI:10.48550/arxiv.2507.07896.
8. J. Vogt and B. Vogt, The structure of carbon monoxide adsorbed on the NaCl(100) surface – A combined LEEDand DFT-D/vdW-DF study, J. Chem. Phys., 2014, 141, 214708.
9. H. C. Chang and G. E. Ewing, Infrared fluorescence from a monolayer of CO on NaCl(100), Phys. Rev. Lett., 1990, 65, 2125–2128.
10. J. D. Beckerle, R. R. Cavanagh, M. P. Casassa, E. J. Heilweil and J. C. Stephenson, Subpicosecond transient infrared spectroscopy of adsorbates. Vibrational dynamics of CO/Pt(111), J. Chem. Phys., 1991, 95, 5403–5418.
11. K. Laß, X. Han and E. Hasselbrink, The surprisingly short vibrational lifetime of the internal stretch of CO adsorbed on Si(100), J. Chem. Phys., 2005, 123, 051102.
12. H. H. Richardson and G. E. Ewing, Infrared spectrqscopy of CO on NaCl film and NaCl(100), J. Electron Spectrosc. Relat. Phenom., 1987, 45, 99–111.
13. L. Chen, J. A. Lau, D. Schwarzer, J. Meyer, V. B. Verma and A. M. Wodtke, The Sommerfeld ground-wave limit for a molecule adsorbed at a surface, Science, 2019, 363, 158–161.
14. S. Corcelli and J. Tully, Vibrational Energy Pooling in CO on NaCl(100):Simulation and Isotope Effects, J. Phys. Chem. A, 2002, 106, 10849–10860.
15. S. A. Corcelli and J. C. Tully, Vibrational energy pooling in CO on NaCl(100): Methods, J. Chem. Phys., 2002, 116, 8079–8092.
16. J. A. Lau, L. Chen, A. Choudhury, D. Schwarzer, V. B. Verma and A. M. Wodtke, Transporting and concentrating vibrational energy to promote isomerization, Nature, 2021, 589, 391–395.
17. A. Nandi, P. Zhang, J. Chen, H. Guo and J. M. Bowman, Quasiclassical simulations based on cluster models reveal vibration-facilitated roaming in the isomerization of CO adsorbed on NaCl, Nat. Chem., 2021, 13, 249–254.
18. A. Choudhury, S. Sinha, D. Harlander, J. DeVine, A. Kandratsenka, P. Saalfrank, D. Schwarzer and A. M. Wodtke, Manipulating tunnelling gateways in condensed phase isomerization, Nat. Sci., 2023, 3, e20230006.
19. S. Sinha, A. M. Wodtke and P. Saalfrank, When carbon monoxide goes “upside down”: vibrational signatures of CO at NaCl(100) from ab initio molecular dynamics, Phys. Chem. Chem. Phys., 2025, 27, 7929–7942.
20. H. Ueba, Vibrational relaxation and pump-probe spectroscopies of adsorbates on solid surfaces, Prog. Surf. Sci., 1997, 55, 115–179.
21. H. Arnolds and M. Bonn, Ultrafast surface vibrational dynamics, Surf. Sci. Rep., 2010, 65, 45–66.
22. M. Bonn, C. Hess, S. Funk, J. H. Miners, B. N. J. Persson, M. Wolf and G. Ertl, Femtosecond Surface Vibrational Spectroscopy of CO Adsorbed on Ru(001) during Desorption, Phys. Rev. Lett., 2000, 84, 4653–4656.
23. H. Zhu, B. Chen, V. V. Yakovlev and D. Zhang, Time-resolved vibrational dynamics: Novel opportunities for sensing and imaging, Talanta, 2024, 266, 125046.
24. S. Sinha, B. Mladineo, I. Lončarić and P. Saalfrank, Multidimensional Neural Network Interatomic Potentials for CO on NaCl(100), J. Phys. Chem. C, 2024, 128, 21117–21131.
25. G. Melani, Y. Nagata and P. Saalfrank, Vibrational energy relaxation of interfacial OH on a water-covered -Al2O3(0001) surface: a non-equilibrium ab initio molecular dynamics study, Phys. Chem. Chem. Phys., 2021, 23, 7714–7723.
26. J. Heidberg, E. Kampshoff and M. Suhren, Correlation field, structure, and phase transition in the monolayer CO adsorbed on NaCl(100) as revealed from polarization Fouriertransform infrared spectroscopy, J. Chem. Phys., 1991, 95, 9408–9411.
27. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized Gradient Approximation Made Simple, Phys. Rev. Lett., 1996, 77, 3865–3868.
28. P. E. Blöchl, Projector Augmented-Wave Method, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
29. G. Kresse and J. Hafner, Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal–Amorphous-Semiconductor Transition in Germanium, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
30. G. Kresse and J. Hafner, Ab Initio Molecular Dynamics for Liquid Metals, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558–561.
31. G. Kresse and J. Furthmüller, Efficient Iterative Schemes for Ab Initio Total-Energy Calculations using a Plane-Wave Basis Set, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
32. S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, J. Comput. Chem., 2006, 27, 1787–1799.
33. A. D. Boese and P. Saalfrank, CO Molecules on a NaCl(100) Surface: Structures, Energetics,and Vibrational Davydov Splittings at Various Coverages, J. Phys. Chem. C, 2016, 120, 12637–12653.
34. S. Batzner, A. Musaelian, L. Sun, M. Geiger, J. P. Mailoa, M. Kornbluth, N. Molinari, T. E. Smidt and B. Kozinsky, The NeQuIP package. 2022, An open source package which can be obtained from https://github.com/mir-group/nequip.
35. R. Dettori, M. Ceriotti, J. Hunger, C. Melis, L. Colombo and D. Donadio, Simulating Energy Relaxation in Pump–Probe Vibrational Spectroscopy of Hydrogen-Bonded Liquids, J. Chem. Theory Comput., 2017, 13, 1284–1292 , PMID: 28112932.
36. T. Ohto, K. Usui, T. Hasegawa, M. Bonn and Y. Nagata, J. Chem. Phys., 2015, 143, 124702.
37. S. A. Corcelli and J. C. Tully, Vibrational energy pooling in CO on NaCl(100): Methods, J. Chem. Phys., 2002, 116, 8079–8092.
38. B. C. Ferrari, M. van Hemert, J. Meyer and T. Lamberts, Vibrational Energy Relaxation in Solid Carbon Monoxide, J. Phys. Chem. C, 2024, 128, 21060–21072.
39. J. Chen, J. Li, J. M. Bowman and H. Guo, Energy transfer between vibrationally excited carbon monoxide based on a highly accurate six-dimensional potential energy surface, J. Chem. Phys., 2020, 153, 054310.
40. S. Sinha and P. Saalfrank, “Inverted” CO molecules on NaCl(100): a quantum mechanical study, Phys. Chem. Chem. Phys., 2021, 23, 7860–7874.
41. I. Andrianov and P. Saalfrank, Theoretical study of vibration-phonon coupling of H adsorbed on a Si(100) surface, J. Chem. Phys., 2006, 124, 034710.
42. S. Sakong, P. Kratzer, X. Han, K. Laß, O. Weingart and E. Hasselbrink, Densityfunctional theory study of vibrational relaxation of CO stretching excitation on Si(100), J. Chem. Phys., 2008, 129, 174702.
43. B. C. Ferrari, G. Molpeceres, J. Kästner, Y. Aikawa, M. van Hemert, J. Meyer and T. Lamberts, Floating in Space: How to Treat the Weak Interaction between CO Molecules in Interstellar Ices, ACS Earth Space Chem., 2023, 7, 1423–1432.

---