The Lowest-Energy Isomer of C2Si2H4 Is a Bridged Ring: Reinterpretation of the Spectroscopic Data Based on DFT and Coupled-Cluster Calculations

The lowest-energy isomer of C2Si2H4 is determined by high-accuracy ab initio calculations to be the bridged four-membered ring 1,2-didehydro-1,3-disilabicyclo[1.1.0]butane (1), contrary to prior theoretical and experimental studies favoring the three-member ring silylsilacyclopropenylidene (2). These and eight other low-lying minima on the potential energy surface are characterized and ordered by energy using the CCSD(T) method with complete basis set extrapolation, and the resulting benchmark-quality set of relative isomer energies is used to evaluate the performance of several comparatively inexpensive approaches based on many-body perturbation theory and density functional theory (DFT). Double-hybrid DFT methods are found to provide an exceptional balance of accuracy and efficiency for energy-ordering isomers. Free energy profiles are developed to reason the relatively large abundance of isomer 2 observed in previous measurements. Infrared spectra and photolysis reaction mechanisms are modeled for isomers 1 and 2, providing additional insight about previously reported spectra and photoisomerization channels.


Introduction
The C 2 Si 2 H 4 molecule has been the subject of very few literature studies, despite being a bicyclobutane/butadiene analog of fundamental interest. Holme et al.
[1] predicted the existence of over a dozen cyclic C 2 Si 2 H 4 stationary points at the Hartree-Fock (HF) level, assigning the lowest-energy isomer as the three-membered ring silylsilacyclopropenylidene. Inclusion of correlation energy using third-order many-body perturbation theory (MBPT(3)) left the isomer ordering unchanged, though it is worth noting that a modest 6-31G(d) basis set was employed. Holme remarked that a three-membered-ring isomer seemed an "unlikely choice" due to the enhanced strain destabilization as compared with competing four-membered-ring structures. Structural compression is energetically favorable when electron-correlation stabilization dominates over strain destabilization. If this is the case for C 2 Si 2 H 4 , one might expect that the maximally compressed isomer, a bridged bicyclobutane analog, should be the global minimum. One such isomer was considered by Holme et al., but it was predicted to be over 10 kcal/mol higher in energy than silylsilacyclopropenylidene.
In the first experimental study on C 2 Si 2 H 4 , Maier et al. reported infrared (IR) measurements supplemented by density functional theory (DFT) vibrational simulations [2]. In their workup a mixture of silane structures was thermalized by pulsed flash pyrolysis at 1500 K, and the resulting mixture was rapidly frozen out in an Ar matrix at 10 K before IR characterization. Signatures of silylsilacyclopropenylidene were detected and upon further irradiation the corresponding ring-opening photolysis products were observed (i.e., (silylethynyl)silylene and ethynylsilylsilylene). Accompanying calculations were performed using the BLYP generalized-gradient approximation (GGA) exchange correlation functional (XCF) and assigned silylsilacyclopropenylidene as the global minimum, in agreement with the prediction of Holme et al. Due to this apparent agreement between theory and experiment, this assignment of the lowest-energy isomer of C 2 Si 2 H 4 has stood for over 20 years.
The purpose of the current study is to apply more sophisticated computational methods to investigate the delicate balance between electron-correlation stabilization and strain destabilization in C 2 Si 2 H 4 . Ground-state DFT methods, and their extension to excited states via time-dependent (TD) DFT [3][4][5], are the most widely-used approaches for describing electron correlation in atoms, molecules, and solids [6]. While superior to the HF self-consistent field (SCF) theory, the popular BLYP and B3LYP approaches have major shortcomings, which often result in qualitatively incorrect predictions [7]. In light of recent advancements, the default choice of XCF needs to be augmented, updated, or supplanted entirely [8]. Since it is unclear which of the hundreds of existing functionals perform best overall, several popular functionals were chosen from the BLYP [9,10], PBE [11,12], and Minnesota functional families [13][14][15][16] for evaluation in the context of the current application.
The quality of XCFs can be ranked as being on one of five rungs of Jacob's ladder, with the fifth rung approaching the "heaven" of chemical accuracy [17]. Here the term "chemical accuracy" takes the colloquial meaning that a method is predictive to within 1 kcal/mol. Functionals including only a local density approximation to the correlation functional [18] represent the first rung, while those incorporating GGA correlation represent the second rung. After density gradients, Laplacians can be included in the mix, resulting in meta-GGA functionals, or the third rung of Jacob's ladder [19]. Hybrid functionals are on the fourth rung of Jacob's ladder, often referred to as hyper-GGAs, while the fifth rung of Jacob's ladder is represented here by double-hybrid (DH) DFT methods [20,21], which semiempirically blend DFT and wave-function theory. Only functionals from the highest three rungs will be considered here.
Electron correlation can be treated within a more hierarchical framework by using wave-function methods, which, unlike DFT, are systematically improvable to the exact solution of the Schrödinger equation. In the ground state, approximations based on the coupled-cluster (CC) hierarchy [22][23][24][25][26][27] are among the most rapidly convergent approaches, with the CCSD(T) method often referred to as the "gold standard" of quantum chemistry. While expensive, CCSD(T) energies extrapolated to the complete basis set (CBS) limit are widely considered benchmark-quality for molecules near their equilibrium geometries. In situations where degeneracy is commonly encountered (e.g., when crossing transition states, when accessing excited states, or when breaking or forming bonds), CCSD(T) can fail [28], and in these cases one can instead turn to new generations of CC methods, such as the left-eigenstate completely renormalized (CR) CC method [29][30][31][32] called CR-CC(2,3). The CR-CC(2,3) approach provides an accurate description of single-bond-breaking processes [33,34], biradical and magnetic molecules [35,36], and transition states [37,38]. Excited states can be accessed in a straightforward manner using the equation of motion (EOM) CC formalism [39][40][41][42], resulting in the CR-EOMCC(2,3) approximation and its size-intensive variant δ-CR-EOMCC (2,3) [31,[43][44][45][46]. Experience shows that reliable energetics are produced for ring-opening [47] and bond-rearrangement processes [48] when CR-CC(2,3) is applied on top of appropriate complete-active-space (CAS) SCF geometries. The steep O(n 7 ) scalings of CCSD(T) and CR-CC methods prohibit their use for all but the smallest systems, but fortunately they are readily applicable to the title molecule.
The motivation for this work spans several fields. Mapping interstellar silicon reaction networks is one application [49], as C 2 Si 2 H 4 may result from bimolecular collisions between the highly-abundant acetylene and its silicon-based analog disilyne. This bimolecular reaction and subsequent isomerization can also be considered as a model system for understanding preferred bonding configurations in coordinate complexes [50], bulk silicon surface-adsorption [51], and defect-inclusion processes [52].
Many of the present modeling decisions were guided by prior benchmarking work on the structure and spectroscopy of hydrogen-free silicon carbide clusters [53,54]. This study continues along these lines, investigating the performance of various modern DFT methods by comparing against spectroscopic values and high-level computational results. The structure of the paper is as follows. The methods employed are detailed in Section 2, the results and discussion are covered in Section 3, while concluding remarks and final recommendations for density functional usage are offered in Section 4.

Energy-Ordering C 2 Si 2 H 4 Isomers
The full potential energy surface of C 2 Si 2 H 4 includes many energetically-similar local minima involving strained three-and four-membered heterocyclic rings. Dozens of unique cyclic species can be imagined, with some resembling cyclobutane or methylcyclopropane, but, unlike their hydrocarbon analogs, C 2 Si 2 H 4 isomers differing by a transannular hydrogen migration are often separated by only a few kcal/mol. Due to the diversity of possible bonding configurations, ranking such structures by energy demands a high-level treatment of the correlation energy and a large basis set. Figure 1 collects structures and dipole moments for 10 low-lying C 2 Si 2 H 4 isomers characterized as PES minima using B3LYP/Def2-QZVP. Dipole moments are also provided, as derived from B2PLYP/Def2-QZVP electronic densities. Saddle-point structures of any order were discarded. Isomer 2 was previously characterized as the global minimum and isomer 6 is reported here for the first time. Isomers 1 and 2 are the focal point of this work, and they have IUPAC names 1,2-didehydro-1,3-disilabicyclo[1.1.0]butane and silylsilacyclopropenylidene, respectively [70].
Optimizations were performed with no symmetry imposed, unlike in previous studies, and consequently many of the four-membered ring structures adopted puckered configurations. This puckering caused insignificant reductions in relative isomer energies accompanied by significant increases in dipole moments. Dipole moments computed using other levels were found to be within   Table 2 collects relative isomer energies for the 10 structures in Figure 1 using the Def2-QZVP basis set and wave-function methods including MBPT at second-, third-, and fourth-orders, as well as CCSD and CCSD(T). All calculations were performed on structures optimized at the MBPT (2) Table 2, were used to generate mean signed errors (MSE) and mean unsigned errors (MUE) for each method. Only MBPT(2) failed to achieve chemical accuracy, i.e., a MUE under 1.0 kcal/mol. Table 3 focuses on DFT-based methods, with pure GGAs, hybrid GGAs, and double-hybrid methods based on BLYP, PBE, and Minnesota XCFs represented. At the BLYP and B3LYP levels, isomer 2 was predicted as the global minimum, consistent with the small-basis BLYP results in Ref. [2]. All other methods tested here ordered isomer 1 below 2. The M06-2X and B2-PLYP+D3 approaches ordered all ten energies correctly, while DSDPBEP86 was also very close, erring by only 0.1 kcal/mol in ordering isomers 8 and 9. Only the DH-DFT methods methods produced MUEs within chemical accuracy, on average. Table 2. Statistical evaluation of wave-function methods for determination of the relative energies of low-lying isomers of C 2 Si 2 H 4 . Isomer energies were computed as ∆E = E(2) − E(1). All values are in kcal/mol.  Figure 2 collects temperature-dependent free-energy profiles computed using MBPT(2), B3LYP, and M06-2X in conjunction with the Def2-QZVP basis set. The study by Maier et al. described sample preparation by flash pyrolysis at 1500 K, and at this temperature all three methods place isomer 2 lowest in energy. It is thus likely that thermodynamic equilibrium of the structures was not achieved within the 10K Ar matrix preceding the IR measurements. Taking into consideration their temperature-independent enthalpic energy displacements, the free-energy profiles for each of the higher-lying isomers in Figure 2 are qualitatively very similar. The notable exception is isomer 4, which has a positive gradient when computed with MBPT(2) and a negative gradient when computed with B3LYP and M06-2X.

Infrared Spectroscopy of Low-Lying C 2 Si 2 H 4 Isomers
In Section 3.1, high-accuracy ab initio methods were employed to determine isomer 1 as the zero-temperature global PES minimum, while free-energy simulations flipped the order of the two lowest-lying isomers at the pyrolysis temperature reported by Maier et al. Assuming that a non-equilibrium mixture of isomers was present during the IR characterization performed at 10K, traces of isomer 1 may have been overlapped with peaks attributed to the dominant species, isomer 2. Broad bands of this type were noted by Maier et al. occurring at 2190 cm −1 and 919 cm −1 . If these peaks are attributable to isomer 1, this is supporting evidence that both isomers were measured in the non-equilibrium mixture, and helps further reconcile the past and present global minimum assignments. Figure 3 compares simulated IR spectra with peaks reported in measurements by Maier et al. Simulated spectra were computed for isomer 1 using the CASSCF(6,6), MBPT(2), and B3LYP levels of theory with the Def2-TZVP basis set. All theoretical spectra exhibit large-intensity frequencies in good agreement with the positions of the measured broad features. Vibrational simulations are known to be insensitive to both basis-set size and correlation effects, so higher-level results were not pursued. Thus, the theoretical vibrational characterizations of isomer 1 qualitatively match the broad peaks at 2190 and 920 cm −1 , and this warrants further analysis. Table 4 presents a quantitative comparison of the strongest IR peaks of isomers 1 and 2. Of the three methods tested, B3LYP produces the most accurate frequencies as compared to the measured IR peaks, with all values being within 50 cm −1 of an observed frequency. Indeed, the spectral fingerprints of the two isomers are quantitatively very similar, supporting the hypothesis that both isomers were present in the original measurements.

Photoisomerization of Low-Lying C 2 Si 2 H 4 Isomers
The experimental analysis by Maier et al. also inferred the presence of isomer 2 by monitoring photolysis products via IR spectroscopy. Irradiation with 313 or 366 nm light was found to isomerize 2 via ring-opening to form (silylethynyl)silylene (H 3 SiCCSiH) and ethynylsilylsilylene (H 3 SiSiCCH). Explicit mapping of isomerization mechanisms can be performed using intrinsic reaction coordinate (IRC) calculations to connect key stationary points on the PES. Photoisomerization mechanisms can be explored by performing vertical excitation energy (VEE) calculations on the ground-state IRC structures, though we note that this does not provide definitive proof of the mechanism since, in general, none of the considered structures are well-defined stationary points on the excited-state surfaces. Figure 4 compares potential energy cuts for the ring-opening of isomer 2, as generated using a variety of wave-function and DFT methods. The underlying ground-state structures were generated from IRC calculations connecting isomer 2 to the observed non-cyclic product H 3 SiCCSiH. IRC calculations were performed using the CASSCF(6,6), B3LYP and M06-2X methods, in Figure 4a-c, respectively. Within each frame, other methods were applied on top of these underlying geometries.  Experimental evidence suggests that irradiation of the reactant structure leads to photoconversion, implying a negative gradient on the excited-state surface at the equilibrium geometry. All theories generated potential cuts consistent with this mechanism. Evidently the second excited state is accessed by absorption of an incident photon and the reaction proceeds by adiabatic relaxation along a negative potential gradient, eventually transitioning through a cascade of avoided crossings to form the product isomer. The wave-function and DFT results differ significantly only in the placement of the avoided crossing between the first-and second-excited states, with the EOMCC methods placing it near the reactant geometry and TDDFT methods placing it near the transition state geometry. Figure 5 collects stationary points characterized using B3LYP/DZP and carefully connected by IRC calculations. Corresponding relative energies are also provided, though, due to the basis set size, discussion is reserved to qualitative aspects of the surface. Isomer 1 was found to be connected to isomer 2 via several transition states and intermediates involving transannular hydrogen migration and bond rearrangements. The largest barrier height separating isomers 1 and 2 is shown to be within 1-2 kcal/mol of that separating isomer 2 and the second observed photolysis product, ethynylsilylsilylene. Thus, the activation energy required to surmount barriers separating isomer 2 and products will also drive the reaction converting isomer 1 to products. The magnitude of the barrier heights excludes the possibility of thermal isomerization in normal laboratory conditions. Maier et al. reported photoisomerization after irradiation at 313 or 366 nm, resulting in the observed non-cyclic photolysis products H 3 SiCCSiH and H 3 SiSiCCH. To confirm that these frequencies will drive photoconversion of both 1 and 2, VEEs were computed at each stationary point shown in Figure 5. Table 5 collects the computed relative ground-state energies and VEEs for each stationary point structure, with labeling preserved from Figure 1. The previously reported irradiation frequencies, converted to 91.3 and 78.1 kcal/mol, are within 10% of the predicted VEEs for isomers 1, 6, and 3. A spectrally broadened light source may photoexcite and interconvert between each of the considered cyclic isomers 1, 2, 3, and 6, eventually driving the photolysis to the observed non-cyclic products H 3 SiCCSiH and H 3 SiSiCCH. Table 5. Relative energies (above) and excitation energies (below) for the C 2 Si 2 H 4 stationary points appearing in Figure 5. All values are reported in kcal/mol.

Conclusions
Ten C 2 Si 2 H 4 isomers were characterized as PES minima and energy-ordered using high-level coupled-cluster calculations. Isomer 1, or 1,2-didehydro-1,3-disilabicyclo[1.1.0]butane, was assigned as the global minimum, in contrast to all existing literature on this system which favored isomer 2, or silylsilacyclopropenylidene. High-level coupled-cluster isomer energies were used to benchmark DFT methods, and BLYP and B3LYP incorrectly predicted isomer 2 as the global minimum, even when very large basis sets were employed. Many DFT methods struggled to order all ten isomer energies correctly with respect to benchmarks. Only the B2-PLYP+D3 DH-DFT method energy-ordered all isomers similarly to the benchmarks, while also producing an MUE within chemical accuracy. Among the pure and hybrid functionals tested, only M06-2X predicted the same energy ordering with respect to the benchmarks but its MUE was slightly above the chemical accuracy threshold.
Assignment of isomer 1 as the global minimum required a reinterpretation of prior experimental data to reconcile with past assignments of isomer 2 as the global minimum. Computations helped show that the frequencies of isomers 1 and 2 overlap, which likely contributed to the broad bands previously observed in the IR difference spectrum. Some PES mapping was also performed to provided insight into the mechanism of photolysis of isomers 1 and 2 to form the observed products H 3 SiCCSiH and H 3 SiSiCCH. These simulations suggested that the ultraviolet frequencies used to irradiate the sample may initiate photolysis in both isomers 1 and 2, leading to the observed products.
Taken together, these simulations provide compelling evidence that the energetically similar isomers 1 and 2 were likely both present as reactants in the measurements reported in Ref. [2]. In light of the present analysis, there exists no contradictory evidence in the literature to dismiss assignment of isomer 1 as the global minimum. In closing, isomer 1 has been identified as a prominent configuration on the global PES of a barrierless acetylene-disilyne collision, and its dipole moment is strong enough for electron binding. Thus, it may be of significance in interstellar silicon carbide reaction networks. This prospect will be investigated in a future study.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: