{'fraction_non_alphanumeric': 0.04155162422878958, 'fraction_numerical': 0.01990660258376794, 'mean_word_length': 4.245485428213839, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Rotation of a single acetylene molecule on Cu(001) by tunneling electrons in STM

We study the elementary processes behind one of the pioneering works on STM controlled reactions of single molecules [Stipe et al., Phys. Rev. Lett. 81, 1263 (1998)]. Using the Keldysh-Green function approach for the vibrational generation rate in combination with DFT calculations to obtain realistic parameters we reproduce the experimental rotation rate of an acetylene molecule on a Cu(100) surface as a function of bias voltage and tunneling current. This combined approach allows us to identify the reaction coordinate mode of the acetylene rotation and its anharmonic coupling with the C-H stretch mode. We show that three different elementary processes, the excitation of C-H stretch, the overtone ladder climbing of the hindered rotational mode, and the combination band excitation together explain the rotation of the acetylene molecule on Cu(100).

Tunneling electrons from a scanning tunneling microscope (STM) form an atomic source of electrons for electronic and vibrational excitations which can be used to manipulate individual atoms and molecules in a controlled manner [1][2][3][4] . The study of acetylene C 2 H(D) 2 rotation on the Cu(001) surface 5 was the first comprehensive and systematic experiment on a single adsorbate manipulation made in combination with STM inelastic electron tunneling spectroscopy (STM-IETS). This method has been established as an indispensable experimental method to gain insight into the vibrationally mediated motions and reactions of single molecules with STM (see, e.g., Ref. 6 and references therein).
The observed in Ref. 5 rotation yield per electron as a function of bias voltage for C 2 H 2 (D 2 ) exhibits a threshold at 358 (266) mV. This corresponds to the excitation of the CH(D) stretch mode which is not a rotational mode. In this respect it is very different from, e.g., the rotation of a single oxygen molecule on Pt(111) surface 4 , where the hindered rotational mode can be directly excited by tunneling electrons 7 . That the onset of a molecular motion/modification is caused by a vibrational mode that is not directly responsible for the observed motion is not unique. Many other examples have been established, e.g., the migration of CO on Pd(110) 8 . The rotation of C 2 H(D) 2 /Cu(001) 5 also demonstrates several peculiar features which have not been previously understood. In this paper we revisit the experimental findings of Ref. 5. Based on DFT calculations and the Keldysh diagram technique we clarify the vibrational modes and elementary physical processes behind the characteristic features of the acetylene 9 rotation on Cu(001).
The experimental data 5 on the C 2 H 2 /Cu(001) rotation is summarized as symbols in Figs. 1 and 2(a). The higher threshold at 358 meV, see Fig. 1, corresponds to excitation of the CH stretch mode. In addition, a lower threshold with a much smaller rotation yield is clearly seen at ∼ 240 meV, we note that this does not correspond to any vibrational energy of the C 2 H 2 /Cu(001), see below. The third feature is the crossover from a single to two-electron process with increasing tunneling current above 10 nA, see Fig. 2(a). This crossover cannot be attributed to a resonant inelastic electron-molecule scattering with rotational excitation as in the case of O 2 /Pt(111) 7 , because of the evident absence of direct excitation of the hindered rotational modes. Moreover, the crossover cannot be attributed to coherent ladder climbing since the reaction order (number of electrons needed) does not approache one at high bias 4,10 .
The identification of the precursor to rotation is not straightforward since the molecule has a number of low frequency vibrational modes. The hindered rotational mode of C 2 H 2 on Cu(001) ∼ 28 meV, whereas the thermally measured barrier for rotation is relatively high ǫ B = 169 ± 3 meV 12 . For this low-energy rotational mode to directly induce rotations, a multi (5-6)-electron processes is expected, whereas in the experiment, as explained above, only a two-electron process is observed at larger tunneling currents.
To analyse the mechanisms of rotation of a single C 2 H 2 /Cu(001) by tunneling electrons we performed DFT cal-  (7), red dashdotted line is the nonlinear process B, Eq.(9), and green dotted line is the combination band process C, Eq.(11). Black solid line corresponds to a sum of all processes. Lines are calculated via Eqs. (3,7,9) with A = 2.5 × 10 −6 , B = 6 × 10 −16 s, C = 1.2 × 10 −8 , and Γ iet (Ω 1 ), the CH mode generation rate Eq.(5), is calculated via the Keldysh technique 3 , the parameters are specified in the text.  . We also estimated the energy barrier ǫ B for C 2 H 2 rotation on Cu(001) to 100 meV using the nudged elastic band method.
The results for the vibrational energies Ω ν and electronhole damping rates γ (ν) eh , ν = 1, . . . 12 are given in Tab. I, in comparison with the previous theoretical 1 and experimen-tal 5,14 results. The last three columns of Tab. I show (in relative units) the components of each modes angular momentum i , where r i , m i are the atomic positions (in the acetylene molecule) and masses, r C.M. the center-of-mass, and δ r (ν) i the atomic displacements in the vibrational mode ν. The high frequency modes (# 1 & 2 in Tab. I) are the symmetric and antisymmetric CH stretch modes. The higher threshold energy of the rotation yield Y(V) = R(V)/I(V) at 358 mV and the corresponding peak in the ∆log(Y)/∆V plot, observed in Ref. 5, indicate that inelastic excitation by tunneling electrons of these modes is a trigger for rotation. Along the reaction pathway there are only two frustrated rotation modes which has angular momentum in the (001) direction, # 6 and # 11.
The full Hamiltonian of the system can be written as a sum of the electronic and vibrational (phonon) parts, where the electronic part H e depends on the normal coordinates of the molecule {δq ν }. Using the Newns-Anderson type Hamiltonian 17 , we write where the indices s(t) and a denote a substrate (tip) and the adsorbate, respectively; the corresponding energy levels are ε s(t) and ε a ({δq ν }). Electronic tunneling matrix elements V t (tipadsorbate) and V s (substrate-adsorbate) give rise to a stationary tunneling current between the tip and the substrate through the adsorbate orbital at applied bias voltage V. The electron occupation functions in the substrate and tip are assumed to be Fermi distributions with the same temperature T but different chemical potentials µ t and µ s , µ s − µ t = eV.
In order to clarify the interaction parts of the full Hamiltonian (1) behind the experimental results on the acetylene rotation, it is convenient to split the rotation rate into three partial processes, where the rates R A (V) and R B (V) are, respectively, the oneand two-electron processes with a higher threshold V ∼ 358 mV, and R C (V) is the one-electron process with a lower threshold V ∼ 240 mV.
To describe the generation of high frequency CH stretch modes, the adsorbate orbital energy in the 1st term of Eq. (2) can be expanded 18 to the first order as where χ is an electron-phonon constant and ε a (0) the unperturbed adsorbate energy, b h is the annihilation operator of the vibrational mode δq h = 2 −1/2 (b † h + b h ) with the frequency Ω h = 358 meV, which is directly excited by the inelastic tunneling current 19 .
The high frequency vibration generation rate then reads 3 , γ (h) eh is the inverse lifetime of the phonon mode h due to electron-hole pair excitation, given in Tab. I. At T = 0 20 Our next task is to find the mechanisms of the energy transfer between the high frequency CH stretch mode and the reaction coordinate (RC) hindered rotational mode. The first possibility is the direct over-barrier rotation excitation due to an inelastic tunneling generation of the stretch mode Ω h 8,21 . The rotation rate corresponding to such single-electron process can be described 22 as a linear function of the CH stretch mode generation rate, A good fit to the experimental data is given by A = 2.5 × 10 −6 , see Fig. 1.. The physical meaning of A is the probability of a CH stretch vibration to excite the C 2 H 2 rotation over the barrier. This value is in reasonable agreement with the well established case of the migration of CO on Pd(110) 8,23 . However, this mechanism cannot describe the two-electron process. In principle, the direct over-barrier rotation excitation from the second excited level of the CH stretch mode provides a twoelectron process. However, this process can be discounted since the short lifetime γ (h) eh −1 ∼ 1 ps gives a nearly two orders of magnitude smaller rate.
Another possibility for the energy transfer is an anharmonic interaction of the CH stretch mode with the reaction coordinate mode. In the simplest case this can be described as a cubic coupling. We expand the vibrational Hamiltonian H ph in Eq. (1) up to cubic terms to account for the coupling between the directly excited (ν = h) and the reaction coordinate (RC) (ν = r) modes, where K h,r,i , K h,r,r are the anharmonic coupling constants, and ν = i is possibly some auxiliary (idler) vibrational mode excited simultaneously with the RC mode. DFT calculation of the anharmonic coupling constants is prohibitively time-consuming. We thererfore construct (see in SM-II&III) a simpler model which takes into account only the pair interactions between the nearest neighbors as springs on rods; the results for the coupling coefficients are shown as colors in Fig 3. The symmetric CH stretch mode #1 couples most efficiently with a pair of the asymmetric rotations #6. In contrast, asymmetric CH stretch mode #2 couples strongly with a pair of non-equivalent phonons, i.e., asymmetric rotation #6 and cartwheel mode #7. The coupling to the symmetric rotation mode #11 is ineffective and the coefficients are approximately two orders of magnitude smaller. Thus, the asymmetric rotation #6 (CH out-of-plane bend) is the most probable candidate for being the acetylene/Cu(001) rotation precursor.
Two different pathways responsible for the two-electron partial process R B are possible. If the barrier height is ε B > 2Ω r , only the second term in Eq. (8) contributes. In this case the rotation starts from the excitation of the high energy symmetric CH stretch mode #1 and two successive decays into pairs of RC mode (# 6) are needed to overcome the barrier. Such a overtone ladder climbing process (see Fig. 4) was first discussed for a desorption of CO molecules in Ref. 4.
If Ω r < ε B < 2Ω r , one pair of the RC stretch excitations is enough to cross the barrier which just gives a correction to the single-electron rate R A . In this case, the only possibility for a two-electron process comes from the third term in Eq. (8). The rotation starts from the excitation of the high energy antisymmetric CH stretch mode (#2), and two successive decays into pairs of RC mode (# 6) and idler mode (#7) are needed to overcome the barrier. This corresponds to the usual singlestep ladder climbing process 20 .
In both cases the generation rate of RC phonons (pairs of phonons) is proportional to the generation rate of the CH stretch phonons, Analysing the Pauli master equations for the one-and two-step RC potential ladder climbing processes shows that the rotation rate is with different formulas for the coefficient B depending on the corresponding anharmonic coupling coefficients K and the ratio of the RC and idler modes linewidth to the detuning The fitted value B = 6 × 10 −16 s reproduce the experimental data above the higher threshold ∼ 360 mV well, see Fig. 1. The hybridizations parameters ∆ s = 200 meV, ∆ t = 12 meV were fixed from the tunneling current I(Ω h ) = 40 nA and the fraction of the inelastic component I inel /I ≈ 1% at V = 400 mV. To simulate higher current levels, the same parameters were used except ∆ t which changes as a function of bias and tunneling current as shown in Fig. 2(b). The calculated rotation rate as a function of tunneling current underlines the importance of the two-electron process at high currents and voltages, see Fig. 2(a). From the dash-dotted lines in Fig. 2(a) which shows the one and two-electron contributions where only the one-electron process, excitation of the C-H mode followed by the anharmonic mode coupling to the RC mode, is responsible for the reaction yield at low current below around 1 nA.
The anharmonic coupling coefficients described above allows for a theoretical extimation of B. If Ω r < ε B < 2Ω r , the estimate is B = 3.5 × 10 −15 s, which is about 7 times larger than the best fit value. However, as noted above, the double excitation of the reaction coordinate mode gives in this case a contribution to R A (V). The resulting value of the coefficient A is then 6 × 10 −3 , i.e., 10 3 times larger than the fitted value 2.5 × 10 −6 . This allows us to exclude the case of a lower rotational barrier from consideration since it would give rise to a much faster one-electron process than seen experimentally. The only possibility left is that ε B > 2Ω r . The theoretical estimate for B is then ∼ 2.4 × 10 −15 s, i.e., approximately 5 times larger than the fitted value. We believe, due to approximations, this value to be in reasonable agreement with the experiment.
As can be seen from Fig. 1, below the higher energy threshold Ω h the rotation rate per electron is very low but non zero above a lower threshold of ∼240 meV. This can be attributed to an inelastic electron tunneling processes which involves a simultaneous combination band 25 generation of two coherent phonons ν = c1, c2. Assuming that the adsorbate energy in Eq. (2) is now a function of these vibrational modes, ε a ({q c1 , q c2 }) , and expanding it in a Taylor series where η = ∂ 2 ε a ({0})/ (∂q c1 ∂q c2 ) and, b c1 and b c2 are annihilation operators of the vibrational modes with frequencies Ω c1 and Ω c2 , and eh damping rates γ (c1) eh , γ (c2) eh . The analysis of this mechanism using the Keldysh-Green's functions shows that the total coherent phonon excitation rate takes the form of a single-phonon excitation rate Γ iet (Ω, V) 3 , Eq. (5), where a single vibrational frequency is simply replaced by the sum of two vibrational frequencies. This gives for the reaction rate, instead of Eq. (7), where Γ iet (Ω c,1 + Ω c,2 , V) = and the coefficient C depends on the parameters of the system as explained in SM-VI. Although it is difficult to make a theoretical estimate of this coefficient, the fitted value to the experimental data is C = 1.2 × 10 −8 . Thus, the combination band single-electron process is about A/C ∼ 200 slower than the process with rotation excitation via CH stretch vibration. This is in reasonable agreement with the fact that the process Eq. (10) occurs in the next order of the perturbation theory compared with the process Eq. (4). The over-barrier rotation mode excitation occurs due to an inelastic tunneling generation of a coherent pair of phonons, e.g. the in-plane bend or wag mode Ω c1 = 141 meV and inplane bend or scissor mode Ω c2 = 118 meV (modes # 4 and 5). Another possibility is that the second of these phonons is the asymmetric out-of-plane rotation Ω c2 = 101 meV (mode #6). The latter, being a hindered rotation mode, may simplify the resulting molecule rotation.
To conclude, we show that the precursor for the acetylene rotation on Cu(001) is the out-of-plane bend (or asymmetric rotation) mode ∼100 meV (# 6 in Table I). Rotation with a higher threshold voltage 358 mV occurs if enough energy stored in the high-frequency CH stretch mode excited by tunneling electrons is transferred to the rotational mode (reaction coordinate). The anharmonic coupling of the CH mode with the hindered rotation mode is found to be responsible for a crossover from a single to multiple electron process for tunneling currents higher than 10 nA. The lower threshold voltage for rotation at 240 mV is attributed to a combination band process of inelastic scattering of tunneling electrons on a pair of lower-energy vibrational excitations of the acetylene molecule.  The details are given for • DFT calculation of the vibrational modes of acetylene molecule on Cu(001) surface and the inelastic tunneling; • Third order anharmonic terms in the vibrational Hamiltonian which are responsible for the energy transfer between the vibrational modes; • Model calculation of the anharmonic coupling constants; • Anharmonic excitation rate of the hindered rotational phonons using Keldysh diagram technique; • Pauli master equation approach to calculate the probability of rotations due to anharmonic decay of high frequency CH stretch mode; • Keldysh-Green's function derivation of the combination band two-phonon rotation excitation.

I. DFT CALCULATION RESULTS
The DFT calculations were performed in SIESTA using a supercell of the Cu(001) with a 4x4 surface and 4 atoms thick slab. Calculational details include a real space cutoff of 200 Ry, Gamma points approximation, double/single-z polarized (DZP/SZP) basis set for the C (DZP), H(DZP), and Cu(SZP) atoms. Computational details for the vibrational frequencies, electron-phonon coupling, and electron-hole pair damping can be found in Ref. 2 . For the nudged elastic band (NEB) calculation more refined computational parameters were necessary, a 5x5 surface, 600 Ry cutoff, 3x3 k-points, and the DZP basis set for all species. Using the refined parameters changed vibrational frequencies by less than 10 %. Although the NEB calculation is fairly converged with computational parameters, the calculated barrier is still approximate because of the basis set superposition error inherent in the SIESTA method. We therefore use the experimental barrier height in the main text.
The equilibrium configuration of C 2 H 2 on Cu(001) and the relaxation of the Cu atoms is shown schematically in Fig. 1. The results for the vibrational energies Ω ν and electron-hole damping rates γ The C-C and CH bond lengths and C-CH bond angle are given in the text.

II. ANHARMONIC COUPLING THE CH STRETCH MODES WITH OTHER VIBRATIONAL MODES OF THE ADSORBED ACETYLENE
In this section we explicitly describe the scheme to calculate the normal coordinates modes of a system (in our case the system involves the acetylene molecule and the surface of copper) and then the anharmonic couplings between them.
Shifting each atom from the equilibrium position { R i } with a displacement vector {δ r i } and expanding the molecule potential energy U({δ r i }) in Taylor series where i, j denote the atomic indices, α, β are the Cartesian indices, U 0 is the potential energy of the molecule in the equilibrium state and coefficients a αβ i j are the matrix elements of the Hessian, Equations of motion for shifted atoms can be written as where m i is a mass of the i-th atom of the adsorbed molecule. For simplicity, we neglect the vibrational modes of the substrate, assuming the copper atoms to be infinitely heavy, and fixed at relaxed equilibrium positions. Substituting δr α i = ξ α i exp (−iωt) into Eqs. (3) one obtain 3N = 12 equations of motion To find the normal modes of the system one needs to solve the secular equation In order to estimate the anharmonic couplings between different modes we expand the potential energy to the next order After rotating to the basis of normal coordinates we obtain where {δ q i } is a set of normal coordinates which represents a solution of the secular equation (5) and K (c) is classical anharmonic coupling between the vibrational modes of the systems. The transformation from the original Cartesian coordinates to the normal coordinates basis has a form where e α m,i is the eigenvector coefficient between normal coordinate δq m and the shift of i-th atom in α direction δr α iδq m = i,α e α m,i δr α i . To quantize the vibrational modes, we introduce a dimensionless displacement vector δq m , so that δq m = δq m (2ω m ) −1 and arrive from coordinates to operators in canonical way δq m → δq m . Then, creation b † m and annihilation b m operators of the corresponding mode can be introduced and the cubic term in Eq. (7) is transformed to where the anharmonic coefficient K m,m ′ ,m ′′ is connected with

III. SPRINGS ON RODS: A SIMPLER MODEL OF ACETYLENE ON CU(100)
In order to estimate the anharmonic couplings K m,m ′ ,m ′′ we introduce the simpler model potential U({ r i }) and imply the scheme discussed in Supplementary material II.
We propose to describe the CH, CC and C-nearest neighbor Cu bonds by springs on rods (the latter to fix the central character of the forces): where ω i j and L i j are the parameters (see in Tab. II) chosen to reproduce the calculated by DFT set of eigenfrequencies. The comparison of the DFT calculated eigenenergies with that within the simpler model are given in table III.  We calculate the anharmonic coefficients K m,m ′ ,m ′′ and find the modes coupled effectively. Figure 3 of the main text shows the dependence of coefficients K 1,n,m and K 2,n,m on mode numbers n and m. These coefficients indicate the most effective ways of decay of the symmetric CH stretch mode #1 and asymmetric CH stretch mode #2. These modes are known to be responsible for the high-energy threshold ∼ 360 meV, they are excited directly by the tunneling electrons.
The symmetric CH stretch mode #1 decays most efficiently via excitation of a pair of equivalent phonons: #4 CH asymmetric in-plane bend or wag, #5 in-plane bend or scissor, #6 out-of-plane bend or asymmetric rotation, #7 cartwheel. Coupling of the symmetric CH stretch mode #1 with the symmetric rotation mode #11 is ineffective since the corresponding anharmonic coefficient is 70 times smaller than that for a coupling with the pair of out-of-plane bend or asymmetric rotation mode #6.
The asymmetric CH stretch mode #2 decays most efficiently via excitation of a pair of non-equivalent phonons, e.g., the pair of the asymmetric rotation #6 and cartwheel mode #7.
This simple estimation of the anharmonic coefficient shows that the rotation of the acetylene molecule is initiated via excitation of the asymmetric rotation mode #6. There are two processes leading to the excitation of the reaction coordinate mode #6, the excitation of a pair of the asymmetric rotation phonons or the excitation of one asymmetric rotation phonon and one phonon of the cartwheel mode #7.

IV. DOUBLE EXCITATION OF HINDERED-ROTATION MODE
In this section we derive explicitly the excitation rate of the reaction coordinate mode. As we have shown in Supplementary material III, the excitation process of the vibrational mode of the acetylene molecule involves two possible pathways, via excitation of the reaction coordinate mode and the auxiliary idler mode or via double excitation of the reaction coordinate mode (terms H ph,1 and H ph,2 respectively, in Eq. (8) of the main text).
We discuss here both of the scenarios using the Keldysh diagram technique. In both cases under consideration the frequencies of the vibrational modes are far from the resonance, Ω h > 2Ω r ∼ Ω r + Ω i and the anharmonic interaction between them can be treated as weak.
In what follows we derive the excitation rate of the RC mode due to the process described by H ph,1 ; the excitation rate due to H ph,2 can be obtained replacing the index i of the idler phonon in all formulas below, with the index r of the RC phonon.
For the description of the effective stationary occupation densities of the RC mode we use the kinetic equation. 3 The anharmonic component of the excitation rate of the RC mode is given by the one-loop polarization operator. Neglecting the temperature corrections, it reads where ρ (ν) ph (ε) is the density of states of the RC and highfrequency modes and n (ν) ph (ε) are the corresponding occupation densities. Formula (11) describes the energy transfer rate to the hindered rotation mode of an adsorbate due to the anharmonic coupling with the CH stretch mode.
We proceed with the calculation of a total excitation rate Γ iet (Ω r , V) of the RC mode After making a substitution of (11) into (12) the total RC excitation rate takes the form (13) The second term in Eq. (13) shows a threshold dependence on bias voltage, because it is proportional to the high-frequency mode occupation numbers n (h) . It can be shown that the first term in (13) can be omitted due to the fact that n (h) ph (Ω a + Ω r ) ≪ n (h) ph (Ω h ). The total excitation rate of the RC phonons due to the anharmonic term H ph,1 is then and due to H ph,2 is Note that both rates are proportional to the CH stretch mode excitation rate Γ iet (Ω h ) and to small phonon densities of the RC and idler vibrational modes far from the resonance.

V. LADDER CLIMBING WITH TWO-AND ONE-STEP PROCESSES
In this section we calculate the excitation rate of the RC phonons R B (V) using Pauli master equations. As discussed in Supplement IV, depending on the height of the rotational barrier, two possible processes can lead to initiation of the rotations of the acetylene molecule on Cu(001). In this section these processes will be considered separately as they involve ladder climbing of two different types.
In case when the rotational barrier height is Ω < ε B < 2Ω, the one-step ladder climbing process takes place 4 . The excitation rate of this process Γ iet,1 (Ω r , V) is given by (14) and the relaxation rate is γ (r) eh , thus the Pauli master equation can be written as The stationary solutions (in respect to P 0 ) for m = 1 states in the localization potential of the RC mode can be written as is defined as a probability rate to overcome the localization potential barrier and in our notations it is the excitation rate from the first excited level If the rotational barrier height is 2Ω < ε B < 3Ω then the excitation rate is assumed to be dominated by the two-phonon anharmonic coupling with the high-frequency mode. The pair of RC phonons excitation rate Γ iet,2 (Ω r , V) is given by (15)  − mγ (r) eh + (m + 2)(m + 1)Γ iet,2 (Ω r , V) P m . (18) The stationary solutions (in respect to P 0 ) for m = 0, 1, 2 states in the localization potential of the RC mode can be written as Reaction rate R B (V) in this case is a sum of the excitation rates from the first excited state 6Γ iet,2 (Ω r , V)P 1 and from the second excited state 12Γ iet,2 (Ω r , V)P 2 , The reaction rate in both cases is a quadratic function of the RC phonon excitation rate which is a feature of the two-step ladder climbing process and differs only in a proportionality coefficient. For a more accurate estimation of proportionality coefficient between the second power inelastic tunneling current and rotation rate we need to take into account the fact that the inelastic tunneling current is a sum of several components which arise from the scattering of tunneling elections on all relevant vibrational modes. In our case we are interested in two CH stretch modes, symmetric #1 and asymmetric #2. According to our DFT calculation the inelastic transmission through these modes are T (1) = 1.3 × 10 11 (s×V) −1 and T (2) = 5.9 × 10 11 (s×V) −1 . In what follows we will take this into account. introducing the probability factors to excite the symmetric CH stretch mode ζ = T (1) T (1) + T (2) −1 = 0.18 and the asymmetric one (1 − ζ) = 0.82.
Substituting the expression for the excitation rate of the RC phonons (14) and (15) into Eqs. (17) and (20) correspondingly and using the expression for the high-frequency phonons occupation densities we obtain the proportionality coefficient B (1) [Eq. (9) of the main text] between the excitation reaction rate and the phonon generation rate R (1) B (V) = B (1) Γ iet,1 (Ω h , V).
Using K h,r,i ≈ 32 meV, γ (a) = 0.2 ps −1 , γ (h) = 0.7 ps −1 , γ r = 0.7 ps −1 , ∆ = Ω h − Ω r − Ω i ≈ 208 meV (Ω i ≈ Ω r ) we obtain B (1) = 3.5 × 10 −15 s. It is 7 times larger than the coefficient B = 6 × 10 −16 s obtained from the best fit of the experimental data. Moreover, in case of a lower reaction barrier Ω < ε B < 2Ω the process of double excitation of the reaction coordinate phonons gives a contribution to the linear part R A (V) of the rotation probability. The rotation rate can be estimated then as R A (V) ≈ Γ iet,2 . Then the impact of this process into the proportionality coefficient A can be written as A (2) = 4πK 2 h,r,r ζτ (h) ρ (r) ph (Ω h − Ω r ).
Using the same parameters as above for the coefficient B (1) , we arrive to the estimated value of A (2) = 6 × 10 −3 , which is three orders of magnitude larger than the value of A obtained from the fitting of the experimental data. Evidently, this makes the case of a lower reaction coordinate barrier Ω < ε B < 2Ω hardly possible . Analogously, the coefficient B (2) between R (2) B (V) and Using K h,r,r ≈ 32 meV, γ (h) ph = γ (1) = 1 ps −1 , γ (r) = 0.7 ps −1 , ∆ = Ω h − 2Ω r ≈ 208 meV we obtain that B (2) ≈ 2.4 × 10 −15 s, 5 times larger than the best fit value B = 6 × 10 −16 s.
Thus the estimated value of the coefficient B is slightly larger than the value obtained from the fit of the experimental data. We believe that our simpler model overestimates the anharmonic coupling coefficients K h,r,r ; for a better agreement with the experimental data we have to take ≈ 1.5 times smaller values K h,r,r ∼ 20 meV.