Condensed-phase isomerization through tunnelling gateways

Quantum mechanical tunnelling describes transmission of matter waves through a barrier with height larger than the energy of the wave1. Tunnelling becomes important when the de Broglie wavelength of the particle exceeds the barrier thickness; because wavelength increases with decreasing mass, lighter particles tunnel more efficiently than heavier ones. However, there exist examples in condensed-phase chemistry where increasing mass leads to increased tunnelling rates2. In contrast to the textbook approach, which considers transitions between continuum states, condensed-phase reactions involve transitions between bound states of reactants and products. Here this conceptual distinction is highlighted by experimental measurements of isotopologue-specific tunnelling rates for CO rotational isomerization at an NaCl surface3,4, showing nonmonotonic mass dependence. A quantum rate theory of isomerization is developed wherein transitions between sub-barrier reactant and product states occur through interaction with the environment. Tunnelling is fastest for specific pairs of states (gateways), the quantum mechanical details of which lead to enhanced cross-barrier coupling; the energies of these gateways arise nonsystematically, giving an erratic mass dependence. Gateways also accelerate ground-state isomerization, acting as leaky holes through the reaction barrier. This simple model provides a way to account for tunnelling in condensed-phase chemistry, and indicates that heavy-atom tunnelling may be more important than typically assumed.


Experimental Methods
shows a schematic diagram of the experimental apparatus used in the current work. The molecular beam dosing procedure used for sample preparation has been described previously. 1 Briefly, a fixed number of molecular beam pulses of CO gas ( 12 C 16 O: Sigma-Aldrich; 99.99% 12 C and 99.95% 16 O, 99.9% chemical purity, 13 (100) crystal surface held at 25K to form the isotopically pure monolayer. After monolayer dosing, samples were cooled to 22K and ~100 layers of 12 C 16 O were added. Cryogenic temperatures were achieved by coupling the sample to a dual-stage helium cryo-cooler (RDK-408D2, Sumitomo), and the sample was held inside a liquid nitrogen-cooled cold shield during measurements to limit warming due to room-temperature surfaces 10 in the laboratory. The sample temperature was monitored and controlled using a Model 335 Lakeshore cryogenic temperature control unit. Sample characterization was performed using a Fourier Transform Infrared (FTIR) spectrometer (VERTEX 70V, Bruker) operating in external mode using either an InSb or an MCT detector cooled by liquid nitrogen.
Following sample preparation and characterization, the sample was cooled to 19K and a pulsed infrared 15 laser fixed at 2138.6 cm -1 was used to induce the v = 0 to v = 1 vibrational excitation in the 12 C 16 O overlayer 2 and the vibrational energy was efficiently transported to the the 13 C 16 O or 13 C 18 O monolayer. 3 For the 12 C 16 O monolayer buried beneath a 12 C 16 O overlayer, vibrational energy transport is not possible; here, the monolayer was pumped directly with IR laser light at 2152.5 cm -1 . A voltage pulse-activated shutter controlled the number of laser shots used to pump the sample. The laser was then blocked and the 20 temperature was raised to the value of interest. The time dependent absorption spectra were then observed using a Bruker FTIR spectrometer equipped with the Rapid-Scan functionality, the light source of which was aligned to be collinear with the excitation laser. From these measurements, the sample absorption spectrum was obtained as a function of time. Changes in the spectrum were determined by using the absorption spectrum of the unexcited sample as the reference channel. These experiments were performed 25 for three different monolayer isotopologues as well as various surface temperatures. The infrared light used to excite the sample was generated by a difference frequency mixing setup wherein the output of a tuneable dye laser pumped by the second harmonic of a seeded 10 Hz Nd:YAG laser mixes with the 1064nm Nd:YAG fundamental light in a LiIO3 crystal.
To ensure that there was no influence of exposure to the FTIR light source on the resultant decay of Odown population, a set of two experiments were carried out where FTIR spectra were acquired at ~90 5 second intervals following excitation of the sample. In one set of experiments, it was ensured that the glow bar was only illuminating the sample during spectral acquisition (~30 seconds per spectrum). Thus, the effective exposure to the FTIR light source was reduced to ¼ of that characteristic for the other measurement, where the FTIR light source was always irradiating the sample. As shown in Fig. S2 the two kinetic traces are indistinguishable, demonstrating that exposure to the glow bar does not influence the 10 measured lifetimes of O-down molecules.

Data Analysis
It is obvious from Fig. S4 that the integral over the absolute value of the C-down feature (area B + C) is considerably larger than area C -B and hence produces kinetic traces with larger signal-to-noise ratios. In the following section we will show that area B + C (method 3) is in fact proportional to the O-down 15 concentration, and analyzing the time dependence from the kinetic traces using B + C gives more reliable time constants than methods 1 and 2, in particular when the rate of flipping is very high (e.g. at > 23K) or the laser excitation efficiency is low (e.g., 12 C 16 O monolayer samples).
In this section we describe in detail how FTIR absorbance spectra are used to obtain the time dependent concentrations of O-down CO on NaCl(100). We first perform a baseline correction (section S2a). We then 20 consider how to extract from the spectral data the peak integrals that are related to the O-down population at each time step (section S2b), and relate the obtained peak areas to O-down population using an exciton model (section S2c).

Baseline correction
The baseline of the obtained absorbance spectra shifts over the course of a lifetime measurement which can last many hours, hence a baseline correction must be implemented. To achieve this, a straight line was fitted through the baseline of spectra obtained at different times in the measured decay and subtracted. Figure S3 shows examples of two experimentally obtained spectra without (black) and with (orange) this 5 baseline correction.

Relative O-Down Concentrations
The relative fraction of O-down molecules in the monolayer can be obtained from the difference absorbance spectrum (see Fig. S4) by determining the area A under the O-down peak (method 1). It can also be calculated from the decrease in the C-down fraction relative to the unexcited sample by integrating 10 the whole C-down feature (method 2). Since its shape arises from a drop in absorption intensity as well as a broadening and red-shift of the band, the integral of the difference spectrum consists of a positive (area B) and a negative (area C) contribution and is given by the area C -B. These two methods give essentially the same result (A ≈ C -B), as the decreased C-down fraction is equal to the increased O-down fraction.

15
It is necessary to show that the absorbance of CO molecules in the sample is linearly proportional to their population. To do this, we simulate absorbance spectra for a simplified matrix of CO molecules with the help of a well-known exciton model. 4,5 The monolayer is taken to have a 1 × 1 unit cell, where CO monomers are either aligned vertically upward (C-down) or downward (O-down). The position of O-down molecules were chosen randomly. From the simulated spectra (Fig. S5) we determine the correlation 20 between the fraction of C-down molecules and the integrated absorbances calculated by the three methods discussed above. We show that for O-down fractions generated in our experiments, all three methods give integrals which are linear in the O-down population (Fig. S6). We find that the integrated absorbances obtained by method 3 discussed above give better signal-to-noise in the kinetic traces (Fig. S7).
The monolayer is considered to be a collection of N molecules, where the ground and first-excited 25 vibrational state of isolated molecule i are given by and ′ , respectively, where the relative energy of the first excited state is given by ℎ̃0. The zeroth-order collective N-molecule states are then perturbed by their pairwise dipole-dipole interactions.
In the ground state, the collective zeroth-order wavefunction is given by The collective single-excitation wavefunctions, where a single adsorbate is in its first excited state, form an 5 N-fold degenerate sub-space with energy ℎ̃0 (in the case that all molecules have the same orientation), so that the l th single excitation exciton wavefunction is given by where l = 1, 2, 3, …, N. The i th component in (2) then corresponds to the zero th -order single excitation state 10 where molecule i is in its first vibrational excited state indicated by a prime.
To determine the true (coupled) excitonic wavefunctions and eigenvalues, we must diagonalize the perturbation Hamiltonian in the basis consisting of the individual terms included in (1) and (2); i.e., 0 and the N single-excitation , where = is the Kronecker delta. Neglecting electrostatic interactions beyond electric dipole-dipole, the perturbation Hamiltonian is written as Here, and are the magnitude and unit vector of the dipole moment of molecule , respectively, is the unit vector from molecule to molecule , and is the intermolecular distance. For a polar diatomic molecule, the dipole moment is parallel to the bond axis of the molecule, so also indicates the orientation 20 of the molecule .
The diagonal element of the perturbation matrix associated with the 0 th vibrational ground state (1) is given by (1) = 〈 0 | (1) | 0 〉 = 00 00 ∑ ∑ < =1 , (5) and those associated with each excited state k represented in (2) may be written as where 01 is = 0 → 1 transition dipole moment (0.105 Debye). 4 We now apply N-fold degenerate perturbation theory to calculate the first-order energy corrections (1) , where (1) = 〈 | (1) | 〉 ( , = 1, ⋯ ) are the matrix elements from Eqs. (6) through (8), and is the Kronecker delta. Diagonalization of (1) gives N eigenvalues, (1) = Δν ( ) ℎ , as well as the associated eigenvectors Ψ , whose total energy is given by = ℎ (ν 0 + Δν ( ) ) = ℎ̃( ) . We note that this applies The integrated absorption cross-section for excitation of the monolayer to the l th single-excitation exciton state obtained from the diagonalization discussed above is given by (10) 20 The absorption spectrum of chosen monolayer geometry (̃) was then calculated as a Gaussian convolution of the discrete stick spectrum [∆̃( ) , ̅ ( ) ] associated with the N coupled single-excitation exciton states. Simulations were performed using a 30×30 grid of CO molecules; increasing the grid size was not found to substantially change the resultant spectra. Resulting difference absorbance spectra, is the spectrum for the 100% C-down monolayer, are shown for 5 various O-down coverages in Fig. S5.
From Fig. S5, we can determine the relationship between integrated absorbances and the percentage of molecules flipped in the monolayer, by calculating for grids with varying fractional O-down populations. The results of applying the three methods to the simulated spectra are shown in Fig. S6A.
Lifetime measurements were limited to samples where not more than ≈30% of the molecules were flipped. An example of a typical lifetime measurement experiment is shown in Fig. S6B, where we can 15 clearly see that the O-down peak area ( ) is ∼0.13 fraction of the full monolayer ( 0 ), indicating that ~30% of the molecules were flipped. The integrated absorbances calculated from method 3 were then used to determine the concentration of O-down molecules for the reasons discussed above.

Table of Lifetimes
Table S1 provides the flipping lifetimes obtained for all samples studied in this work. The lifetimes are obtained by fitting the O-down population data, calculated as described in the previous section, to a single exponential decay and then averaging the resultant fit parameters over several experiments. The uncertainties are estimated with 95% confidence interval. Rates were fit to an Arrhenius model to extract 5 the parameters given in Table S2.

Transition State Theory
The rate constant for the conversion of O-bound to C-bound CO as a function of temperature can be obtained using transition state theory (TST), where 0 is the classical barrier to reaction from a previously reported two-dimensional DFT calculation 5 for a coverage of ¼ monolayer. 6 (See also Sec. S6 below.) The partition functions of the transition state and the O-down reactants, ‡ and , are obtained by treating all motion (internal stretch, frustrated rotation, and frustrated translation) as simple harmonic oscillators, ‡ = ∏ 1 where frequencies ‡ and were obtained from Ref. 6

WKB Theory
A commonly employed model used to consider tunneling through a barrier is the semi-classical Wentzel-Kramers-Brillouin (WKB) approximation. The tunneling rate in the WKB in mass-weighted coordinates is defined by the following equation Here, is the attempt frequency taken here to be the frequency of the frustrated rotation, ( ) is the potential energy along a reaction coordinate s; the integration limits 1 and 2 are the turning points at energy E on each side of the barrier. Note that we have used the RHS of Eq. 13 where the Z coordinate is parametrized w.r.t the θ coordinate (see Fig. S8 (a)). Here, = C + O and = C O / are the CO total and reduced mass, respectively, and eq ≈ 1.14Å is the CO bond equilibrium distance. Using the 2D 10 PES from Ref. 6  We also used the WKB model to simulate thermally activated tunneling from 5 to 31 K on the same 20 PES used above. This relies on a thermal average of the state specific WKB tunneling rate constants ( ).
25 All vibrational eigenstates m dominantly localized in the "O-bound" well, up to the classical isomerization barrier, were included in the WKB thermal averaging. In this way we include the influence of vibration on the WKB tunneling rate. The results are shown in Fig. S9. The isotope-dependent thermal WKB rate constants were fit to an Arrhenius form-those fitted parameters are reported in Table S2.
Quantum Rate Theory: One-phonon Assisted Thermal Rates 5 Fermi's Golden Rule (FGR) combined with system-bath approaches is a widely used formulation for the study of vibrational relaxation dynamics in condensed-phase systems, where FGR gives a transition rate between a chosen initial and final state of the system. This rate depends on the coupling between the two states as well as the density of final states. Here, we use such a system-bath approach to develop an FGR model of isomerization rates for CO at the NaCl (100) surface; here, the "system" is defined by the "O-10 bound" / "C-bound" initial/final states of the CO molecule on the NaCl (100) surface, and the "bath" corresponds to the phonons of the CO/NaCl interface (i.e. motion of the NaCl surface atoms plus the neighboring CO molecules). System and bath modes are coupled by a vibration-phonon coupling operator which enters the FGR rate expressions as specified below. The basic physics is illustrated in Fig. S10.
Below, we describe our system, the chosen initial and final states, and the coupling model in more

Model and FGR Rate Expressions
Total Hamiltonian In our system-bath approach, our system (an inverting CO molecule on a NaCl (100) surface) is treated in two dimensions (2D), considering the important vibrational adsorbate modes involved in the isomerization reaction, namely and , which are coordinates for the frustrated rotation and perpendicular 25 CO-surface motion, respectively. The system modes are coupled to a harmonic bath comprising the NaCl (100) surface phonons as well as the low frequency vibrational modes of the neighboring CO adsorbates up to the Debye frequency of the surface, i.e., 222 cm -1 . The model Hamiltonian of the total system is then expressed as where ̂=̂2 is the 2D adsorbate Hamiltonian, ̂ describes the phonon bath, and ̂ the coupling between adsorbate (system) and bath.

System Hamiltonian and Eigenstates
The 2D system Hamiltonian is given by 10 where and = 2 are the total CO mass and moment of inertia, respectively, and fixed coordinates are indicated in the potential, . The definitions of these fixed coordinates can be found in Ref. 6 .
As mentioned earlier, this 2D Hamiltonian describes the perpendicular CO-surface motion and frustrated CO rotation for an unburied layer of CO on a NaCl(100) surface. For the potential, ( , ), we adopt a function which was calculated from gradient-corrected density functional theory (DFT) with 15 dispersion corrections for a coverage of ¼ ML as obtained in Ref. 6 . Details about the anharmonic potential are described in that reference. In particular, one finds a separation of "O-bound" and "C-bound" configurations of about 613 cm −1 , with a classical barrier for "O-bound" to "C-bound" reaction of 574 cm −1 .
Note that for the experimental coverage (one ML, buried under an overlayer), slightly different potential parameters are expected, an effect which we neglect in the present theoretical model. 20 In Ref. 6 , corresponding 2D system eigenfunctions ( ) were calculated by diagonalizing the system Hamiltonian, where = ( = ) will be used to represent "O-bound" ("C-bound") localized eigenstates. Up to energies very close to the barrier, all eigenstates are either localized in the "O-bound" or "C-bound" well. 25 These 2D system eigenstates correspond to hindered rotation and CO-surface stretching vibration levels, which when excited, will either relax (downward transitions) or become further excited (upward transitions) via system-bath coupling. Our model accounts for only transitions occurring from "O-bound" to "C-bound" vibrational states, i.e., transitions between states and in different wells of the doublewell potential, while transitions within the same well are neglected. All eigenstates corresponding to either 5 "O-bound" or "C-bound" configurations below the classical barrier for the isomerization reaction have been considered for the rate calculations (see table S3 below). These eigenstates have been computed for all three CO isotopologues studied in experiment as well as for 12 C 18 O. 89 Since the transition rates in our model are sensitive to the system eigenfunctions, we found it necessary to increase the grid size for the system DOFs compared to Ref. 6 when obtaining the 2D system eigenfunctions in order to obtain converged rates. To extend our calculations to higher grid sizes than possible with ordinary direct-product grids, we utilized a Potential Optimized Fourier Grid Hamiltonian Here min[ 2 ( , ); ] represents the value of which minimizes the 2D potential for a given value of . Further, is the -value of the minimum of the 2D potential in the "C-bound" configuration ( ≈ 3.33 Å). The 1D zero-order states are then given by: The 1D eigenvalue problems can be solved with the FGH method using large, spatially uniform grids. 9 To ensure orthonormality of the eigenfunctions, a Gram-Schmidt reorthogonalization routine is also performed on the final 2D eigenfunctions. Further, the number of zero-order product states are chosen such that the total energy does not exceed certain energy cutoff values. Specifically, we selected only those 1D 10 zero-order states of the coordinate whose energies were below the desorption limit (~ 1400 cm -1 ) and the total energy ( + ) was limited to 4000 cm -1 . Details concerning convergence with respect to grid sizes will be discussed below.

Bath Hamiltonian and System-Bath Coupling
Next We assume a bilinear system-bath interaction, where coupling is linear in system modes , and in bath modes, , which allows only for single-phonon excitations / relaxations in the bath. Then, the systembath coupling Hamiltonian can be written as Here, ( , ) are system-bath coupling functions, taken to be linear in both system modes, whereby coefficients / represent the coupling constants of the corresponding degrees of freedom (DoFs) to the bath modes and is the equilibrium CO bond length (1.141 Å). 0 and 0 are reference coordinates 5 which, again, need not be specified because they have no influence on the matrix elements entering the rate expressions below.
The coupling constants / could, in principle, be calculated from multi-dimensional system-bath potentials and Taylor expansion to first order 10 . Instead, we make two simple approximations. First, we assume that both system modes couple equally to the phonon modes, i.e., we set = = . Following 10 a well-known general procedure, 11 can be obtained from the so-called bath spectral density, ( ), in a discretized form Here, = = / are now equally spaced ( ) frequencies of the phonon bath up to a cutoff frequency , taken to be the Debye frequency of the NaCl(100) surface, 222 cm -1 . The second simplifying 15 approximation assumes the bath spectral density is proportional to the vibrational density of states (VDOS), where is a proportionality constant.
The parameter serves to scale the thermal rates to be approximately within the range of experimental 20 values. It is worth noting that since appears as a multiplicative factor in the FGR rate expression given below, it has no influence on the kinetic isotope effect.
The VDOS is evaluated for a monolayer of CO molecules on a NaCl (100) surface, with one "O-bound" CO (see Fig. S11a) from ab initio molecular dynamics (AIMD) calculations at 30 K, using periodic DFT with the Perdew-Becke-Ernzerhof (PBE) functional and Grimme's D2 dispersion correction Specifically, a 2√2×2√2×3 slab model was used with four CO molecules per unit cell, and only the uppermost NaCl layer and the CO molecules were allowed to move in the AIMD calculations. The VDOS is then extracted 5 from the velocity-velocity autocorrelation function following an approach similar to ref. 12

10
We utilize the FGR formulation for transition rates where, initial and final states, Φ and Φ respectively, are products of the anharmonic system states The "upward" ("downward") rates lead to the system being excited (relaxed) from an "O-bound" state with a lower (higher) energy, , to a "C-bound" state with higher (lower) energy, , following absorption (emission) of a single phonon of frequency , cf. In the following rate expressions, we use a continuum limit, such that the FGR rates take the form The advantage here is that no broadening factor is needed, which is usually required when the δ distribution is taken as a Lorentzian in the FGR rate expression. Moreover, the Lorentzian broadening procedure appears numerically less robust and more sensitive at the low temperatures for tunneling rates and is therefore not 20 further considered here. We finally note that the continuum expressions (34)  The zero-temperature ground state tunneling rates obtained in this way are contrasted in Table 1 of the main text with corresponding semiclassical WKB rates. In Table S2, the trends in Arrhenius parameters derived both from FGR and WKB (and also TST) are compared. It must be noted that a one-to-one comparison 20 cannot be made because of the different approximations underlying the WKB and FGR models and the different dimensionality of the tunneling path, and also as we will see below, due to the temperature range considered for obtaining the Arrhenius parameters in both cases. Qualitative differences are striking, however.
In practice, absorbed the partition function giving a modified ' = / . (The weak T-dependence of was neglected.) We set ' = 1 × 10 −12 atomic units throughout. This choice of ' gives reasonable thermally activated rates when compared to experiment. For ground state tunneling, the rate is 2.6 × 10 −5 s -1 for the heaviest isotopologue, 13 C 18 O. The experimental upper limit for this isotope for the buried (overlayer) system is 3 × 10 -8 s -1 (see section S7). Since our model system is free of CO overlayers, the 5 tunneling rate is expected to be faster compared to the overlayer system.

Limitations of the FGR Model
We would like to emphasize that the FGR model developed for this work is highly simplified. For example, unlike the experimentally-studied buried monolayer samples; the theoretical FGR model has no overlayers. The presence of overlayers would restrict motion of the CO adsorbate in the direction, which 10 would of course strongly affect the overall isomerization rate. This is why we do not aim to compare the absolute FGR and experimental rates for isomerization, but rather, demonstrate that the large non-intuitive isotope effects observed experimentally could be explained by a state-specific tunneling picture.
Other important approximations represent directions for future work and are summarized as follows: (i) The spectral density treats every phonon mode as being essentially equally capable to couple 15 to both system modes. Of course, in reality different phonons will couple differently to the individual system modes, which could strongly impact the overall tunneling rates.
(ii) No two-phonon and higher contributions to rates are included, excluding therefore also elastic scattering events which could arise from absorption and emission of the same phonon, connecting initial and final system states of the same energy, i.e., connect pairs of degenerate 20 system states. However, note that the DOS used in our bath model is not zero-valued near zero frequency owing to the frustrated motions of neighboring CO molecules. Hence nearly degenerate states can couple and give rise to non-negligible rates in our model.
(iii) A bilinear coupling model, Eqs. (24) and (25), was used, which is certainly only approximate along the coordinate due to the asymmetric double-well potential along this coordinate. 25 (iv) The system states are considered only in 2D.
Despite these and further simplifications, our minimal model already accounts for the qualitative effects seen in experiment. Specifically, it predicts: an Arrhenius dependence of tunneling rates with effective activation energies lower than the barrier height for a bare monolayer; the occurrence of (few) "tunneling gateways" which dominate the process; and large isotope effects with non-intuitive mass dependence.
Specific results (beyond those in the main text) are described in the next subsection. As noted earlier, the 5 results discussed below are restricted to the isotopologues studied experimentally, though the same treatment was also applied to 12 C 18 O.

Arrhenius Parameters
We performed FGR total rate calculations for all three isotopologues from 5-31 K. See Fig. 2c of the 10 main text. Two regimes are seen: a deep tunneling regime (T < 20 K), where rates are only weakly dependent on T, and a thermally activated regime exhibiting Arrhenius behavior. As in experiment, the mass-ordering of rates is counter-intuitive both in the deep tunneling and in the activated tunneling regimes.
The FGR model shows an activated tunneling regime that is shifted slightly to higher temperatures compared to experiment. For this reason, we compare Arrhenius values taken from a shifted temperature 15 range of 26.8-31 K (see table S2 and inset in Fig. 2c). This chosen range is 4.2 K wide, comparable to experiment. The shift of our activated tunneling regime and therefore also the onset of deep tunneling regime to higher temperatures relative to experiment can be understood by the fact that our theoretical model is for a non-buried adsorbate layer (at a coverage of 1/4). The experiment, on the other hand, probes a buried adsorbate with a higher expected barrier (due to restricted motion) and therefore reduced tunneling.

20
Of course, quantitative deviations between experiment and our simple model emerge. Most notably, the mass-ordering of activation energies differs; whereas for 12 C 16 O is 28 and 52 cm -1 below those of 13 C 16 O and 13 C 18 O, respectively, the experimentally-derived result for the lightest isotopologue lies between the two heavier ones. However, as demonstrated by Figure 2b, the correlation between Arrhenius parameters is captured by the FGR results. Finally, while the FGR model captures qualitatively the fact that 25 heavier particles can tunnel faster than lighter ones, the precise ordering of isotopologues seen in experiment is not reproduced. This is not surprising given the simplicity of the model.

Importance of Tunneling Gateways
State-to-state ( ) and initial state-selected ( ) isomerization rate constants were calculated for all three experimentally studied isotopologues according to Eqs. (34), (35) and (36). The dominance of 5 tunneling gateways for the heaviest isotopologue, 13 C 18 O, at T=27K is demonstrated in Fig. 3a of the main text, by plotting the state-to-state breakdown of total rates. In Fig. S12, we present similar state-to-state breakdowns of total rates for all three experimentally studied isotopologues at T = 27 and 30 K. These plots  States localized around = ±180° (upper panels) correspond to "O-bound" and states around = 0° 20 (lower panels) to "C-bound" configurations. We emphasize that only states with a high degree of excitation along the direction of hindered rotation participate in gateways.
Gateways are also important for deep tunneling at low T. These are evident in Fig. S16 as sticks near Em=700 cm −1 , in particular tor 12 C 16 O and 13 C 16 O. These gateways act as leaky holes through the barrier, leading to a dramatically weaker T-dependence than seen for WKB tunneling and a non-intuitive mass 25 dependence.

Convergence of System Wavefunctions and Ground State Tunneling Rate
The calculation of initial and final-state wavefunctions ( , ) and ( , ) and corresponding tunneling matrix elements turned out to be a numerical challenge, because the latter depend sensitively on the grid chosen to calculate and represent the wavefunction. In this work, we systematically increased the grid size of our basis for evaluation of the 2D eigenfunctions in order to obtain converged tunneling matrix 5 elements, , for the computation of transition rates. The different grid sizes used are shown in Fig. S17.
Figures S18 through S22 show results of our convergence study for system state energies and FGR matrix elements for selected tunneling gateways. Notice that the energies of the states must be converged to far better than 1 cm −1 to obtain converged matrix elements. Without this level of convergence problems arise. For example, Fig. S18

Attempted Measurement of the Ground State Tunneling Lifetime for 13 C 18 O
To determine the rate of ground-state flipping from O-down to C-down, the O-down sample was prepared as described in section S1 at a surface temperature of 7K and monitored for about 120 hours (Fig.   S23). The change in integrated absorbance within this time is only ∼3%. Based on this, we calculate a lifetime on the order of 1 year, which should be considered a lower limit.     , that is localized in the "O-bound" well of the potential can reach a higher-lying state in the "C-bound" well by either relaxing / absorbing a surface phonon with frequency (+), leading to an upward rate ↑ (see the corresponding FGR expression later), or a lower-lying state by exciting / emitting a surface 5 phonon with (possibly another) frequency (−), leading to a downward rate ↓ . Boltzmann averaging over all possible initial states and summing over all possible final states gives a T-dependent, total isomerization rate which is strongly isotope dependent.  . State-to-state breakdown of total FGR rates for three isotopologues at T = 27 K (left) and T = 30 K (right). Each vertical bar indicates the fraction of the total rate represented by the state-to-state upward rate and downward rate (color coded), describing transitions from "O-bound" states at energy Em to "Cbound" states at En. All energies are with respect to the classical minimum of the "C-bound" configuration. Figure S13. Ratio of initial state-selected (upward plus downward) rates km to total rates at T =27, 30 K for ¼ ML (a) 12 C 16 O, (b) 13 C 16 O, (c) 13 C 18 O on NaCl(100) vs. initial "O-bound" states | ⟩. The numbers on the x-axis give the energies of the initial states | ⟩ relative to the classical minimum of the "C-bound" configuration. The barrier height is 1187 cm -1 on this scale and the lowest energy "O-bound" 5 states are at 685, 684 and 682 cm -1 for the three isotopologues. Note: States involved in giving dominant contributions to rates have been highlighted in the figure with m → n. When more than one final state contributes to the initial-state selected rate km, it is indicated by m → n1, n2, … and so on.  In each panel, the angle of CO with respect to the surface normal (θ) is the vertical axis, where θ=0° represents the C-O bond parallel to the surface normal in the "C-bound configuration". The horizontal axis is the distance of the CO molecule from the surface (Z). 5 State energies are given relative to the classical "C-bound" minimum in the PES. State numbers are w.r.t to the "C-bound" ground state being State 0.  = 1400, = 1800 (orange box)) was used for calculating the 2D eigenfunctions used in the rate computation for all three isotopologues discussed above.    The points show results of measurements and the line is a best fit exponential corresponding to 2 year lifetime. As we cannot rule out that this is noise, the measurements are only useful to put a lower limit on the lifetime of the ground state O-down isomer. 5