Spin-Forbidden Carbon–Carbon Bond Formation in Vibrationally Excited α-CO

Fourier transform infrared spectroscopy of laser-irradiated cryogenic crystals shows that vibrational excitation of CO leads to the production of equal amounts of CO2 and C3O2. The reaction mechanism is explored using electronic structure calculations, demonstrating that the lowest-energy pathway involves a spin-forbidden reaction of (CO)2 yielding C(3P) + CO2. C(3P) then undergoes barrierless recombination with two other CO molecules forming C3O2. Calculated intersystem crossing rates support the spin-forbidden mechanism, showing subpicosecond spin-flipping time scales for a (CO)2 geometry that is energetically consistent with states accessed through vibrational energy pooling. This spin-flip occurs with an estimated ∼4% efficiency; on the singlet surface, (CO)2 reconverts back to CO monomers, releasing heat which induces CO desorption. The discovery that vibrational excitation of condensed-phase CO leads to spin-forbidden C–C bond formation may be important to the development of accurate models of interstellar chemistry.

. Polarization-dependent FTIR spectra. Absorption features assigned to (left) CO and (right) C 3 O 2 /CO 2 are shown. Red and black traces correspond to s-and p-polarized spectra, respectively. Data correspond to the ca. 300 layer sample represented in Fig. 1 but do not include the background subtraction illustrated by Fig. S2. Figure S2. Illustration of background subtraction of FTIR spectra. The black trace is an uncorrected p-polarized FTIR spectrum of a ~300 layer 12 C 16 O sample acquired after the sample had been exposed to 60000 laser shots. Orange lines show the portions of the spectrum that were chosen to form the background interpolation, resulting in the subtracted spectrum shown in red. Figure S3. Initial and final p-polarized FTIR spectra for three samples of α-CO subjected to the excite-probe procedure. The initial CO coverage and isotopic purity of the dosing gas are indicated in the top-right corner of each panel. The weak features at ca. 2000 cm -1 in (a) are due to iron carbonyl impurities that were not fully removed from the dosing gas by the pentane ice bath. Table S1. Energies (in Hartrees) calculated for the singlet and triplet pathways. In addition to the species involved in reactions [1][2][3], results are provided for the MECP and IM geometries identified in closer examinations of the reaction mechanism. The CCSD(T) total energies include the ZPE correction obtained from the ωB97M-V treatment. Additionally, the relative total energies of the two spin states of each species are provided, as well as the T 1 diagnostic values, ; the values are all < 0.02, indicating that the single-reference approach is sufficient here. 1

S1. Spectral Assignments
The presence of CO 2 as an impurity in the dosing gas used to prepare the samples considered here is evidenced by the appearance of the CO 2 antisymmetric stretch fundamental transition in the pre-excitation spectra in Figure S3. The spectra reported here show a broad feature centered at ca. 2341 cm -1 with a sharper feature at 2346 cm -1 . The sharper feature lies between two frequencies observed for a pure monolayer of CO 2 on a NaCl(100) surface (see Table S5); 2 the broader feature lies in a range that is in fair agreement with previous reports of IR spectra of CO 2 embedded in CO matrices. [3][4][5][6][7][8] Figure S4. Lineshape of the shot-dependent feature at ca. 2250 cm -1 . The difference absorption spectrum for a ~ 300 ML sample exposed to 55000 laser shots is shown in black. Fitting this to a two-peak Lorentzian lineshape [CO] 0 results in the red traces, where the component peaks are also shown. Dashed lines indicate previously reported absorption frequencies assigned to C 3 O 2 and C 3 O in CO matrices. 7 Figure S4, the shot-dependent feature at ca. 2250 cm -1 shows a lineshape consistent with two underlying transitions at ~2242 and 2248 cm -1 . The infrared absorption spectrum of C 3 O 2 diluted in a CO matrix shows a strong absorption at ~2242 cm -1 , as well as a weaker combination band at ~2400 cm -1 ; 4 these features have been clearly identified in CO ices following irradiation with protons or UV photons. [3][4][5] Though fairly weak, there is some evidence for the ca. 2400 cm -1 combination band in our S6 spectra (Fig. S5, blue area). We can thus confidently assign at least the 2242 cm -1 component as arising from formation of C 3 O 2 in our sample. Figure S5. Demonstration of the absence of spectral evidence for formation of C 3 O and C 2 O. The difference absorption spectrum for a ~300 ML sample is shown in gray, with a smoothed version shown in black for clarity. The red trace shows the spectrum simulated by assuming that the side-peak at 2248 cm -1 represents the formation of C 3 O, leading to the appearance of a weak feature at 1913 cm -1 that is not seen experimentally. The blue trace is a spectrum simulated with the assumption that the area under both ~2250 cm -1 peak components contributes to the calculation of [C 3 O 2 ]. Additionally, the C 2 O region of the spectrum is highlighted in green.

C 3 O 2 , C 3 O. As shown in
In previous reports of C 3 O 2 formation in CO ices, the lineshape of the ca. 2250 cm -1 band observed here always shows some degree of asymmetry, with some reports achieving sufficient resolution to reveal a sidepeak at ~2248 cm -1 consistent with the weaker component of Fig. S4. In some cases, the evolution of this side-band was noted to correlate with the increase of a lower-intensity transition below 1950 cm -1 assigned to C 3 O. 3,7 However, as noted by Trottier et al., 5 it is also possible that this side-peak corresponds to a matrix perturbation of the C 3 O 2 band. Figure S5 shows an experimental post-excitation difference spectrum in the regions where C 3 O 2 and C 3 O are known to absorb. If the 2248 cm -1 side-band corresponds to C 3 O, the absorption coefficients reported by Jamieson and coworkers 7 may be used to simulate the intensity of the ca. 1913 cm -1 transition of C 3 O. As shown in the red trace of Figure S5 (pink area), this would likely be below our detection threshold. However, we may also simulate the expected intensity for the C 3 O 2 combination band at ca. 2400 cm -1 ; if we consider the area under the 2248 cm -1 feature to contribute to the value used in calculating [C 3 O 2 ] (Section S2a), we find better agreement for this weak transition (blue curve, blue area in Fig. S5).
Additionally, we note that the formation of C 3 O would occur through the reaction For this to take place in our samples, the reactive products of reactions [1-2] must be generated in close proximity, or diffuse together before reacting with a CO molecule. Previous studies of CO/NaCl(100) have shown that the highest-energy adsorbates produced through VEP are quite dilute, 13 so the likelihood of two nearby CO molecules achieving sufficient excitation to proceed through [1] is quite low. Because of this low density of reaction centers, as well as the limited diffusion inherent in a 7K ice, we do not consider C 3 O as a significant contributor to our observed spectra, and we take the dominant carrier of the ca. 2250 cm -1 signal in our spectra to be carbon suboxide.

C 2 O .
In high-energy irradiation of CO ices, the formation of C 2 O is typically monitored through its absorption at ca. 1980-1990 cm -1 . 7 This region where C 2 O would absorb is highlighted in Figure S5 (green area), and shows that this molecule is not accumulated in sufficient quantities to observe.

S2. Quantitative Analysis of FTIR Spectra a. Calculation of Column Densities
The integrated absorbances of the spectral features summarized in Table S5 are  where is the gas-phase integrated absorption coefficient (in cm molecule -1 , see Table S5), is the = 1 number of surfaces probed, and is the angle of incidence of the FTIR light source. The 2D density = 45°a ssociated with a single layer within the α-CO crystal is given by 6.27 10 14 molecules cm - where is the lattice constant for α-CO at 8K. 14 0 For CO, the integrated absorbance is taken to be the sum of the LO and TO peak areas, obtained by numerically integrating the background-subtracted spectra in the 2135-2145 cm -1 region. The CO 2 column density was obtained in a similar manner by taking the numerical integral of the 2330-2355 cm -1 spectral region, so that the results include contributions from both components of the observed band. The absolute S8 CO 2 column densities obtained in this manner reflect not only that formed due to interaction with the laser, but also the CO 2 initially present in the samples, as well as the roughly linear increase arising from background adsorption. The reported column densities in Figure 2 have been corrected for these effects (see Fig. S6). Figure S6. Growth of CO 2 column density over the time scale of the excite-probe procedure. The excite-probe program was applied to a ~150 ML sample of CO with (red) and without (black) laser excitation, and column densities were calculated as described in Section S2a to assess the increase in CO 2 that occurs solely due to adsorption of residual gas in the UHV chamber. This linear increase is indicated by the blue line, which was subtracted from all CO 2 column densities to correct for this background signal.
For C 3 O 2 , the overlap of the spectral feature with the CO phonon band complicated evaluation of the integrated absorbance; results were calculated from the ~2240-2280 cm -1 range of the difference spectra (where the spectrum of the initial sample is subtracted from the post-excitation result). This range was selected to minimize contributions from the CO phonon side-band, but may not fully account for the integrated absorbance arising from carbon suboxide.

b. Yield Analysis of Concentration Data
As can be seen from Figure 2, the extent of product formation is far too small to fully explain the observed CO depletion. Defining the yield of CO depletion as where , we find that ~33, 53, and 58% of the initial CO coverage was lost after 60000 laser shots for samples prepared with ~ 150, 300, and 500 layers, respectively. To determine the portion of this CO loss that [CO] 0 can be attributed to chemical transformations, we consider CO 2 and C 3 O 2 to be the only products and CO the only reactant, so that the net reaction is given by With this, and noting that in all three cases the calculated amount of CO 2 formed is greater than or equal to C 3 O 2 , the maximum amount of CO consumed by chemical reactions can be estimated as .

× Δ[CO 2 ]
Dividing by shows that for all samples, less than 0.2% of the total CO lost can be accounted for by production of CO 2 and C 3 O 2 . These percentages, as well as the absolute changes in concentrations and percent product yields (given by ), are provided in Table S6.
Assessing the dependence on thickness is complicated by the fact that the thickest sample was prepared with natural abundance CO, whereas the two thinner samples were prepared using isotopically pure 12 C 16 O; the enhanced presence of isotopic impurities likely impacts the pooling dynamics in the ca. 500 ML sample. 15 However, some trends may be identified, in particular when the column densities are converted to more intuitive units. The number of layers in the sample (given by the sum of the column densities of all three species) may be multiplied by the thickness of a layer of α-CO (given by the ratio with = 2 / 3 3 2.22 10 22 molecules cm -3 ) to yield the total thickness of the sample (Fig. S7a). 14 Column densities × (molecules cm -2 , obtained by multiplying the data in Fig. 2 by ) may then be converted to concentrations 2 (molecules cm -3 ) by dividing by this thickness; the volumetric concentrations are used to obtain the fractional composition of our sample over the course of the excite-probe procedure (Fig. S7b).
Considering these results, we note that while increasingly thick samples show larger percent depletion of CO (Table S6), the overall change in sample composition decreases. In all three cases the product formation saturates when the photoproducts constitute ca. 0.02% of the total sample, though this slightly decreases as sample thickness increases.

S3. Comparison of Different Levels of Theory
The geometries of the reactants, TS, and products for reactions [1][2][3] were also obtained at the B3LYP/aug-cc-pVTZ and PTPSS-D 3 /6-311+G* levels of theory. These results are compared to the ωB97M-V/6-311+G* results in Figure S8a, showing similar energetic pictures for reactions [1-3] as is given in Figure 3 of the main text.
As noted in the main text, the calculated barrier heights for reaction [1] in our work show that the SP(S 1 ) geometry is lower in energy than the CO 2 + C products, in contrast to the results of Barreto and coworkers. 16 To investigate this discrepancy, repeated calculations at the MP2/aug-cc-pVTZ level of theory were performed to obtain the geometries and frequencies of all species involved in reaction [1], followed by single-point CCSD(T)/aug-cc-pVQZ calculations; this was the approach used in the prior work. The results of these theoretical methods, as well as the B3LYP results mentioned above, are shown in Figure S8b. Only the MP2 level of theory gives a singlet SP energy greater than that of the products to [1]; the CCSD(T) energies of the MP2 geometries give a similar energetic picture to that obtained from DFT for the singlet pathway, although the triplet pathway shows that the conversion to CO 2 + C is effectively barrierless.
Further, an IRC calculation on the MP2/aug-cc-pVTZ singlet SP shows the same behavior represented in Figure 4B, indicating that for all model chemistries used here, the SP(S 1 ) geometry is in actuality a transition state to the formation of IM(S 1 ) from 2CO molecules. showing that IM(S 1 ) 2CO decomposition is a barrierless process. → Figure S9. FSM calculations regarding IM(S 1 ). The FSM results for (blue) 2CO IM(S 1 ) and (black) IM(S 1 ) → → CO 2 + C( 1 D) are shown. Energies do not include ZPE effects and were obtained at the ωB97M-V/6-311+G* level of theory. The "reactant" and "product" geometries used as inputs for the FSM calculations are also shown.
One could also imagine that the IM(S 1 ) geometry leads to formation of cyclic C 2 O, c-C 2 O, which is identified as a stable isomer on the singlet surface. This could occur through dissociation, IM(S 1 ) c-C 2 O → + O, and the remaining O atom would then react with a CO molecule to form CO 2 . Alternatively, the cyclic structure could be formed by direct reaction of IM(S 1 ) with CO, IM(S 1 ) + CO c-C 2 O + CO 2 . These → S14 possibilities were explored using our ωB97M-V/6-311+G* framework, yielding the energy diagram of Figure S10. In both cases formation of the cyclic structure requires surmounting a reaction barrier; given the results shown in Figure S9 it seems far more likely that IM(S 1 ) dissociates to form 2CO molecules, rather than proceeding to products through a cyclic intermediate.

b. Molecular Dynamics Simulations
While the calculations summarized above indicate that the lowest-energy singlet pathway does not provide a clear connection between IM(S 1 ) and products, it should be noted that IM(S 1 ) formed from SP(S 1 ) will at first possess considerable energy. This could provide access to alternative, higher-energy pathways capable of forming the CO 2 + C( 1 D) products. To determine this possibility, Born-Oppenheimer molecular dynamics (BOMD) simulations were performed in Q-Chem using the ωB97M-V/6-311+G* framework.
These trajectories propagated the IM(S 1 ) geometry with an initial kinetic energy (KE) given by one of the two values illustrated in Figure S11. From this diagram it can be seen that the lower-KE choice (KE 1 ) represents formation of IM(S 1 ) at the ZPE of SP(S 1 ), which is nearly degenerate with the potential energy of the CO 2 + C( 1 D) products. We also consider a higher-KE choice (KE 2 ), which would correspond to formation of vibrationally-excited SP(S 1 ) formed from CO molecules at the uppermost energy levels known to be populated through VEP. For both KE choices, fifty sets of initial velocities were generated which randomly distributed this kinetic energy along the six normal modes of IM(S 1 ) illustrated in Fig. S11. S15 Dynamics were simulated for 500 time steps of duration 20 a.u. (0.484 fs), using the microcanonical (NVE) ensemble so that the total energy remains constant.
In all cases, the structure dissociated into 2 CO fragments moving away from each other. This can be seen by the representative atom-atom displacement plots in Fig. S12-S13, where at long times, the only two atom-atom separations within bonding tolerances are the two separated CO monomers (O 1 -C 2 / O 3 -C 4 ). All other displacements diverge, representing increasing separation between the CO centers-of-mass (COMs).
Defining as the first time point where the separation between COMs of the CO fragments is > 3 Å, we can construct the histograms in Figure S14, showing that even when given sufficient energy to exceed the barrier to form CO 2 + C( 1 D), trajectories initialized on the singlet pathway rapidly lead back to 2CO molecules.
The constant-energy conditions used to propagate the trajectories represented in Figures S12-S13 differ from the experimental situation, where the molecular system is embedded in a 7K environment that can serve as a heat sink. Thus, by assuming energy conservation, the BOMD simulations performed here should provide an enhanced probability of observing conversion to CO 2 + C( 1 D), particularly in the E 2 simulations where IM(S 1 ) is initialized with sufficient energy to energetically access this geometry. Thus, the results here clearly show the unlikelihood of forming CO 2 on the singlet pathway within our theoretical framework.

S5. Curve-Crossing in Reaction [1]
The reaction pathway calculations represented by Fig. 3 suggest a crossing of the singlet and triplet surfaces of reaction [1] in the vicinity of the transition state. To study this in more detail, the structures generated in the triplet IRC (Fig. 4A) were used to simulate geometries sampled at even intervals along the reaction coordinate of [1]. Single-point CCSD(T)/def2-qZVP calculations on these geometries in the singlet and triplet states yield the energy diagram of Figure S15a, which represents the singlet and triplet diabatic surfaces through the conversion of 2CO to CO 2 + C. The crossing point for these non-optimized surfaces is found to occur shortly before the triplet transition state geometry; this was used as an initial structure to determine the minimum-energy crossing point (MECP) along the reaction pathway of [1].
The MECP optimization was performed at the ωB97M-V/6-311+G* level of theory, where the branching-plane method was used to identify the crossing point for the (reference) singlet state with the first SF-TDDFT triplet state. 18 The CCSD(T)/def2-qZVP energy of this species, including ZPE correction, is shown in Figure S15b. As shown in Fig. S16, the MECP geometry largely resembles SP(S 1 ), but has a torsional angle that is more consistent with that of a CO dimer in an α-CO crystal (Table S7).   19 The "dimer" structure formed by considering the two highlighted molecules is shown in the bottom-right, alongside the ωB97M-V/6-311+G* geometries of SP(S 1 ) and the MECP for [1] (middle-right). All dimer structures may be described in terms of the six geometrical parameters shown in the top-right; these parameters are presented in Table S7, showing that the three angular coordinates are strikingly similar between the α-CO and MECP cases. Table S7. Geometrical parameters for the CO dimer structures considered here. Parameters as defined in Figure  S16 are shown for the singlet/triplet SP, the ωB97M-V/6-311+G* MECP, as well as a CO dimer in α-CO. The parameters for the crystal structure are obtained using the same parameters specified in Fig. S16 -i.e. an equilibrium CO bond distance of 1.128 Å and a lattice constant of 5.646 Å -where the average bond length is given by obtained from the Morse oscillator wavefunctions associated with the ground singlet state of CO. 20

S6. Intersystem Crossing Rates a. Fermi's Golden Rule Model
We consider non-radiative transitions of the type where i and f are two electronic states of molecule A with differing spins. The N-dimensional vectors and contain the N quantum numbers specifying the vibrational levels of the initial and final states. We will restrict our considerations to exoergic transitions, so that as defined in Eq. (S1) is always positive. Δ We consider that the energy released by Eq. (S1) is absorbed by the phonons of the surrounding CO crystal, as discussed in Section S6b.
Within the Born-Oppenheimer approximation, the initial and final states may be expressed as Here, replaces the typically-employed delta function, due to the ability of the CO phonon bath to (Δ ) account for the energy discrepancy in Eq. (S1) (see Section S6b).
Given the observation that the singlet and triplet surfaces in Fig. 3 cross, it seems plausible that the SOCs of interest in the current system will depend on the nuclear configuration. We therefore will also consider a spin-vibronic (SV) mechanism for ISC. 22 In this framework, the vibronic states given in Eq. (S2) are coupled through the spin-vibronic Hamiltonian, , expressed as an expansion of the purelyelectronic around a reference geometry, , where is the th normal coordinate of . With this, the SV-ISC rate is where the SV matrix elements are given by From the state-to-state rates given by Eqs. (S3-S4), the total FC-and SV-ISC rates for a particular initial state i are then calculated by summing over the manifold of final states f,

b. Treatment of the Phonon Density of States
As noted above, we consider the surrounding CO matrix and its vibrational density of states (DOS) to serve as a means of ensuring energy conservation. To do this, rather than using a delta function in our rate expressions, we use the phonon DOS of the surrounding crystal at the energy difference between (Δ ), the initial and final vibronic levels of . Given the low temperature of our system, we restrict our considerations to exothermic processes where as defined in Eq. (S1) is strictly positive. Δ Previous studies of the phonon DOS for α-CO indicate that the phonon spectrum is highly structured and cuts off above 150 cm -1 . 23 For simplicity, we use a three-dimensional Debye DOS in Eqs. (S3-S4), with a cutoff frequency of cm -1 . = 200 Calculations of the rates discussed below were also tested using the NaCl(100) vibrational DOS reported previously 13 as well as an approximation to the α-CO phonon DOS; 23  where and are the reduced mass and frequency (in cm -1 ) associated with mode in state .
We note that evaluation of the nuclear integrals involved in the SV-and FC-ISC rates requires one to express the vibrational wavefunctions given above in terms of the same set of coordinates. However, the this prevents the expression of the nuclear integrals as products of one-dimensional results.
To avoid this difficulty, we instead invoke a parallel approximation. We assume that the normal modes of the singlet and triplet states are similar enough that the same coordinates can be used to express vibronic wavefunctions in both states, so that the normal coordinates of the singlet state are given by . = + Δ S22 Figure S17 shows the mode-matching convention used; as vibrational modes with imaginary frequencies are not appropriate for treatment with a harmonic oscillator model, we neglect in our rate evaluations. 1 We then calculate the coordinate shifts , which are obtained from the equilibrium geometries and Δ as well as the (normalized) atomic displacement vector associated with mode , , by With these approximations we obtain our final expressions for our vibrational wavefunctions,  Table S1. Mode-specific parameters were obtained by our computational work and are summarized in Table S8.
We then express the SV matrix element which couples the level of ⟨Ψ │ │Ψ ⟩ = ( 2 , 3 , 4 , 5 , 6 ) SP(S 1 ) to the level of SP(T 0 ) as = ( 2 , 3 , 4 , 5 , 6 ) where the T superscripts have been removed from the normal modes for clarity. The EOM-CCSD/6-311G formalism was used to evaluate the electronic contributions, where the SOC matrix element ⟨ │ │ ⟩ and derivative couplings were taken to be those obtained at the initial (singlet) state geometry. The derivative terms were obtained from the SOC curves given in Figure S18; values were obtained by taking the derivatives at the appropriate values of given in Table S8. The SOC calculations at the equilibrium Δ geometry gives = 33.86 cm -1 . We note that similar results are obtained when the triplet ⟨ │ │ ⟩ /ℎ geometry is chosen as the reference point.

S23
In practice, calculating these rates requires restricting our consideration to a limited range of vibrational levels, which we do here by only considering SP(S 1 , ) and SP(T 0 , ) states with .    respectively) as well as the corresponding time constants (given by the inverse of the corresponding rate) for transitions where we restrict our calculations to include all final states SP(T 0 , ) with . The total SP(S 1 ,0)⇝SP(T 0 ) ∑ ≤ number of final states is also given, as is the energy (in cm -1 ) of the highest level relative to SP(T 0 , ). Further 0 increasing was not found to change the resultant rate (for up to = 15). total states max.

b. Other CO Dimer Structures
While the above treatment shows that rapid ISC is possible between vibrational levels of the saddle points identified for reaction [1], there are several other points on the singlet surface which may be considered to provide access to the triplet pathway. In particular, the MECP identified in Section S5 seems an intuitive choice. However, a frequency analysis at the ωB97M-V/6-311+G* level of theory yields two imaginary frequencies for this geometry in the singlet state (Table S3); thus, application of the ISC frameworks considered here would require further limiting the dimensionality of the nuclear wavefunctions.
This would most likely lead to an over-estimation of the Franck-Condon factors describing the nuclear contributions to the overall rates. Given the strong resemblance between SP(S 1 ) and the singlet MECP, we consider the rates calculated above to be a reasonable estimate of ISC transitions out of dimer structures resembling the singlet MECP, and note that the DFT treatment of the gas-phase model should serve only as a qualitative guide to the processes occurring in the crystal samples studied here.
The identification of the cyclic IM(S 1 ) provides an additional possible means by which the triplet surface may be reached, via relaxation to the linear IM(T 0 ) OCCO geometry. However, as noted by Jamieson and coworkers, 7 this structure has not been observed experimentally, and our work indicates that if IM(S 1 ) is formed it should rapidly reconvert to the 2CO reactants (Fig. S16, Section S4b). Additionally, the large difference in geometries for the singlet/triplet structures indicates that the Franck-Condon overlap would be poor, leading to low ISC rates. Thus, we do not consider IM(S 1 ) as a likely point of access for the productforming triplet pathway.