Prediction Through Quantum Dynamics Simulations: Photo-excited Cyclobutanone

Quantum dynamics simulations are becoming a standard tool for simulating photo-excited molecular systems involving a manifold of coupled states, known as non-adiabatic dynamics. While these simulations have had many successes in explaining experiments and giving details of non-adiabatic transitions, the question remains as to their predictive power. In this work, we present a set of quantum dynamics simulations on cyclobutanone, using both grid-based multi-configuration time-dependent Hartree (MCTDH) and direct dynamics variational multi-configuration Gaussian (DD-vMCG) methods. The former used a parameterised vibronic coupling model Hamiltonian and the latter generated the potential energy surfaces on-the-fly. The results give a picture of the non-adiabatic behaviour of this molecule and were used to calculate the signal from a gas-phase ultrafast electron diffraction (GUED) experiment. Corresponding experimental results will be obtained and presented at a later stage for comparison to test the predictive power of the methods. The results show that over the first 500 fs after photo-excitation to the S$_2$ state, cyclobutanone relaxes quickly to the S$_1$ state, but only a small population relaxes further to the S$_0$ state. No significant transfer of population to the triplet manifold is found. It is predicted that the GUED experiments over this time scale will see s signal related mostly to the C-O stretch motion and elongation of the molecular ring along the C-C-O axis.


I. INTRODUCTION
Molecular Quantum dynamics (QD) has become a powerful computational tool for understanding fundamental reactivity.By solving the time-dependent Schrödinger equation these simulations can follow the time evolution of nuclei after a molecular system is prepared in a particular way, for example in a molecular beam or pump-probe experiment.The key feature is that the quantum nature of the nuclei is taken into account, which is essential for systems involving either tunnelling or non-adiabatic effects.
A present question, that is to be addressed in this paper, is whether QD simulations have predictive power, which is necessary for them to become a mature and trusted tool.This work addresses a challenge issued last year and aims to predict outcomes of experiments on cyclobutanone conducted at the SLAC Megaelectronvolt Ultrafast Electron Diffraction facility.These experiments commenced in late January 2024.The molecules will be photoexcited at 200 nm with a short (approximately 100 fs wide) pulse and diffraction images will be recorded over many picoseconds.Up to 200 fs, images will be recorded every 30 fs, and then the delay time will be increased.The challenge to simulation is to produce signals that can be compared to the experimental outputs and aid in their interpretation.Ultrafast electron diffraction is proving its value by providing data on structural dynamics in photochemistry 1,2 .Examples are the dissociation of water 3 , the ring-opening of cyclohexadiene 4 and imaging the passage of CH 3 I as it passes through a conical intersection 5 .
Cyclobutanones are an interesting target as they exhibit distinct chemical reactivity compared to cyclic ketones with larger rings, primarily owing to their inherent ring strain of around 105-120 kJmol −1 . 6,7Alongside the considerable ring strain, there is a heightened electrophilicity of the carbonyl carbon atom, a ring puckering induced by steric interaction of substituents at C2 and C4, and significant photochemical reactivity.Despite this, cyclobutanones were long considered as academic curiosities. 80][11][12] Reviews of relatively standard unsubstituted cyclobutanone reactions can be found in 6,8 , and, for example, in a more recent publication it has been reported that an isonitrile-based, four-component Ugi reaction, 13 utilising cyclobutanone, effectively produced an aspartame analogue. 14he earliest synthetic protocol for the production of cyclobutanone, from cyclobutanecarboxylic acid, authored by N. Kirscher 15 , suffered from inefficiency, yielded low quantities and featured several sequential reaction steps.Since then, several high-yield synthetic methods have been devised.One synthetic scheme involves the epoxidation of methylidenecyclopropane to 1-oxaspiro [2.2]pentane, leading to a subsequent lithium-catalyzed rearrangement. 16An alternative synthetic scheme entails a dialkylation process (utilizing 1-bromo-3-chloropropane) of 1,3-dithiane.This is followed by deprotection of the ketone through treatment with a mercury salt and cadmium carbonate. 17ere has been extensive research into the photo-induced behaviours of cyclobutanone, and it's derivatives, both theoretically and experimentally.An early study of cyclobutanone showed two channels to photolytic decay; the C2 channel, whereby ketene and ethylene are produced, and the C3 channel, whereby cyclopropane/propene and CO are produced, with a ratio of 2:3. 18rther theoretical and experimental studies have shown that the ratio of photoproducts is dependant upon the wavelength of excitation.Upon excitation between 340-240nm the (nπ * ) S 1 state is accessed, with the absorption maximum at around 280nm. 19,20Upon excitation to this state, three vibrational modes are activated corresponding to CO stretching, CO out-of-plane wagging and ring puckering. 21Studies have indicated that photoexcitation at shorter wavelengths (<315nm range) an IC process via α-cleavage occurs resulting in the production of photofragments with a ratio (C2:C3) of approximately 2:1. 22,23In the 345-315nm range, an ISC process occurs via a triplet state, followed by the production of the photofragments, with a ratio (C2:C3) of approximately 1:2 at 326.3 nm 23,24 and approximately 1:7 at 343.7 nm. 23This indicates an activation barrier to the breaking of the CC bond, a ring opening process, which has been found to be between 9.6 kJmol −121 and 29 kJmol −1 . 25,26By comparison, the barrier to ring opening of cyclopentanone and cyclohexanone is in excess of 63 kJmol −1 . 21Experimentally, the formation of hot ketene fragments from cyclobutanone, in a cyclohexane solution occurred, at the time delays of 0.25 ps, indicating an ultrafast ring-opening pathway.The subsequent relaxation of the vibrationally hot photoproducts and the growth of the ketene fundamental band were fitted to biexponential functions with shared time constants of 7 ps and 550 ps. 25 A theoretical study using the Ab-Initio Multiple Spawning method estimated the S 1 lifetime to be around 484 fs with an average time taken for the α-cleavage occurring in 176.6 fs. 27It should be noted that theoretical findings are in the gas phase.
Experimental and theoretical studies have also explored the second excited (S 2 ) state of cyclobutanone, and other cyclic ketones.In the Franck-Condon region the S 2 state exhibits Rydberg 3s character and undergoes an IC process, facilitated by the low frequency ring puckering mode, to the S 1 (nπ * ) state. 28The rate at which this 3s state decays was found, experimentally, to be around 0.74 ps.A complementary theoretical study, running quantum dynamics on potentials generated from a linear vibronic coupling Hamiltonian using five degrees of freedom, found the decay of the 3s state to be around 0.95 ps . 29,30It was also shown that in addition to the ring puckering mode, the CO out-of-plane deformation was significant.By comparison, the experimentally obtained ratio of the relative rates of the 3s→nπ * process in cyclobutanone, cyclopentanone and cyclohexanone were determined to be 13:2:1. 28e calculations presented here will use a standard work-flow for the study of a photoexcited molecule, in this case cyclobutanone, with different techniques building up a picture of the molecular dynamics.The work-flow has four stages.The first is to choose an appropriate level of quantum chemistry to describe the excited-state potential energy surfaces and couplings.The second is to build a model Hamiltonian using the vibronic-coupling scheme.This will provide information on which vibrational modes are excited on photo-excitation and which modes provide non-adiabatic coupling between the electronic states.The third step is to run grid-based QD simulations using the Multi-Configuration Time-Dependent Hartree (MCTDH) method.This will provide an accurate description of the short-time dynamics, along with absorption spectra.Finally, direct QD will be performed using the variational Multi-Configuration Gaussian (vMCG) method.In these calculations the potential energy surfaces are calculated on-the-fly, allowing the molecule to undergo long-range motion such as fragmentation, which is not possible in the MCTDH calculations due to the nature of the model Hamiltonian.The vMCG description also has an underlying trajectory nature which will enable a simulation of the experimental signal and relate it to the evolving molecular geometry.
The MCTDH and DD-vMCG simulations used a development version of the Quantics Suite 34,35 , with the Vibronic Coupling Hamiltonians built using the VCHam package 36

A. Quantum Chemistry Calculations
The starting point is to choose the level of electronic structure theory able to describe the system of interest, here cyclobutanone, balancing accuracy with cost.We start by finding and characterising the ground-state equilibrium structure.For this the coupledcluster singles and doubles (CCSD) method was used with a 6-311++G** basis sets.The next step is to choose the level of theory needed to describe the excited-states.To describe potential long range motion, a multi-configurational method is best and as cyclobutanone is quite small with only five heavy atoms, a complete active space self-consistent filed (CASSF) wavefunction is possible.After trying a number of different complete active space (CAS) sets of orbitals, it was decided that a CAS with 8 electrons in 10 orbitals would be suitable.The orbitals include the π orbitals and low lying Rydberg orbitals.The CAS orbitals are shown in Fig. 1.After further tests calculating excitation energies at the C 2v structure, it was decided that state-averaging over six states and the large, 6-311++G** basis set was required for stable results.Additionally, the inclusion of a second-order perturbation theory correction (CASPT2), using the RS2C method implemented in Molpro, was used to provide improved energies.

B. Vibronic Coupling Model Hamiltonian
Assuming a diabatic electronic basis, the Hamiltonian can be written in matrix form as where T is the kinetic energy operator, 1 a unit matrix and W the diabatic potential matrix.
Thus in this representation the kinetic energy operator is diagonal in the electronic states and the nuclear-electronic, non-adiabatic, coupling is represented by potential-like functions as the off-diagonal elements of the potential matrix.The key feature is that, unlike the usual adiabatic representation, the non-adiabatic couplings do not contain singularities 37,38 .
The Vibronic Coupling Model represents the diabatic potentials using Taylor expansions around the Franck-Condon (FC) point and truncated to a low order.
In the standard Linear Coupling (LVC) scheme emplayed here, the zero-order Hamiltonian is taken to be the ground-state Hamiltonian in the Harmonic approximation.
Here, the coordinates used are mass-frequency scaled normal modes with ω α representing the ground-state frequencies.The zero-order diabatic matrix is a diagonal matrix with the energies of the states at the Franck-Condon point.The first-order expansion matrices contain the linear on-diagonal coupling parameters, κ α along with the off-diagonal parameters, λ (ij) α .The former are the gradients of the potential surfaces at the Franck-Condon point, while the latter are the non-adiabatic coupling between states i and j.Higher-order matrices contain only on-diagonal terms as required to model the diabatic surfaces, hence the name linear vibronic coupling.Terms up to fourth-order were required for some modes due to the anharmonicity of the surfaces.
The parameters are only non-zero if they obey the symmetry rule that the product of normal modes and states involved in the expansion term must contain the totally symmetric representation.For example, the linear on-diagonal terms in terms of the related Hamiltonian matrix element is where H el is the usual electronic Hamiltonian with ψ i , the electronic wavefunction at the Franck-Condon geometry, and the derivative is also evaluated at this point.This matrix element is only non-zero if Γ i × Γ i × Γ α ⊃ A 1 , where Γ i , Γ α are the irreps of the electronic states and normal mode, respectively, and A 1 the totally symmetric irrep.Thus only totally symmetric vibrations have non-zero κ parameters.In a similar way, λ parameters are only non-zero if Γ i × Γ j × Γ α ⊃ A 1 and thus only vibrations with the correct symmetry can couple states.These rules are easily extended to higher order terms.

C. Multi-Configuration Time-Dependent Hartree (MCTDH)
The MCTDH method provides a variational solution to the time dependent Schrödinger equation (TDSE) using the wavefunction ansatz i.e. a full direct product expansion of the wavefunction in terms of p sets of low-dimensional basis functions, φ i (q i , t), known as single-particle functions (SPFs).These functions are time-dependent and depend on a set of physical coordinates which provides an underlying grid.The functions |s⟩ are time-independent vectors that define the electronic state populated.This is usually referred to as the single-set formulation of MCTDH as a single set of SPFs are used to describe all states.
The method is thus a contraction scheme from the full basis set to a variational one.
The resulting equations of motion for the time evolution of the SPFs and the expansion coefficients A j 1 ...jps are well described in the literature 39,40 and will be not be given here.
For large systems, the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) variant must be used [41][42][43] .In this, the SPFs are expanded in the form of the MCTDH ansatz Eq. ( 8).These new basis functions (SPFs of the first layer) can in turn be expanded in this form.This procedure is continued to provide layers of functions with a set of primitive grid functions forming the lowest layer.In this way, a full tensor contraction scheme is set up and used to variationally solve the TDSE.
The ML-MCTDH method was used in this work to allow simulations including all degrees of freedom.The Hamiltonian is provided by the vibronic model described above.The primitive basis functions used for all coordinates were harmonic oscillator discrete variable representations (DVRs) 39 .Calculations need to be converged with respect to the basis sets.
In order to achieve this, the size of the primitive basis is checked to ensure the wavepacket does not significantly populate the end grid points (population less than 10 −6 ).The lowest natural populations of the SPFs are also kept below 10 −3 .The latter was ensured throughout the calculation by dynamically growing the basis in any layer whenever the population of the least important SPF approaches this limit 44 .
The MCTDH calculations using the vibronic coupling Hamiltonian are able to simulate the absorption spectrum from the Fourier Transform of the autocorrelation function The function g(t) is used to remove artefacts from the Fourier Transform.The chosen form is where the first term is to ensure that the autocorrelation function goes to zero at the end of the simulation, time T , and the second term a damping function.The damping time used was 150 fs.
In addition, diabatic state populations are directly obtained from the wavefunction which give a rate of electronic relaxation as the molecule changes electronic configuration.The short-time dynamics may also be analysed for the evolution of the molecular geometry in terms of excitation of the vibrational coordinates.

D. Direct Dynamics variational Multi-Configuration Gaussian (DD-vMCG)
To create a more flexible wavefunction that has a direct connection to molecular structures, the wavefunction ansatz can be written as a superposition of Gaussian functions with each Gaussian function a separable product of one-dimensional Gaussian functions j (q 1 , t) . . .g having the form where ζ jκ , ξ jκ and η j are quadratic, linear and scalar parameters respectively.
Solving the TDSE using this ansatz, and the Dirac-Frenkel variational principle, leads to the variational Multiconfigurational Gaussian (vMCG) method 45,46 .The expansion coefficients evolve with equations similar to those for MCTDH.However, the non-orthogonality of the GWP basis must be taken into account with the overlap matrix, S, Hamiltonian matrix, H and overlap time-derivative, τ In the standard approach, the widths of the Gaussians, ζ, are kept fixed (frozen Gaussians) and the scalar parameter, η is fixed by the requirement for the Gaussians to be normalised and for the phase to be kept zero.The time-dependence is then carried by the linear parameters, ξ.Collecting the set of parameters for a multi-dimensional Gaussian into a vector, Λ jα = ξ jα , the equations of motion (EOM) for the Gaussian functions can be written as Using the relationship between the general form of Eq. 14 and the Gaussian wavepackets of Heller 47 , the linear parameters can be written in terms of the coordinate and momentum at the centre of the Gaussian Λiα = ξiα = −2ζ i qα + i ṗα (20)   and the vector X in Eq. ( 19) is related to the classical equations of motion Thus the Gaussians move along trajectories that have a classical component with an additional variational coupling between leading to faster convergence of the wavefunction.
Due to the localised nature of a Gaussian function, it is reasonable to calculate integrals of the potential energy using a Local Harmonic Approximation (LHA) in which the potential is expanded to second order around the centre of a Gaussian Solutions to the TDSE no longer converge on the exact result, but the integrals can all be performed analytically and the algorithm can be used for direct dynamics simulations.
Direct dynamics simulations use this formalism to calculate the potential from quantum chemistry calculations on-the-fly.The LHA requires the energies, gradients and Hessians at the Gaussian central coordinate, and these can all be provided by a quantum chemistry calculation.In the standard DD-vMCG protocol, these quantities are calculated and stored in a database.New points are added to the database only when a GWP has coordinates that are significantly different from any point stored in the database.The surfaces experienced by the evolving GWPs are provided by Shepard interpolation between the points in the database.Calculations are run in the diabatic picture, and the potentials are diabatised using the propagation diabatisation scheme.For full details see Ref. 48 .
The major effort in the vMCG method is the inversion of the C matrix in Eq. ( 19) which has the size of (N f × n) 2 , where N f is the number of degrees of freedom and n the number of basis functions.This effort can be reduced by partitioning the system and making the wavefunction a multi-configurational product from two (or more) sets of functions.
This is the G-MCTDH ansatz 45 .The EOM are the same as Eqs.(15, 19) but now the function G i is a configuration, i.e. a product of basis functions, G i = G (1) i 2 . . .and the propagation of the Gaussian parameters include mean-field operators connecting the different partitions.

E. Gas Phase Ultrafast Electron Diffraction (GUED) Signal
GUED uses high energy pulses of electrons to provide structural information in the form of a diffraction pattern.Electron scattering by a potential field (a molecule) can be described by measuring the momentum transfer from an electron to the target on collision.For a set of atoms under the independent atom model, this can be formulated to give an intensity as a function of scattered momentum transfer s that has two parts: an "atomic" due to scattering directly off an atom and a "molecular" due to interference from scattering off neighbouring atoms 49 .The atomic part is and the molecular part is where N is the number of atoms, f i are the "atomic form factors" accounting for the structure of an atom in terms of nuclei and electrons and r ij denotes the inter-atomic distance between the i th and j th atom.The total intensity is simply the sum of these two terms.Actual measurements record the scattering signal as a modified scattering intensity sM (s, t) = s I mol (s, t) I at (s) (26)  where the time dependence has been added to indicate that the signal is obtained by scattering off the sample at time t.It is more useful to record the signal as the difference between the signal at t and an initial (reference) signal ∆sM (s, t) = s ∆I mol (s, t) where α is a smoothing factor.Changes in the molecular geometry can thus be seen in the PDF as losses and rises as atoms move apart.It is, however, summed over all the pairs in a molecule so each inter-nuclear distance may be related to a number of different atom pairs, making interpretation with a molecular simulation difficult.From a simulation, the PDF can be simulated by using the atomic coordinates along a trajectory to get first the molecular intensity, Eq. ( 25), and then obtain the modified scattering intensity before calculating the difference PDF of Eq. ( 28).
Finally, it needs to be taken into account that a molecule is not a classical object moving along a single trajectory.The scattering takes place from the molecular density given by the evolving wavepacket.This can be taken into account by integrating over the molecular configurations where Ψ 2 (q) is the probability of the molecule being in configuration q.
In vMCG, a useful approximation for calculating the expectation values of operators that depend only on the atomic coordinates is where r j is the molecular configuration defined by the centre coordinate of the GWP, and O(q j ) the value of the operator at that point.The expectation value is thus approximated by a weighted sum along the GWP trajectories, with the weighting given by the Gross Gaussian Populations where S ij is the GWP overlap matrix.The final difference PDF from the signal calculated over the molecular density can thus be approximated in this way as a weighted sum where I mol,i is the molecular scattering intensity from vMCG trajectory i.In fact, this can be carried through the procedure for calculating the difference PDF and one can provide simply a weighted sum of difference PDFs.
The modified scattering intensity for each trajectory at each time step was calculated using a code to simulate electron diffraction signals provided by Wolf et al. 50  of 0.0087 eV.The frequency of the transition mode, ν 1 means it has a zero point energy level at 0.0062 eV.Thus even though the barrier is very small, the low frequency of the transition mode means that it is a stable structure and at 0 K the equilibrium structure will be C s .However, at room temperature, a significant population of molecules will be in excited vibrational levels and the ground-state structure will be closer to C 2v .It was also found that at the CASSCF and CASPT2 levels of theory used in the dynamics calculations the C 2v structure is the minimum energy on the ground-state.We therefore will need to examine the dynamics starting from both structures.

B. Excitation Energies and Absorption Spectra
The Excitation energies and characters from both structures to the lowest two singlet and three triplet states are listed in Table I.Calculations used a CAS with 8 electrons in 10 orbitals and a 6-311++G** basis.The orbitals are shown in Fig. 1.
Little difference is found for the excitation from the two structures.In agreement with previous work 21,26,30 , the lowest singlet state is an excitation from the in-plane π * y to the π * orbital at around 4.2 eV.The second state is an excitation out of the same orbital into a diffuse Rydberg orbital lying around 6.3 eV.Interestingly the CASPT2 correction shifts this state up in energy.The T 1 and T 2 states are the same configurations as S 1 and S 2 , but lie slightly lower in energy.The T 3 state lies above S 2 at 6.7 eV and is a π to π * excitation.

C. Vibronic Model Hamiltonian
The vibronic coupling model Hamiltonian of Sec.II B was obtained using the massfrequency scaled normal modes of the ground-state C 2v structure, with the imaginary frequency of ν 1 taken as real.The ground-state and excitation energies were calculated at the CASPT2 corrected SA6-CAS(8,10)/6-311++G** level of theory at twenty one points along each vibrational mode, treating the singlet and triplet manifolds separately.The parameters of the model were then chosen to optimise the fit between the adiabatic potential surfaces of the model and the calculated adiabatic energies 52 .Due to the anharmonicities of the potentials, quartic potentials (expansion to fourth order) were used for all modes.The diabatic coupling, however, was truncated to first order (linear vibronic coupling).
The singlet and triplet manifolds were then merged into one operator, and the spinorbit couplings calculated at the Franck-Condon point added as constant values between the singlet and triplet states, using the magnitude of the complex components to give an effective coupling treating each triplet state as a single state.The excitation energies are given in Table I.The spin-orbit couplings, given as a normalised average over the three components of the triplet state are given in Table II.The model provides information on the key vibrational modes: those that are directly vibrationally excited on electronic excitation and those that provide coupling between the states.Tables III and IV list the linear κ and λ parameters.Coupling strengths are defined as κ α /ω α , where ω α is the frequency of the mode.Using these strengths as a criteria, a subset of 12 key modes can be identified.ν 5 , ν 16 , ν 19 , ν 20 and ν 21 are the totally symmetric vibrations with the highest κ coupling strengths, i.e. they have significant gradients on the potential surfaces at the Franck-Condon point.Modes ν 4 and ν 9 , with A 2 symmetry carry the coupling between S 0 and S 1 .Modes ν 3 and ν 17 , with B 2 symmetry couple between S 0 and S 2 , while modes ν 2 and ν 25 , with B 1 symmetry carry the coupling between S 1 and S 2 .Finally, ν 1 has symmetry B 1 and dominates the coupling between T 1 and T 2 .The coupling between T 1 and T 3 is dominated by ν 4 , and between T 2 and T 3 by ν 11 .This latter coupling though is quite weak and ν 11 is not included in the key modes as the high energy T 3 state is unlikely to be important for the short term dynamics.The key vibrations are plotted in Fig. 3.
The parameters for the full final Hamiltonian are given as the Quantics operator file in the supplementary datasets.Cuts through the potential surfaces for the full singlet and triple manifold are also given in the SI.In the triplet manifold, the major features are due to the near degenerate T 2 and T 3 states, which are also near degenerate with the S 2 state.
The spin-orbit coupling is however weak.In this manifold, ν 2 , ν 3 and ν 16 show double well structures in T 2 , while in ν 21 the T 2 and T 3 states cross at the Franck-Condon point and T 2 looks to be dissociative.ν 21 is the C-O stretch vibration and is also the most interesting mode in the sing manifold.The S 2 / S 1 crossing is seen to negative values along this modes and in the downhill direction from the FC point,.The S 1 state is also long and fairly flat, but not dissociative, out to positive values.The mode ν 5 , which is the ring stretching along Mode Symmetry κ (2) (eV) κ (3) (eV) κ (4) (eV) κ (5) (eV) κ (6) (eV)

D. Importance of Triplet States
ML-MCTDH simulations were run for 200 fs on the full singlet and triplet vibronic coupling operator of Sec.III C. The tree defining the layering along with the number of functions is given in the SI.Simulations started in both the C 2v and C s structures.At the end of this, less than 0.1 % population flows into the T 1 state, and the populations of the T 2 and T 3 are even smaller.This is not surprising given the small size of the coupling -for both the S 1 and S 2 states there is zero SOC with the close lying triplet (T 1 and T 2 , respectively).
Of course the vibronic coupling model treats the SOC with only a constant value taken at the Frank-Condon point.To see if more ISC occurs if the nuclear motion is taken into account, DD-vMCG simulations were run including all 6 states and calculating the SOC at each point added to the DB.All properties were calculated at the CAS(8,10) level.The wavepacket was partitioned into 2 parts with the 12 key modes in one partition and treated with 16 GWPs and the remaining modes in the second partition treated with 8 GWPs.This is a small basis set but with 128 trajectories allows a reasonable searching of the main  Mode Symmetry λ (1,2) /ω (eV) λ (1,3) (eV) λ (2,3) (eV) λ (4,5) (eV) λ (4,6) (eV) λ (5,6) (eV) At the end of the simulation, all three triplet states had accumulated approximately 0.1 % population.Thus, including the nuclear motion has indeed increased the ISC, but it is still not likely to be significant for the short-term dynamics (over the first 0.5 ps).As a result, the triplet states will be ignored in the final analysis.The experimental spectrum for the S 1 band is a broad, featureless band from 3.8 eV to 5 eV with the maximum at around 4.2 eV 21 .This is a dark state and the intensity is low.
The calculated spectrum from the C 2v structure is found to be in the correct energy region, but the maximum is too low in energy by 0.4 eV, the band is too narrow and too structured.
Starting from the C s structure, the structure is much reduced, but the band is too broad and high in energy.
The experimental S 2 band is again broad, but with a weak progression of peaks.It starts at 6 eV and runs to 6.8 eV, with a maximum at 6.5 eV 21 .The calculated spectrum from the C 2v structure lies in the correct energy region, perhaps 0.2 eV too low, and has too much structure.The spectrum from the C s structure is again too broad and high in energy.
The spectra indicate that the spectra are better reproduced by excitation from the C 2v structure even though this is a transition state.This is, however, not a surprise due to the low barrier that means that at room temperature the ground state wavefunction will not be purely in the ground-state C s minimum and will spread across the equivalent structures to be more C 2v like.

IV. PREDICTION AND CONCLUSION
We have performed quantum dynamics simulations to predict the dynamics of cyclobutanone as it will be observed by a gas-phase ultrafast electron diffraction experiment over 500 fs after initiating the dynamics by excitation to the S 2 electronic state.Cyclobutanone is a non-trivial molecule for simulations.It been recently shown that the population dy- namics for the related cyclopropanone molecule are very sensitive to the electronic structure method used 55 .The fact that the ground-state equilibrium structure is a shallow well, and the molecule easily converts between puckered structures through a planar transition state means that choosing the initial conditions is also not straightforward as the geometry will be temperature dependent.The potential involvement of triplet states brings a further challenge.
The final results from the simulations indicate that triplet states are not significantly populated over the first 500 fs.However, it is also found that while the initial relaxation from S 2 is fast, taking place in around 50 fs, subsequent relaxation from S 1 is slow.Thus at later time-scales there will be crossing to the triplet manifold and this can result in photofragmentation 56 .Our simulations also show that in the singlet manifold over 500 fs the molecule does little more than vibrate along the central C-C-O axis.The C-O moiety is
) as the change in signal gives the structural changes in the molecule without needing to extract I mol .This signal can be related to the change in "pair distribution function" (PDF) by integrating over s to give the signal as a function of inter-atomic distances ∆P (r, t) ≈ r smax s min ∆sM (s, t) sin(sr)e −αs 2 ds.( with subsequent weighting according to Eq.(33) to get the final time-dependent signal.The time-dependent difference pair distribution function ∆P (r, t) was then obtained using s min = 0 Å−1 , s max = 10 Å−1 and α = 0.036 Å2 matching the experimental spatial resolution of 0.6 Å.III.RESULTSA.Ground-State Equilibrium StructureUsing CCSD/6-311++G**, two stable structures of cyclobutanone were found.The first is planar with C 2v symmetry.The second is bent and has C s symmetry.These are shown in Fig.2and the coordinates, along with bond lengths and angles, are given in the Supplementary Information.The ring puckering angle between the C 1 C 2 C 3 and C 4 C 2 C 3 planes was found to be 171.7 • , and the angle of the C-O bond to the C 1 C 2 C 3 plane was 12.3 • .Frequency calculations showed that the C s structure is a minimum and the C 2v structure a transition state.The frequencies of the vibrational modes of both structures are listed in the Supplementary Information.With one exception, the transition mode, changes in the normal modes and frequencies are minimal between the two structures.The transition mode with an imaginary frequency of 80 cm −1 is the ring puckering motion, ν 1 and at C s it has a real frequency of 100 cm −1 .The difference in energy between the structures, ∆E = E C 2v − E Cs is the barrier height to the out-of-plane bend and has a value

FIG. 2 .
FIG. 2. The stable ground-state structures of cyclobutanone with the atom numbering used throughout the paper.The angle shown on the C s structure is the ring-puckering angle calculated at the CCSD/6-311++G** level of theory.

FIG. 3 .
FIG. 3. The 12 key normal modes of cyclobutanone, calculated at CCSD/6-311++G** level on the Ground-state C 2v structure.Arrows indicate the motion of the atoms during each vibration.

E 4 .FIG. 4 .
FIG. 4. State populations from MCTDH calculation on the singlet manifold vibronic coupling model of cyclobutanone.Calculations (a), (c) start with the C 2v ground-state structure while (b), (d) start from C s .(a), (b) include only the key 12 modes, while (c), (d) are the full 27D calculations.

FIG. 5 . 3 (
FIG. 5. Spectra from ML-MCTDH calculation on the singlet manifold vibronic coupling model of cyclobutanone.(a), (c) The S 1 spectrum starting from either the (a) C 2v or (c) C s structure.(b), (d) The S 2 spectrum starting from either the (b) C 2v or (d) C s structure.

50 FIG. 6 .FIG. 7 .
FIG. 6. Pair Distribution Function (PDF) from 512 DD-vMCG trajectories calculated on CASSCF potential surface.(a) Reference PDF for the initial structure.(b) Difference PDF as a function of time.The blue peaks denote loss with respect to the reference spectrum of (a), while red peaks denote gain