Using a multistate mapping approach to surface hopping to predict the ultrafast electron diffraction signal of gas-phase cyclobutanone Open Access

The photochemistry of cyclic ketones has been studied extensively.^1–6 Despite their small size, they have a rich photochemistry with multiple competing pathways that include dissociation, fluorescence, and intersystem crossing.^7,8 The relative importance of the different pathways is closely linked to the size of the organic ring and the associated ring strain.⁷ The photodissociation of cyclic ketones has been of particular interest ever since the pioneering work of Norrish.⁹ 

In this class of molecules, a target of recurring interest is cyclobutanone, (CH₂)₃CO.^{1,2,7,10–12} Early experimental studies identified the photoproducts of cyclobutanone as propylene (cyclopropane) and carbon monoxide (C₃H₆ + CO) or ethylene and ketene (CH₂CH₂ + CH₂CO). The corresponding quantum yields were found to be 40% and 60%,¹ although these were later shown to have a strong dependence on the excitation wavelength.⁸ Recent experimental work includes ultrafast transient absorption spectroscopy on cyclobutanone in solution.⁷ Upon excitation to the S₁ (nπ*) state by UV pulses in the range 255–312 nm, Kao et al. found that singlet dissociation pathways are dominant at shorter wavelengths, predominantly via Norrish Type-1 α cleavage, with minor contributions from ketene formation. At higher excitation energies, such as 200 nm, the n-3s Rydberg state comes into play.^13–15 In this regime, using time-resolved mass spectrometry and photoelectron spectroscopy, Kuhlman et al. found that a ring puckering mode effectively couples the S₂ and S₁ states, allowing for rapid internal conversion.¹⁶ 

The experimental studies have been accompanied by theoretical work. Electronic structure calculations using complete active space self-consistent field, CASSCF(10,8), and multistate complete active space second-order perturbation theory, MS-CASPT2(10,8), predicted a small barrier for ring-opening on the S₁ state, with the possibility of minor contributions from barrierless dissociation in the T₁ state.¹⁷ Building on this, simulations of cyclobutanone using the ab initio multiple spawning (AIMS) method found a ring-opening mechanism to be dominant with a small portion of trajectories yielding CH₂=CH₂ + CH₂CO and CH₂CH₂ + CH₂=CO fragmentation.¹⁸ We also note that Kuhlman et al. constructed a five-dimensional linear-vibronic Hamiltonian (LVH) model to simulate the S₂/S₁ decay in cyclobutanone using the multiconfigurational time-dependent Hartree (MCTDH) method, which succeeded in replicating the short time constants of 0.95 ps.¹⁹ 

This paper is motivated by the prediction challenge based on an ultrafast electron diffraction (UED) experiment carried out at the SLAC MeV-UED facility.²⁰ UED has emerged as a powerful experimental technique for observing molecular dynamics in the recent decade,^21–25 alongside ultrafast x-ray scattering.^26–29 Notably, both techniques have made advances in extracting information that extends beyond structural dynamics.^30,31 We use this challenge as an opportunity to benchmark the electronic structure calculations for cyclobutanone. This comparison includes non-standard selected-configuration interaction (SCI) electronic structure methods that overcome the issue of large active space selection, constitute a black-box alternative to complete active space (CAS) methods, and provide complementary benchmarks.^32,33 The nonadiabatic simulations are carried out using the newly developed multistate mapping approach to surface hopping (multistate MASH) method.³⁴ Furthermore, we use the simulations to predict the UED signals and analyze these to determine the main contributions to the scattering.

Many computational methods have been developed to simulate nonadiabatic dynamics over the past decades. Some, such as MCTDH^35,36 and multi-layer MCTDH,³⁷ can provide numerically exact results. However, these methods require high-accuracy precomputed potential energy surfaces (PESs), which, combined with the exponential scaling of quantum mechanics, imposes severe limitations on the number of degrees of freedom that can be treated. To circumvent this problem, the nuclear wavefunction can be represented as a linear combination of traveling Gaussian basis functions or classically guided trajectories, and the required electronic structure quantities can be calculated on-the-fly, also known as direct dynamics. Such approaches can significantly reduce the phase space explored and, therefore, the computational cost per degree of freedom. On-the-fly methods include examples where the nuclei are propagated using equations of motion derived from the variational principle, notably variational multiconfigurational Gaussian (vMCG),^38–40 or semiclassically, which includes methods such as full or ab initio multiple spawning (FMS/AIMS)^41–44 and multiconfigurational Ehrenfest (MCE),^45–47 or classically, such as trajectory surface hopping (TSH).⁴⁸ In addition to the above-mentioned methods, mapping methods treat the electronic and nuclear degrees of freedom in an equal manner. Examples include the Meyer–Miller–Stock–Thoss (MMST) mapping Hamiltonian, symmetrical quasi-classical windowing, and generalized spin mapping.^49–52 Ideas from this last category informed the development of the surface hopping variant used in this paper (see Sec. II A 2 below). For an in-depth overview of different methods for nonadiabatic dynamics, the reader is directed toward the edited volume in Ref. 53.

For nonadiabatic molecular dynamics, TSH is perhaps the most commonly employed approach to support the interpretation of experiments.^54–61 In TSH, an ensemble of individual trajectories is used to represent the propagation of the nuclear wavepacket. In each trajectory, the nuclei are treated classically and the electronic amplitudes are used to predict hops to other electronic states. For a full discussion of TSH, please see Ref. 62, while we continue by highlighting the key aspects of the fewest switches surface hopping (FSSH) variant of TSH.

The algorithm presented above is robust but also causes one of the most well-known shortcomings of TSH, the so-called overcoherence issue. A number of, more or less ad hoc, decoherence correction schemes have been developed in response.^67,68 The method presented next aims to alleviate the need for such ad hoc corrections.

To start trajectories, multistate MASH determines the nuclear initial conditions in the same manner as TSH, mapping the ground-state wavefunction onto a classical phase space. However, as in other mapping approaches, MASH also considers the electronic coefficients c as a phase space variable. Therefore, rather than describing the initial state as pure, such as done in FSSH, multistate MASH represents the initial state as a distribution over c = (c₁, …, c_N). In practice, this is done by random assignment of the complex-valued coefficients, c_i = x_i + iy_i, for each electronic state from Gaussian distributions with zero mean and standard deviation one for both the real and imaginary components. The resulting set of coefficients is normalized such that $1=∑i=1N|ci|2$⁠. Coefficients for which the initial state of interest does not have the largest absolute value are resampled until suitable sets are found. It is worth noting that other choices of initial distributions are possible, but in most cases, have been found to yield similar results.³⁴ Therefore, we use this method of selecting initial conditions due to its ease of implementation.

In the two-state case, Mannouch and Richardson have shown that this prescription overestimates the rate of transfer compared to the Landau–Zener theory. It does, however, maintain a physical (non-negative) population estimate with a low number of trajectories. Note that this approach differs from what was previously used in both MASH and multistate MASH. Further information on alternative population estimators is given in the supplementary material.

One potential side-effect of multistate MASH is hopping between states that are not coupled. It is possible to nevertheless accept all hops, regardless of the coupling between the states, which was the approach taken in Ref. 34. However, due to the risk of unphysical behavior, we introduce an energy criterion for all hops, such that |E_a − E_b| < 0.055 Hartree (≈1.5 eV) to prevent erroneous hops between uncoupled states. A value of ≈1.5 eV was selected based on previous dynamics with TSH; however, further testing of this value is still required. Hops that do not meet this criterion are treated on the same footing as reflected hops, with the total velocity of all atoms reflected. Note that this direction of rescaling differs from the previous work with MASH and is chosen due to simplicity.

We note that further improvements in the scattering signal are straightforward. For instance, the alignment effects resulting from the linear polarization of the pump laser can be accounted for Ref. 76 and corrections to the form factors due to relativistic effects can be made;⁷⁷ however, for the latter, we note that this effect is comparatively minor for the electron energies at the SLAC MeV-UED source. It is also possible to calculate the scattering cross sections from ab initio electronic wavefunctions, which makes it possible to accurately account for the effects on scattering from chemical bonding, electron correlation, and change in the electronic state due to excitation by the pump pulse.^78–80 

A ground-state minimum energy structure of cyclobutanone was obtained using state-averaged CASSCF,⁸¹ SA(3)-CASSCF(12,12)/aug-cc-pVDZ, i.e., state-averaging over the three lowest singlet states, an active space with 12 electrons in 12 orbitals with the aug-cc-pVDZ Cartesian basis set.⁸² Further details on the active space are given in the supplementary material. All electronic structure calculations were performed in OpenMolcas version v23.02.^83,84 The ground-state minimum energy structure was subjected to a frequency calculation yielding no imaginary frequencies, confirming the structure to be a true minimum. Further excited state structures were optimized using SA(3)-CASSCF(12,12) to obtain the S₂ minimum, two S₂/S₁ minimum energy conical intersections (MECIs), and two S₁/S₀ MECIs. The selected structures were subjected to linear interpolation in internal coordinates (LIIC) to yield a progression of intermediate structures. Potential energy curves along the LIICs were calculated using SA(3)-CASSCF(12,12) and extended multistate complete active space second-order perturbation theory, XMS-CASPT2, with the same (12,12) active space as for the CASSCF calculations (see the supplementary material).⁸⁵ For stability, an imaginary shift of 0.5 a.u. was applied for all XMS-CASPT2 calculations.⁸⁶ 

For further comparison, the potential energy curves were also calculated using a modified version of adaptive configuration interaction (ACI) and Monte Carlo configuration interaction (MCCI).^87–89 The state-averaged MCCI calculations were performed using MCCI V4^89,90 with the SA error-controlled calculations performed using a locally modified development version of GeneralSCI 1.0.⁹¹ Further details are given in the supplementary material.

The calculated ground-state frequencies were used to generate a Wigner distribution from which 1000 geometries were sampled. At each geometry, excitation energies and oscillator strengths were calculated using SA(3)-CASSCF(12,12). The photoabsorption cross section was constructed using the nuclear ensemble approach (NEA), convoluted by Lorentzian functions with a phenomenological broadening parameter of 0.05 eV.⁹² Initial conditions for the multistate MASH trajectories were sampled from this Wigner distribution using an implicit laser pulse of width 0.2 eV, centered at 6.2 eV. All trajectories were initiated from the bright S₂ (n-3s Rydberg) state using multistate MASH (see Sec. II A 2 for more details)³⁴ implemented in a developmental version of SHARC.⁹³ The electronic structure calculations were carried out using SA(3)-CASSCF(12,12)/aug-cc-pVDZ with all three electronic states available in the dynamics. The trajectories were propagated using a 0.5 fs time step, with 25 substeps, for 300 fs using a local diabatization scheme. Previous applications of multistate MASH have employed time steps in the region of 0.25 fs, and due to the similarities of the algorithm to TSH, similar time steps are expected to be adequate.³⁴ For trajectories in which the total energy changes by more than 0.25 eV, a value selected based on previous studies,^57,94–96 the propagation is stopped, and only time points with the correct total energy are used in the analysis. In addition, we note that a subset of trajectories becomes “trapped.” This occurs when the trajectories enter a region where a reflected hop occurs but the reflected hop fails to remove the trajectory from the “forbidden” region. Such trapped trajectories are only included in the analysis until the point that they get trapped. For a summary of the trajectories removed from the analysis, see the supplementary material. Populations of the ensemble of trajectories are estimated according to Sec. II A 2. Finally, rotationally averaged total scattering patterns were calculated using the geometries from the multistate MASH trajectories and the independent atom model (see Sec. II B) for the initial 200 fs of the dynamics.

The character of the two lowest-lying S₁ and S₂ electronic states of cyclobutanone has previously been identified as nπ* and n-3s, respectively.^13,16 The next set of excited states are the higher-lying 3p Rydberg states, which appear a full 0.7 eV above the n-3s Rydberg state.¹⁴ Triplet states have been shown to be important for the photodynamics upon excitation with long wavelengths.⁷ However, due to the short excited state lifetime of cyclobutanone following excitation into the S₂ state, triplet states are less likely to make a significant contribution to the current photodynamics. For these reasons, we have chosen to focus our attention on an accurate description of the two lowest-lying singlet states S₁ and S₂.

A summary of the low-lying excited singlet states available to cyclobutanone at the equilibrium geometry is presented in Table I. Both SA(3)-CASSCF(12,12) and XMS-CASPT2(12,12) yield an S₁ state with nπ* character at the Franck–Condon (FC) geometry and an S₂ state with the n-3s Rydberg character. Furthermore, the excitation energies for SA(3)-CASSCF(12,12) are in excellent agreement with XMS-CASPT2(12,12), with SA(3)-CASSCF(12,12) energies showing a constant difference of 0.04 eV for both the excited states. Excitation energies for SA(3)-ACI and SA(3)-MCCI are also presented in Table I. These methods give a good description of the low-lying singlet states of cyclobutanone with excitation energies deviating by ∼0.05 eV for the valence nπ* state and ∼0.1 eV for the n-3s Rydberg state compared to the XMS-CASPT2(12,12) benchmark. In the ACI simulation, a large orbital window, corresponding to CAS-CI(12,36), has also been tested. The large orbital window gives results analogous to the ones obtained with CASSCF(12,12), providing further support for this choice of active space.

A theoretical photoabsorption cross section of cyclobutanone calculated using SA(3)-CASSCF(12,12) is shown in Fig. 1, allowing for a direct comparison with an experimental absorption spectrum. Good agreement can be observed between the two spectra, further validating our electronic structure method along with confirming the assignment of states using the insets shown in Fig. 1.

After the photophysics of cyclobutanone has been studied in the FC region, the different deactivation pathways available upon photoexcitation can be characterized. A series of potentially important optimized geometries for these pathways are shown in Fig. 2, optimized using SA(3)-CASSCF(12,12)/aug-cc-pVDZ. The S₀ minimum energy geometry (S₀, upper left) can be characterized by a single carbon atom breaking the planar symmetry of the molecule, along with a C=O bond length of 1.19 Å and an α–CC bond of 1.58 Å. When the molecule is excited by a 200 nm pulse, the S₂ excited electronic state is populated at the FC region and can potentially relax into the S₂ minimum. The S₂ minimum energy geometry (S₂, upper center) displays subtle changes in bond lengths, along with a planarization of the carbon ring. Two S₂/S₁ MECIs were successfully optimized, with the S₂/S₁ MECI (S₂/S₁, lower left) already reported previously.¹⁹ This MECI displays a rather severe compression of the C=O bond. The second S₂/S₁ MECI (S₂/S₁, upper right) exhibits a dissociated α–CC bond alongside a compressed C=O bond. Two further MECIs are shown in Fig. 2, corresponding to MECIs between the ground and first excited state. Both these S₁/S₀ MECIs have been reported previously by Liu et al.¹⁸ The first of these MECIs (S₁/S₀, lower center), termed CI-3 in Ref. 18, breaks the ring structure to produce CH₂=CH₂ and O=CCH₂. The second S₁/S₀ MECI (S₁/S₀, lower right), termed CI-1 in Ref. 18, is closely related to the second S₂/S₁ MECI (S₂/S₁, upper right), with a dissociated α–CC bond. However, in contrast to the S₂/S₁ MECI, the S₁/S₀ MECI contains a slightly stretched C=O bond relative to the S₀ geometry. It is obvious from this static perspective of excited state geometries, that cyclobutanone possesses a rich photochemistry. However, the relative importance of each of these structures is still unknown, as are the final products of each conical intersection (CI).

To gain a better understanding of the PESs of cyclobutanone and to benchmark SA(3)-CASSCF(12,12) away from the FC region, an LIIC was performed on some of the key geometries shown in Fig. 2. Four structures were selected: the S₀ minimum (Fig. 2, S₀, upper left), the S₂ minimum (Fig. 2, S₂, upper center), the dissociative S₂/S₁ MECI (Fig. 2, S₂/S₁, upper right), and finally, the α–CC bond dissociation S₁/S₀ MECI (Fig. 2, S₁/S₀, lower right). Here, the two CIs were selected due to being the most commonly encountered CIs during the dynamics presented in Sec. IV C.

Upon photoexcitation to the S₂ state, cyclobutanone in the S₀ minimum structure can relax to the closeby S₂ minimum via planarization of the carbon ring (the green line shown in panel 1 of Fig. 3). Here, cyclobutanone undergoes a dissociation of an α–CC bond, resulting in a barrier height of 0.21 eV for SA(3)-CASSCF(12,12), somewhat below the 0.49 eV predicted by XMS-CASPT(12,12), with the barrier heights given as energy differences relative to the S₂ minimum. Note that the absolute values for transition state barriers in LIICs should be treated with care as they are known to be overestimated. Regardless, as a result of this smaller reaction barrier, SA(3)-CASSCF(12,12) is likely to estimate a faster rate of decay from S₂ to S₁. A steep drop in energy is observed for both the methods after the transition barrier leading down to the S₂/S₁ MECI (see Fig. 2, upper right and inset of Fig. 3). This S₂/S₁ MECI is peaked and has a single pathway according to the criteria set out in Ref. 98. This allows for efficient funneling onto the S₁ state (the blue line shown in Fig. 3) where a barrier-less decay is observed to the structurally similar S₁/S₀ MECI (see Fig. 2) and where the S₁/S₀ MECI is observed with a peaked bifurcating branching space.

It can be deduced from Fig. 3 that there is excellent agreement between SA(3)-CASSCF(12,12) and XMS-CASPT2(12,12), with slight differences observed for the approximate transition barrier height. Nevertheless, overall the SA(3)-CASSCF(12,12) results offer an excellent comprise between speed and accuracy for the on-the-fly nonadiabatic dynamics of gas-phase cyclobutanone. Additional LIICs for the same geometries are given in the supplementary material (Figs. 3 and 4) for both SA(3)-ACI and SA(3)-MCCI. These further support our assessment. We also note that this indicates that good agreement can be obtained with methods that potentially offer a more black-box approach to the calculation of multireference wavefunctions than the CAS family of methods.

A total of 229 trajectories were initiated on the S₂ state. Trajectories were removed from the statistics if either a change of 0.25 eV or more was observed in the total energy or if consecutive reflected hops indicated that the trajectory was trapped in a “forbidden” region of the PES. This resulted in a significant number of discarded trajectories. For a full summary of discarded trajectories, see the supplementary material (Fig. 2). To briefly summarize, after 50 fs, the number of trajectories had dropped to 146, and by 200 fs, only 18 trajectories of the initial 229 trajectories had survived, while at 300 fs, this was down to 10 trajectories. At 300 fs, a total of 71 trajectories had been discarded due to consecutive reflected hops, while 148 trajectories had been discarded due to issues with the energy conservation. We acknowledge that this could potentially cause artifacts in the simulations, especially if certain processes are systematically removed or enhanced by the removal of trajectories.

The adiabatic populations for the three electronic states S_0–2 used in our simulations are shown in Fig. 4 for the initial 200 fs of the dynamics. A rapid decay of the S₂ state (green line) can be observed after ∼8 fs. This decay is correlated with a population increase in S₁ (teal line), indicating that the fast initial dynamics on S₂ rapidly reaches an S₂/S₁ CI. For the early time dynamics (t < 30 fs), the ground state can be seen to be inaccessible from the S₂ state (panel 1 shown in Fig. 3). After 30 fs, there is a steady increase in the ground-state S₀ (purple line) population, showing passage through the S₁/S₀ CI and, therefore, formation of the photoproducts linked to each CI (Fig. 2). After 100 fs, we observe a decreased rate of population transfer from S₂ to S₁, resulting in depletion of S₁ and an overall decrease in the population of the two excited states matched by a continued increase in the population of the ground state. After $∼175$ fs, the population transfer is largely over, although errors in this region will be high due to the small number of trajectories. The populations shown in Fig. 4 indicate ultrafast non-radiative deactivation pathways for cyclobutanone, there is the likely possibility of further non-equilibrium ground-state chemistry as a result of very hot molecules in the ground state following the decay dynamics.

As shown previously, the deactivation pathways available to cyclobutanone imply the cleavage of one or two bonds in the ring structure. To identify these different pathways, the bond lengths for all multistate MASH trajectories are plotted, separated into α- and β–CC bonds according to the inset shown in Fig. 5. Four distinct outcomes are observed, as follows:

• There is no bond breaking (purple), and the molecule remains in the FC region. These are characterized by short α and β–CC bond lengths (R < 2.2 Å). It is important to note that the trajectories included in this group do not progress over 90 fs.

• There exists a single α–CC bond breaking (teal). This pathway gives rise to a ring-opened structure, i.e., CH₂CH₂CH₂CO. This is indicated by an elongation of an α–CC bond, but due to the carbon backbone still being intact, the distance between the two carbon atoms is limited to R ∼ 4.5 Å. This pathway is a result of passing through the S₁/S₀ CI (lower right, Fig. 2).

• Two α–CC bonds break to liberate CO, where there is a large increase in α–CC bonds to values above 4.5 Å (blue). Attempts to optimize an MECI using SA(3)-CASSCF(12,12) involving CO dissociation were unsuccessful; however, trajectories were observed with CO dissociation after having decayed to the ground state via the α–CC dissociation S₁/S₀ CI (lower right, Fig. 2). This means the CO dissociation takes place on the ground state rapidly after passing through a CI. It is worth noting that for option (iii), a single trajectory also displays ground state dissociation of CH₂CH₂CH₂ into CH₂CH₂ and CH₂, which is shown in Fig. 5 as the blue line displaying increases in both α and β–CC bonds.

• An α–CC bond and β–CC bond breaking (pale green) showing the stepwise ring-breaking mechanism where ethene and ketene are produced. Here, an α–CC bond breaks first, followed by a β–CC bond breaking. This can be linked to the S₁/S₀ MECI shown in Fig. 2 (S₁/S₀, lower center).

UED cross section percent differences shown in Fig. 6 present four distinct features at s ≈ 1.25, 2.4, 4.6, and 7.4 Å⁻¹. These features constitute a fingerprint of the structural dynamics taking place in cyclobutanone; therefore, they can be interpreted with the help of additional static signals for the key geometries discussed in Sec. IV B (purple line shown in Fig. 7). A static signal for a fixed time [80 fs for panels (a) and (b) and 150 fs for panels (c) and (f)] is also shown in Fig. 7 for each pathway (green lines), showing some broadening of the peaks due to vibrational dispersion, along with subtle shifts of peak maxima. In addition, the electron scattering cross sections are calculated for the four individual outcomes using groups of trajectories in the dynamics simulation [(i), (ii), (iii), (iv) in Sec. IV C] and shown in Fig. 8 to aid the interpretation of the results. Due to significant similarities between the photoproducts and the high degree of symmetry in cyclobutanone, there is a large number of overlapping features in all static and individual cluster signals (see Figs. 7 and 8). Based on these similarities, we postulate that the strong negative signal at ∼1.25 Å⁻¹ and the strong positive feature at ∼2.4 Å⁻¹, common to all photoproducts is the breaking of one or two α-bonds to produce a ring-opened structure, CO and CH₂CH₂ CH₂, or ethene and ketene. These two features, which reach a maximum at t ∼ 50 fs, are then related to the electronic state population transferring through the S₂/S₁ CI and subsequently through the S₁/S₀ CI [see Figs. 7(d), 7(e) and 7(f)]. This timescale is consistent with what is observed in the populations of S₁ and the S₀ ground state, as shown in Fig. 4, but significantly faster than that observed in Ref. 19. It is possible that our simulations underestimate the time in S₂ due to a smaller barrier, as discussed in the context of Fig. 3 for SA(3)-CASSCF(12,12) compared to the matching XMS-CASPT2, resulting in somewhat fast decay and product formation.

The other UED feature common to the total signal (Fig. 6) and all individual pathways (Fig. 8) is an oscillatory positive signal at s ≈ 7.4 Å⁻¹. For the initial 50 fs of simulation, in this s region, one can see strong coherent motions with an oscillatory period of $∼15$ fs, indicating a possible coherent vibrational motion. All the bonds between heavy atoms (C and O) are plotted for an exemplary trajectory in the supplementary material (Fig. 6). Only the CO bond length has a period that matches these oscillations, suggesting that the carbonyl stretch may be responsible for this feature.

An additional feature can be observed at s ∼ 4.6 Å⁻¹ with weak positive intensity for the initial 80 fs of the time-resolved UED becoming weakly negative after 80 fs. This positive weak feature (4%–6%) only shows in those pathways where β-bond breaking does not occur [Figs. 8(a) and 8(c)] and it is negative in the one where this happens [Fig. 8(b)]. When looking at the static signals shown in Fig. 7, one can see that the CO dissociation channel [Fig. 7(f)] yields a positive signal at s ∼ 4.6 Å⁻¹ and that this feature is also present in the ring-open structures [Figs. 7(c) and 7(d)] at s ∼ 4.9 Å⁻¹. Therefore, one can relate the sign and strength of this feature to the separation between the two carbons forming the β bonds in the molecule (see the inset in Fig. 5). The change of sign for this feature at t ≈ 80 fs, together with the great influence path (iv) has on the dynamics, suggests that ethene and ketene formation can be directly observed by looking at this region of the time-resolved UED signal.

It is important to note that the FC group presents oscillatory behavior in all the aforementioned UED features. These oscillations can be explained by taking into account that the molecule is vibrationally excited (i.e., hot) when the FC region is accessed, giving rise to high-amplitude motions. This behavior is not seen in the averaged signal due to the small contribution of this feature to the overall dynamics.

In summary, motivated by an experiment at the SLAC MeV-UED facility and the associated prediction challenge, we have applied multistate MASH and SA(3)-CASSCF(12,12)/aug-cc-pVDZ to identify the gas-phase photodynamics of cyclobutanone after photoexcitation by 200 nm pump pulses into the n-3s Rydberg state. Mechanistic details have been identified using both static analysis and nonadiabatic dynamics, finding a fast decay from the S₂ state into the S₁ via two CIs. A further two CIs allow for a similarly fast access to the ground state. This results in the same set of photoproducts previously observed, namely, α ring-opening and ethene + ketene production, the former of which can further dissociate on the ground electronic state to liberate CO.

For direct comparison with the experimental results, we have also predicted the gas-phase isotropic UED signal of cyclobutanone (see Fig. 6). We find that there is significant overlap between several of the products, meaning that the UED signal contains multiple features. However, several distinct features correspond to the formation of reaction products and loss of cyclobutanone in the equilibrium geometry. The timescale of the reaction can also be inferred from the UED signals, which, unsurprisingly, closely matches that of our simulated reaction dynamics. However, we do note that the reaction, in our simulations, proceeds significantly faster than previous observations.^16,19

Significant caveats remain for our predictions. The instability of the active space means that we could potentially be blind to certain reaction channels. In addition, the value obtained for each channel’s quantum yield is susceptible to significant errors from the choice of the electronic structure method and the associated instabilities. Furthermore, the loss of trajectories means that the predictions for later times suffer from poor statistics and should be treated with caution. Notably, the convergence of properties such as quantum yields can be quite slow with respect to the number of trajectories.⁹⁹ Nevertheless, experimental signals, especially for scattering, can converge for comparatively small numbers of representative trajectories.^26,100–103 This, combined with our reasonably accurate electronic structure (known to be important for the validity of the trajectories^96,99,104) lends credence to our predictions. In a broader context, this highlights the need for robust “black-box” electronic structure methods compatible with on-the-fly nonadiabatic dynamics simulations. Methods such as ADC(2) or TDDFT already exist but are not universally applicable as they have well-documented flaws. A potential way forward could be the selected CI family of methods,^89,90 some benchmarking of which has been already been included in our paper. Further development and benchmarking of these methods is required to confirm whether they indeed offer an accurate electronic structure in a more robust, black-box fashion. In addition, it is important to add that the truncation of trajectories due to uncoupled hops affects the statistics of the dynamics and, therefore, needs to be refined in future work. Finally, the multistate MASH used in the simulations is still being developed, and we anticipate that refinement of the algorithm will improve the results further.

The supplementary material describes the active space employed in our calculations and presents the LIICs for the ACI and MCCI calculations along with further computational details for these methods. It also includes the total isotropic UED signal calculated with a tighter total energy threshold together with a description of the trajectories discarded due to the energy threshold. Finally, a summary of population estimators and a CO vibrational analysis is given.

The authors thank the organizers of the cyclobutanone prediction challenge for arranging this special edition. A.K., M.J.P., L.H., and A.P. acknowledge funding from the Leverhulme Trust (Grant No. RPG-2020-208), and A.K., A.M.C., and M.J.P. further acknowledge EPSRC EP/V006819 and EP/T021675/1. A.K. also acknowledges EPSRC EP/V049240 (together with M.S.), COSMOS Programme Grant (No. EP/X026973/1), and Grant No. DE-SC0020276 from the U.S. Department of Energy. J.E.R. was funded by a mobility fellowship from the Swiss National Science Foundation under Award No. P500PN_206641/1.

The authors have no conflicts to disclose.

Lewis Hutton: Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Andrés Moreno Carrascosa: Conceptualization (equal); Data curation (supporting); Formal analysis (equal); Funding acquisition (equal); Writing – original draft (equal); Writing – review & editing (lead). Andrew W. Prentice: Data curation (supporting); Formal analysis (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Mats Simmermacher: Data curation (supporting); Formal analysis (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Johan E. Runeson: Data curation (supporting); Formal analysis (supporting); Writing – review & editing (supporting). Martin J. Paterson: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Writing – review & editing (supporting). Adam Kirrander: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (lead); Writing – review & editing (equal).

All data and codes are available at reasonable request.