Prediction through quantum dynamics simulations: Photo-excited cyclobutanone Open Access

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 required for a complete description of systems involving either tunneling or non-adiabatic effects.

A present question, which 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 photo-excited at 200 nm with a short (∼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,³ the ring-opening of cyclohexadiene,⁴ and imaging the passage of CH₃I as it passes through a conical intersection.⁵ 

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 kJ mol⁻¹.^6,7 Alongside 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.⁸ Since the mid-1970s, the role of cyclobutanones as synthetic reactants and intermediates has expanded, enabling the production of a wide variety of compounds with diverse applications; a small selection can be seen in the following articles.^9–12 Reviews of relatively standard unsubstituted cyclobutanone reactions can be found in Refs. 6 and 8; for example, in a more recent publication, it has been reported that an isonitrile-based, four-component Ugi reaction,¹³ utilizing cyclobutanone, effectively produced an aspartame analog.¹⁴ 

The earliest synthetic protocol for the production of cyclobutanone, from cyclobutanecarboxylic acid, authored by N. Kirscher,¹⁵ 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 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.¹⁶ 

There has been extensive research into the photo-induced behaviors of cyclobutanone, and its 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.¹⁷ 

Further theoretical and experimental studies have shown that the ratio of photoproducts is dependent on the wavelength of excitation. Upon excitation between 340 and 240 nm, the (nπ*) S₁ state is accessed, with the absorption maximum at around 280 nm.^18,19 Upon excitation to this state, three vibrational modes are activated corresponding to CO stretching, CO out-of-plane wagging, and ring puckering.²⁰ Studies have indicated that following photoexcitation at shorter wavelengths (⁠$<$315 nm range) an IC process via α-cleavage occurs resulting in the production of photofragments with a ratio (C2:C3) of ∼2:1.^21,22 In the 345-315 nm range, an ISC process occurs via a triplet state, followed by the production of the photofragments, with a ratio (C2:C3) of ∼1:2 at 326.3 nm^22,23 and ∼1:7 at 343.7 nm.²² This indicates an activation barrier to the breaking of the CC bond, a ring-opening process, which has been found to be between 9.6 kJ mol⁻¹²⁰ and 29 kJ mol⁻¹.^24,25 By comparison, the barrier to ring opening of cyclopentanone and cyclohexanone is in excess of 63 kJ mol⁻¹.²⁰ Experimentally, the formation of hot ketene fragments from cyclobutanone in a cyclohexane solution occurred with a time delay 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 and 550 ps.²⁴ A theoretical study using the ab initio multiple spawning method estimated the S₁ lifetime to be around 484 fs with an average time taken for the α-cleavage occurring in 176.6 fs.²⁶ It should be noted that theoretical findings are in the gas phase.

Experimental and theoretical studies have also explored the second excited (S₂) state of cyclobutanone and other cyclic ketones. In the Franck–Condon region, the S₂ state exhibits Rydberg 3s character and undergoes an internal conversion (IC) process, facilitated by the low frequency ring puckering mode, to the S₁ (nπ*) state.²⁷ The 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.^28,29 It 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 was determined to be 13:2:1.²⁷ 

The calculations presented here will use a standard work-flow for the study of a photo-excited 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 straightforward simulation of the experimental signal and relate it to the evolving molecular geometry.

To provide a reference for the challenge, the original manuscript has been deposited on the arXiv at DOI 10.48550/arXiv.2402.09933.

All quantum chemistry calculations were performed using the Molpro 2022 package.^30–32 The MCTDH and DD-vMCG simulations used a development version of the Quantics Suite,^33,34 with the vibronic coupling Hamiltonians built using the VCHam package.³⁵ 

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 characterizing the ground-state equilibrium structure. For this, the coupled-cluster singles and doubles (CCSD) method was used with a 6-311++G** basis set.

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 field (CASSCF) wavefunction is possible. To select the complete active spaces (CASs), a number of spaces were tried: (6, 6), (6, 5), (8, 8), and (8, 10). State averaging over either five or six states was also investigated. The character of the states and their energies were compared to experimental results and previous calculations.^20,21,37 After this, it was decided that a CAS(8, 10) with eight electrons in ten orbitals was required.

The CAS orbitals include the π orbitals and low lying Rydberg orbitals and are shown in Fig. 1. After 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.

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 $Aj1…jps$ are well described in the literature^38,39 and will be not be given here.

For large systems, the Multi-Layer Multi-Configuration Time-Dependent Hartree (ML-MCTDH) variant must be used.^40–42 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 for 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).³⁸ 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 that the wavepacket does not significantly populate the end grid points (population less than 10⁻⁶). A measure of the quality of the time-dependent basis set is given by the natural populations. These are the eigenvalues of the reduced density matrix associated with each set of basis functions.^38,42 These eigenvalues are related to the importance of the functions, and if the lowest eigenvalue is small, adding further functions to the basis set will not significantly change the wavefunction. To keep a reasonable quality of wavefunction, the lowest natural populations are kept below 10⁻³. 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.⁴³ 

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 analyzed for the evolution of the molecular geometry in terms of excitation of the vibrational coordinates.

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 diabatized using the propagation diabatization scheme. For full details, see Ref. 47.

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.⁴⁹ 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 $smin=0Å−1$⁠, $smax=10Å−1$⁠, and $α=0.036Å2$ to match the experimental spatial resolution of $0.6Å$⁠.

The nuclear coordinates used in the simulations are mass-frequency scaled normal modes. As a result, displacements are dimensionless and have the scale that 1 unit of displacement stretches the normal mode to have a potential energy equal to the zero-point energy $(12ℏω)$⁠. In the following, distances are simply recorded as x units of displacement or when it is clear from the context that we are referring to a displacement, x units.

Using CCSD/6-311++G**, two structures of cyclobutanone were found at critical points on the ground-state potential surface. The first is planar with C_2v symmetry. The second is bent and has C_s symmetry. These are shown in Fig. 2, and the coordinates, along with bond lengths and angles, are given in the supplementary material. The atom numbering was chosen to emphasize the proximities of the carbon atoms to the oxygen. Consequently, atom 1 is the closest carbon to the oxygen, carbons 2 and 3 are the symmetrically equivalent central atoms, and carbon 4 is the furthest away.

The ring puckering angle between the C¹C²C³ and C⁴C²C³ planes was found to be 171.7°, and the angle of the C–O bond to the C¹C²C³ plane was 12.3°. Frequency calculations show that the C_s structure is a minimum and the C_2v structure is a transition state. The frequencies of the vibrational modes of both structures are listed in the supplementary material. 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⁻¹ is the ring puckering motion, ν₁, and at C_s, it has a real frequency of 100 cm⁻¹. The difference in energy between the structures $ΔE=EC2v−ECs$ is the barrier height to the out-of-plane bend and has a value of 0.0087 eV. The frequency of the transition mode, ν₁, means it has a zero point energy level at 0.0062 eV. Consequently, while 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 most probable 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.

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 eight electrons in ten 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,^20,25,29 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₁ and T₂ states are the same configurations as S₁ and S₂, but lie slightly lower in energy. The T₃ state lies above S₂ at 6.7 eV and is a π to π* excitation.

The vibronic coupling model Hamiltonian of Sec. II B was obtained using the mass-frequency scaled normal modes of the ground-state C_2v structure, with the imaginary frequency of ν₁ 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 21 points along each vibrational mode, treating the singlet and triplet manifolds separately. The parameters of the model were then chosen to optimize the fit between the adiabatic potential surfaces of the model and the calculated adiabatic energies.⁵¹ 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 spin–orbit couplings calculated at the Franck–Condon point were 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 manifold as a single state. The excitation energies are given in Table I. The spin–orbit couplings, given as a normalized 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 by the ratios $κα(i)/ωα$ and $λα(ij)/ωα$⁠, where ω_α is the frequency of the mode. Using these strengths as a criteria, a subset of 12 key modes can be identified. The ν₅, ν₁₆, ν₁₉, ν₂₀, and ν₂₁ modes 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 ν₄ and ν₉, with A₂ symmetry, carry the coupling between S₀ and S₁. Modes ν₃ and ν₁₇, with B₂ symmetry, couple between S₀ and S₂, while modes ν₂ and ν₂₅, with B₁ symmetry, carry the coupling between S₁ and S₂. Finally, ν₁ has symmetry B₁ and dominates the coupling between the T₁ and T₂ states. The coupling between T₁ and T₃ is dominated by ν₄ and between the T₂ and T₃ states by the ν₁₁ mode. This latter coupling is, however, quite weak. Consequently, ν₁₁ is not included in the key modes as the high energy T₃ 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 material. Cuts through the potential surfaces for the full singlet and triple manifold are also given in the supplementary material. In the triplet manifold, the major features are due to the near degenerate T₂ and T₃ states, which are also near degenerate with the S₂ state. The spin–orbit coupling, however, is weak. In this manifold, ν₂, ν₃, and ν₁₆ show double well structures in T₂, while in ν₂₁, the T₂ and T₃ states cross at the Franck–Condon point and T₂ appears to be dissociative. Mode ν₂₁, the C–O stretch vibration, is the most interesting mode in the singlet manifold. The S₂/S₁ crossing is seen to be at a negative value along this mode and in a downhill direction from the FC point. The S₁ state is also long and fairly flat out to positive values, but not dissociative. The ring stretching mode along the C–C–O axis, ν₅, is also potentially interesting as it has a strong gradient at the FC point in all excited 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 supplementary material. Two simulations were started with the initial wavepacket centered either at the C_2v transition state geometry or the C_s minimum energy structure. In both cases, the wavepacket was a separable product of Gaussian functions with a width of $1/2$ along each normal mode, i.e., the harmonic ground-state wavefunction in the mass-frequency scaled normal modes of the model. The wavepackets were then placed in the S₂ state in a vertical excitation. At the end of this, less than 0.1% population flows into the T₁ state, and the populations of the T₂ and T₃ are even smaller. This is not surprising, given the small size of the coupling. In both the S₁ and S₂ states, there is zero spin–orbit coupling with the closest lying triplet state (T₁ and T₂, respectively).

Of course, the vibronic coupling model treats the spin–orbit coupling with only a constant value taken at the Frank–Condon point. In order to discern the effect of the inclusion of nuclear motion on the spin–orbit coupling and the amount of ISC, DD-vMCG simulations were run, including all six states, calculating the spin–orbit coupling at each point added to the DB. All properties were calculated at the CAS(8,10) level. The wavepacket was partitioned into two parts with the 12 key modes in one partition and treated with 16 GWPs. The remaining modes in the second partition were treated with 8 GWPs. While this is a small basis set, 128 trajectories allows for a reasonable searching of the main configuration space to find any regions where ISC occurs.

At the end of the simulation, all three triplet states had accumulated ∼0.1% of the population. Consequently, while the inclusion of the nuclear motion has indeed increased the ISC, 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.

ML-MCTDH simulations were run for 200 fs on the singlet manifold vibronic coupling model of cyclobutanone. Two sets of simulations were carried out as in the simulations of Sec. III D, including the triplet states, i.e., the initial wavepacket was centered at either the C_2v or C_s structures. In each case, two simulations were performed, including either all modes or only the key 12 modes: ν₁, ν₂, ν₃, ν₄, ν₅, ν₉, ν₁₆, ν₁₇, ν₁₉, ν₂₀, ν₂₁, and ν₂₅. The tree structures and basis set sizes for converged calculations are given in the supplementary material.

Diabatic state populations for the first 100 fs after a vertical excitation to the S₂ state are shown in Fig. 4 as there is little change during the subsequent period up to 200 fs. Starting from the planar C_2v structure and including all modes [Fig. 4(c)], a fast step-wise decay of population is seen to S₁ and 90% of the population is lost from S₂. Only a small population rise is seen in S₀, with ∼5% having reached the ground-state after 100 fs. A similar picture is seen starting from the C_s minimum energy structure [Fig. 4(d)]. However, there is more population transfer to the ground-state, with ∼15% transfer to S₀ by 100 fs. Simulations including only the key 12 modes [Figs. 4(a) and 4(b)] show very similar population dynamics, indicating that these are the key modes representing the short-time dynamics. Interestingly, the decay to the ground-state is slightly less in the full, 27D calculations. These populations are similar to previous MCTDH dynamics on a simpler vibronic coupling model.²⁹ The system retaining a significant population in S₁ on this time-scale is supported by the fact that cyclobutanone is known to fluoresce.⁵² 

The minimum energy conical intersections (MECIs) between the S₂/S₁ and S₁/S₀ states from the vibronic coupling model were obtained by first minimizing the energy gap between the states of interest and then minimizing the energy along the intersection seam using a Lagrange constraint. The geometries in terms of the normal modes are given in the supplementary material. The S₂/S₁ intersection was found to be at 5.4 units of displacement away from the FC point. Large displacements are made along the ν₅, ν₁₆, ν₁₉, ν₂₀, and ν₂₁ modes. The energy of the intersection is at 5.97 eV and is below the FC point, at 6.32 eV. In a recently introduced classification of conical intersection types,⁵³ this is a direct sloped intersection, which provides the step-wise fast transfer observed. The S₁/S₀ intersection is further away from the FC point, at 11.5 units of displacement at an energy of 4.58 eV, and is thus energetically accessible. However, the displacements required to reach it involve large-scale motions along ν₁ and all the B₂ vibrations; hence, this intersection is not directly accessed and the transfer to the ground-state is limited.

Absorption spectra calculated from the Fourier transform of the autocorrelation function are shown in Fig. 5. Similar results are obtained with the full and 12D simulations. In order to remove artifacts due to the finite propagation time, T, the autocorrelation function was multiplied by an exponential damping of 150 fs. The spectra are shifted by the ground-state zero-point energy. Spectra for absorption to S₁ were also calculated from simulations starting with a vertical excitation to S₁.

The experimental spectrum shows a broad, featureless band from 3.8 eV to 5.0 eV, with the maximum at around 4.2 eV,²⁰ representing the S₁ band. This is a dark state with low intensity. The calculated spectrum from the C_2v structure is found to be in the correct energy region. However, the maximum of this band is lower in energy; by 0.4 eV, it is too narrow and too structured. Starting from the C_s structure, the fine structure of the band is much reduced, although the band is too broad and high in energy.

The experimental S₂ band is also broad, but with a weak progression of peaks between at 6.0 eV and 6.8 eV, with a maximum at 6.5 eV.²⁰ The calculated spectrum from the C_2v structure lies in the correct energy region (around 0.2 eV lower) and has too much structure. The calculated spectrum from the C_s structure is, again, too broad and high in energy.

These results indicate that the spectra are better reproduced by excitation from the C_2v structure despite it being a transition state. This is unsurprising due to the low transition barrier, meaning that at room temperature, the ground state wavefunction will not be purely in a ground-state C_s minimum. It will also access the C_2v structure.

Finally, the expectation values of the normal modes can be examined to see which modes are excited. These are shown in Sec. S4 of the supplementary material. Starting from the C_2v structure, the modes ν₅, ν₁₆, ν₁₉, ν₂₀, and ν₂₁ all show large oscillations, particularly ν₂₁ and ν₅. All are totally symmetric vibrations, and there is no significant loss of symmetry, or out-of-plane motion, over the 200 fs, although mode ν₂ is starting to be activated after 150 fs. Mode ν₂₁, corresponding to the C–O stretch, is seen to extend and oscillates around a positive displacement. Starting from the non-planar C_s structure, the expectation values for the non-symmetric modes ν₁ and ν₂ now also undergo significant dynamics. Mode ν₁ starts with a value around 2 units, and this decreases, showing that the ring becomes more planar over time. ν₂, moves to a negative value, showing the oxygen moves further out of plane.

The final step to give a more complete molecular picture of the dynamics, and obtain the GUED signal, was to perform direct dynamics simulations using the DD-vMCG algorithm. This calculates the potential surfaces on the fly and provides a better description of any long range motion than can be provided by the vibronic coupling model used in the MCTDH calculations of Sec. III E. The surfaces were calculated only at the CASSCF level in order to save time. While this will provide less accurate energetics than the CASPT2 calculations used in the model Hamiltonian, the topography of the surfaces and couplings should be good enough to capture the major atomic motions.

The coordinates used for the simulation are the mass-frequency scaled normal modes of the C_2v transition state. Points were added to the database keeping the C_2v symmetry of the surface by generating appropriate symmetry replicas of each point. An initial calculation of 100 fs was also run starting with the C_2v geometry, and 8 GWPs, so as to ensure the symmetry of the initial points in the database.

The sample in the experiments will be cooled to below room temperature, and hence, the main simulations were started from the C_s minimum. The wavepacket was divided into two partitions with the key 12 modes in one partition and the remaining modes in the second. The database was grown by running an initial simulation for 500 fs, with 16 GWPs in the first partition and 8 in the second. A second simulation collecting further points was then run for 500 fs with 32/16 GWPs in the two partitions. Points were added whenever a structure had an atom over 0.2 bohr away from any structure already in the database. The final database had 22 391 structures. A production run was then carried out for 500 fs, again with a 32/16 GWP basis set, utilizing only the points existing in the database, i.e., not collecting new points. The population dynamics of this final simulation is similar to the model Hamiltonian results shown in Fig. 4. However, less crossing out of the S₂ state occurs, resulting in 46% of the population in S₁, 42% in S₂, and 12% in S₀.

The MECIs between the S₂/S₁ and S₁/S₀ surfaces, provided by the direct dynamics quantum chemistry database, were obtained using an optimization procedure in the same way as those in the vibronic coupling model. The MECI geometries, in both normal mode and Cartesian coordinates, are given in the supplementary material. In both cases, the MECI is slightly closer to the FC point, with the S₂/S₁ MECI now 3.2 units of displacement away and the S₁/S₀ MECI 8.3 units away. The energies have also changed, with the S₂ MECI moving down to 5.2 eV and the S₁/S₀ MECI up to 4.90 eV.

The simulation provides 512 trajectories (the centers of the GWP configurations). Using the procedure in Sec. II E, these can then be used to simulate a PDF, as would be produced by a GUED experiment. The modified scattering intensity, the simulation of the original signal, is shown in the supplementary material. The difference PDF as a function of time is shown in Fig. 6(b), with the reference PDF in Fig. 6(a). The reference is the pair distribution function of the initial, C_s minimum energy structure. The peak at 1.5 Å is due to all of the nearest neighbors. The peak at 2.3 Å is due to the distances C¹–C⁴, C²–C³ (both across the ring) as well as C²–O and C³–O. The peak at 3.2 Å is due to C⁴–O, and the weak peak at 4.3 Å is due to the distance from the protons attached to C⁴ to O. Atom numbering in this section follows the numbering in Fig. 2.

In Fig. 6, loss of intensity is observed in the neighboring atom peak at 1.5 Å with a signal gain around 1.8 Å This is due to the bonds vibrating. The major feature, however, is the strong loss of the 2.3 Å peak and gains at 3 Å and 3.6 Å occurring with a periodicity of 65 fs. The question arises as to what the two new distances are. Plots of the four distances that make up the 2.3 Å peak as a function of time, averaged over the trajectories using the GGP weighting, are shown in Fig. 7. The across ring distance C²–C³ is seen to undergo a small amplitude vibration. The across ring distance C¹–C⁴, however, undergoes a large amplitude vibration, stretching to 2.5 Å, which, along with the C²–O and C³–O distances that extend out to 3 Å, is responsible for the gain just below 3 Å. The gain at 3.5 Å is due to the C⁴–O vibration and the loss of the original signal masked by the gain at 2.5 Å.

The overall picture, hence, is that the molecule stretches along the C–C–O axis. Examination of the bond lengths shows clearly that the C²–C⁴ distance changes minimally, and the stretch is dominated by the C–O group moving away from the other carbon atoms. The dynamics divides into two regimes. In the first 300 fs, the vibrations are coherent and stay in a tight grouping, whereas at later times, these distances spread, and in most, the vibrational motion dies away. An exception is the C–O bond, which appears to be increasing in amplitude. A final note on the dynamics is that there is little change in the ring puckering mode, with an average angle ranging from 166° to 173°. The oxygen out-of-plane motion is a little stronger, with the out-of-plane angle ranging from 12° to −7°. In contrast to earlier AIMS simulations starting from S₁²⁶ and static calculations of the S₁ surface,²⁵ there are no signs of ring opening over the 500 fs.

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 following the initiation of the dynamics by excitation to the S₂ electronic state. Cyclobutanone is a non-trivial molecule for simulations. It has recently been shown that the population dynamics for the related cyclopropanone molecule is highly sensitive to the electronic structure method used.⁵⁴ The ground-state equilibrium structure is a shallow well, and the molecule easily converts between puckered structures through a planar transition state. Consequently, the choice of 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 the triplet states are not significantly populated over the first 500 fs. However, it is also found that while the initial relaxation from S₂ is fast, taking place in around 50 fs, the subsequent relaxation from S₁ is slow. Thus, at later time-scales, there will be crossing to the triplet manifold, and this can result in photo-fragmentation.⁵⁵ 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 increasing in energy toward the end of this time, and it is to be expected that the molecule will fragment to C₃H₆ + C–O at later times. No ring puckering dynamics is observed. The pair distribution function was obtained, which would be measured by the electron scattering. It is seen that this plot is not easy to interpret without the simulations. The simulations give a molecular picture to the pair distances as each peak is due to multiple pairs, as well as overlap. This function, however, clearly reflects the dynamics seen in the simulation.

The simulations are of course not definitive. While a nuclear wavepacket, including quantum effects, was used, the vMCG basis set was not very large, and so this lack of convergence may mean some features are not seen. The potential surfaces used in the final simulations were also only at the CASSCF level, which may not be sufficient for a good description away from equilibrium geometries. The approximations used in the dynamics, e.g., use of the LHA for integrals, and the propagation diabatization procedure may also bring errors. While it is not presently feasible, future work could remove these uncertainties by making it possible to calculate the signal using the full accuracy of the MCTDH method run on the flexible surfaces provided by the DD-vMCG simulations. The integrals required for the GUED signal are not trivial in a grid-based method but are, in principle, possible. In addition, finally, it should be noted that the independent atom model used in calculating the GUED signal may not fully account for changes in the electron density brought about by the non-adiabatic dynamics that couple the electronic motion to the nuclei. A comparison with the experimental result will be a great test of the present capabilities and will show what is still in need of improvement.

The supplementary material includes the coordinates of the ground state structures, vibrational frequencies, cuts through the vibronic coupling model potentials, basis sets used in the ML-MCTDH calculations, information on conical intersections in the singlet manifold, and a plot of the difference scattering spectrum that would be obtained from a GUED experiment. Details of available datasets containing the files from the simulations presented are also given.

This work was funded by the EPSRC under the program Grant COSMOS (No. EP/X026973/1). O.B. also acknowledges UCL for funding. A.F. acknowledges financial support from the Cluster of Excellence “CUI: Advanced Imaging of Matter” of the Deutsche Forschungsgemeinschaft (DFG)—EXC 2056—Project No. 390715994, from the International Max Planck Graduate School for Ultrafast imaging & Structural Dynamics (IMPRS-UFAST) and from the Christiane Nüsslein-Vollhard Foundation.

The authors have no conflicts to disclose.

Olivia Bennett: Conceptualization (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Antonia Freibert: Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). K. Eryn Spinlove: Conceptualization (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Graham A. Worth: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article and its supplementary material.