Ultrafast electron diffraction of photoexcited gas-phase cyclobutanone predicted by ab initio multiple cloning simulations Open Access

Ultrafast electron diffraction (UED) has evolved into a powerful method for structural dynamics.^1,2 Although UED and the closely related method of ultrafast x-ray scattering^3,4 arguably provide the most direct access to structural dynamics in photoexcited molecules, the interpretation of experiments is nontrivial. Despite significant progress in the development of inverse methods, which aim to produce a (time-dependent) molecular model commensurate with the experimental data,^5–9 the gold standard for interpreting ultrafast experiments remains the comparison to high-quality simulations of the photoexcited target molecule. However, such simulations are challenging, and their veracity depends keenly on numerous methodological choices. As a much needed step toward surveying and evaluating good practice, the Journal of Chemical Physics recently announced the Prediction Challenge: Cyclobutanone Photochemistry, to which this paper is a response.

There have been a number of different methods put forth under the collective term of nonadiabatic dynamic methods. Methods can broadly be partitioned into grid-based and on-the-fly methods. Grid-based methods, such as multiconfigurational time-dependent Hartree (MCTDH) and its multilayer variant (ML-MCTDH), are capable of giving high-accuracy results but at the cost of expensive precomputed potential energy surfaces (PESs).^10,11 To alleviate this need, on-the-fly trajectory-based methods are often employed, ranging from trajectory surface hopping (TSH), where nuclei are treated classically and electronic quantities are calculated from electronic structure methods, nonadiabatic events are then simulated with hopping events between electronic states.¹² To increase the quantum nature of the nuclei, full multiple spawning (FMS) or ab initio multiple spawning (AIMS) can be used. FMS/AIMS uses a series of Gaussian basis functions to recreate a nuclear wavepacket; each of these Gaussians is propagated classically by their centroid position, with spawning events happening between surfaces in areas of configuration space with high nonadiabatic coupling.^13–16 Another technique using classically guided Gaussians is multiconfigurational Ehrenfest (MCE). More details on MCE with cloning, which is known as Ab Initio Multiple Cloning (AIMC), can be found in Sec. II. However, briefly, multiple Gaussians are propagated according to Ehrenfest forces, and a cloning algorithm is employed in a similar way to spawning in AIMS.^17,18 For a more detailed summary of the methods available, see Ref. 19. With such a large number of possible methods, questions naturally arise about the accuracy of these methods when calculating experimental observables.

Cyclobutanone [(CH₂)₃CO] is a small cyclic ketone shown both experimentally and theoretically to have a rich photochemistry.^19–25 After cyclobutanone is excited into its S₁ (nπ*) state, both C₃H₆/CO and C₂H₄/CH₂CO are produced in a 2:3 ratio although this ratio has been shown to be strongly wavelength dependent.²³ Further studies into the photodynamics of cyclobutanone after excitation into S₁ found an α ring-opening reaction to be the dominant deactivation pathway and can result in the formation of ketene (CH₂CO). It should be noted that this study takes place in a solvent environment.²⁶ These reaction pathways have also been observed from static calculations and AIMS simulations.^27,28 In addition, the role of triplet states is expected to be minimal unless the excitation energy is low enough.²⁶ The second absorption band arises from an n-3s Rydberg transition (S₂).^29–31 Studies probing the dynamics initiated from the S₂ state using time-resolved mass spectrometry and time-resolved photoelectron spectroscopy indicate that an out-of-plane ring puckering motion is present, allowing for the rapid S₂/S₁ decay.³² This mechanism was supported by a five-dimensional linear vibronic Hamiltonian used for MCTDH simulations.³³ It has been suggested that after decay to the S₁ state, the deactivation mechanism of cyclobutanone should be ultrafast, producing a product of m/z = 42.³² 

The challenge is motivated by an experiment at the SLAC Megaelectronvolt ultrafast electron diffraction (SLAC MeV-UED) facility, where a gas-phase sample of cyclobutanone is irradiated with a 200 nm laser pulse and time-resolved UED signals are recorded. At this excitation energy, a low-lying n → 3 s (S₂) Rydberg state in cyclobutanone is excited.^29–32 

In this work, the photodynamics of cyclobutanone is simulated using the Ab Initio Multiple Cloning (AIMC)^34–36 approach, which is, in principle, a fully quantum, formally exact methodology based on using Gaussian coherent states propagated by Ehrenfest trajectories as a basis for the quantum dynamics of nuclear wave functions. AIMC was successfully applied before^37–41 to simulate the process of the photodissociation of a number of heterocyclic molecules. Based on AIMC dynamics results, the isotropic gas phase time-resolved UED pattern of cyclobutanone photoexcited using a 200 nm pulse is calculated for the initial 200 fs of dynamics, allowing for direct comparison to experimental data.

The total contribution of two clones into the wave function |Ψ(R, r, t)⟩ is exactly the same as that of the original basis function. However, cloning increases the size of the basis set, creating additional flexibility as two clones can now move in different directions. The cloning is applied when the magnitude of the breaking force $FI(br)=|aI|2∇VI−∑J|aJ|2∇VJ$ exceeds a threshold, and, at the same time, the magnitude of non-adiabatic coupling is below a second threshold. Cloning is an extremely important part of the AIMC method, as it allows AIMC to reproduce the bifurcation of the wave function at conical intersections.

The trajectories in the AIMC approach can be calculated independently using potential energy forces and non-adiabatic coupling vectors calculated “on the fly” by an electronic structure code. Then, the time-depended Schrödinger equation for amplitudes c_n is solved in post-processing on the precalculated trajectory-guided basis (1).

In practice, to achieve good convergence, a number of sampling techniques have to be used. Swarms of coupled trajectory guided Gaussians, as well as their trains guided by the same trajectories, are among those techniques.⁴⁴ It has been demonstrated that MCE can produce results that are well converged,⁴⁵ and AIMC, its ab initio direct dynamics version, is more accurate than surface hopping or Ehrenfest dynamics.^46,47 A technique that allows for taking into account the pulse shape and dynamics that occur during the excitation has been developed.⁴⁸ In its simplest form, which is used in the present paper, AIMC can yield a qualitative or semiquantitative picture of the process, similar to that given by surface hoping, AIMS,¹⁹ and many other popular techniques.

For the modeling of the UED signals, we anticipate that the experimental data in this challenge can be modeled reliably using the independent atom model (IAM). This approximates the scattering signal as a coherent sum of scattering from isolated atoms centered at the positions of the nuclei in the target molecule. Notably, this model excludes the contribution of the bonding electrons and the characteristics of the electronic states of the molecule. Should the quality of the experimental data necessitate that these effects are accounted for, then numerical codes capable of this exist,^49–52 albeit at significantly higher computational costs.

The trajectories were calculated using our own AIMC code, where potential energies, forces, and non-adiabatic couplings were given “on the fly” by the MOLPRO⁶³ electronic structure package at the SA(3)-CASSCF(12,12)/aug-cc-pVDZ level of theory. We note that the electronic structure method has been benchmarked in another paper submitted to the same challenge by one of the co-authors (AK). In brief, three electronic states were taken into consideration: a ground state and two lowest singlet excited states. Higher energy Rydberg states have been shown to exist (3p character) but are unlikely to be important for dynamics after excitation into S₂; therefore, we do not include them in our simulations.³³ We also do not take triplet states into consideration in this work, as they have been shown to only play a role in dynamics upon excitation with long wavelengths.²⁶ The initial positions and momenta for all trajectories are randomly sampled from the ground state vibrational Wigner distribution using vibrational frequencies and normal modes calculated at the same level of theory. This ground state wavepacket is then simply lifted to the second excited state within the Condon approximation. As in our previous simulations,^34,37–41 the cloning thresholds were taken as 5 × 10⁻⁶ a.u. and 2 × 10⁻³ a.u. for the magnitude of breaking force and non-adiabatic coupling, respectively. In order to prevent uncontrolled growth of the number of branches, a maximum of three cloning events per trajectory are allowed.

In ab initio dynamics, the number of trajectories is severely limited by the high cost of electronic structure calculations (especially for larger molecules). When the initial multi-dimensional wavefunction is randomly sampled with a small number of Gaussians, these Gaussians, due to the high dimensionality of the system, will be located far away from each other with no coupling between them. Applying constraints to run Gaussians closer together would be an inefficient use of computer time unless we need to reproduce a particular quantum effect of the nuclear motion. In this work, we use a simplified semiclassical version of AIMC, where we do not consider the coupling between the trajectories. Instead, each branch simply gets its amplitudes at the time of cloning, and this amplitude determines the statistical weight of that branch.

We initially run an ensemble of 39 Ehrenfest trajectories with the number of branches growing in the process of cloning. All trajectories were propagated for ∼200 fs with a 0.06 fs (2.5 a.u.) timestep. The relatively small number of trajectories and short duration of the dynamics are due to the strict deadline for this work. Nevertheless, despite not very good statistics, our calculations show clear UED patterns for the cyclobutanone photodynamics.

Figure 1 presents the dynamics of the populations for S₀, S₁, and S₂ electronic states. For the first 25 fs of dynamics, the molecules stay in the S₂ state, and then S₂ → S₁ population transfer starts. The growing population of the S₁ state immediately initiates the next step of population transfer, from S₁ into the ground state. Within the next 10 fs, the S₁ state population reaches the equilibrium level of about 20%, when the rates of S₂ → S₁ and S₁ → S₀ transfers are about the same. For the rest of our dynamics, the S₁ state population is fluctuating around this level, while the S₂ state exhibits exponential decay into the ground state S₂ → S₁ → S₀.

The transfer of population between electronic states initiates the process of cloning, which starts at about 50 fs in time and continues uniformly for the duration of the dynamics. By 200 fs, cloning gives rise to 121 branches from the initial 39 trajectories. The small number of initial conditions limits the sampling of the phase-space in the current simulations and is unlikely to result in fully converged calculations. This aspect will be further considered and evaluated in the planned subsequent work.

Figure 2 shows the dynamics of C–C bonds breaking in the cyclobutanone ring. The bond is considered broken when the distance between two atoms exceeds 3 Å. The process of ring opening starts at about 25 fs, simultaneously with the beginning of the non-adiabatic decay of the S₂ state by breaking β-CC bonds. The analysis of trajectories shows that the length of one of β-CC bonds in the molecule often increases toward dissociation from the very beginning of the dynamics, while the lengths of other CC bonds oscillate with a period of about 20–30 fs. Within the next 50 fs, 30% of β-CC bonds break, which corresponds to 60% of the opened rings. At this stage of the dynamics, α-CC bonds are starting to break. In the vast majority of cases, α-CC bonds break in already opened rings, creating ethylene (CH₂)(CH₂) and ethenone (CO)(CH₂) molecules. In some of these ethenone molecules, C=C bonds also later break, creating CO and CH₂ radicals.

After about 100 fs, some opened rings are beginning to close again (or, at least, their ends approach each other to less than 3 Å). Later, the ring can open again, creating an oscillatory behavior in the number of broken bonds.

By the end of the dynamics, the yield in the (CH₂)(CH₂) + (CO)(CH₂) dissociation channel is 40.6%, the yield in the (CH₂)(CH₂) + (CO) + (CH₂) dissociation channel is 3.5%, and the yield of ring opening is 31.0%; also, 17.2% of molecules have remained in the closed ring form. The remaining 7.7% are found at various other intermediate configurations at the end of our 200 fs dynamics; the longer-term dynamics will be a subject of our future work.

To aid our understanding of individual reaction pathways, Fig. 3 shows the bond lengths of the carbon backbone as a function of time for each reaction pathway discussed above; this yields four distinct pathways. Figure 3(a) shows no reaction taking place; all bonds stay sub 2 Å, indicating cyclobutanone remains in the Franck–Condon geometry for the duration of our simulations. There is still the potential for reactions to happen at longer time scales than have been simulated. Panel (b) shows a single bond breaking in the β position of the four-membered ring in a ring-opening reaction. Here, there is an immediate increase in the bond length of a β-CC bond in the S₂ state, confirming what is observed in Fig. 2. A maximum distance of ∼4.5 Å is reached where velocities are reversed, reducing the bond distance to a minimum of ∼3 Å. It is noteworthy that very few simulations underwent an α-CC ring opening process, in contrast to what has been previously observed as the dominant reaction pathway.²⁶ We have, therefore, optimized a minimum energy conical intersection (MECI) using SA(3)-CASSCF(12,12)/aug-cc-pVDZ, shown in the center inset in Fig. 2 displaying the β-CC ring opening. The main reaction pathway observed in our simulations is displayed in Fig. 3(c), with an initial β-CC ring opening reaction observed, followed by a subsequent break of an α-CC bond to yield CH₂CO and C₂H₄ in a stepwise mechanism. From our simulations, we see that this process always takes place in the order of β-CC bond breaking, likely decaying via the MECI shown in Fig. 2, then α-CC bond breaking, which is linked to the S₁/S₀ MECI shown in Ref. 28. Further fragmentation is possible with CH₂CO being able to fragment to produce CO and a CH₂ radical; however, this appears to be a minor channel in our simulations. The trajectory shown in Fig. 3(d) also displays the “tethered” motion shown in Fig. 3(b), causing a rebound effect observed in panel (d) until the full dissociation event has taken place.

To aid in the interpretation in Fig. 4, we have also calculated the static signal, shown in Fig. 5, for all reaction products observed in the AIMC simulations, guided by the reaction pathways in Fig. 3. Representative structures were taken from the trajectories at arbitrary times, showing each reaction observed, and the IAM was then used to calculate the static UED signal for each structure, briefly comprising (a) α-CC bond breaking; (b) β-CC bond breaking; (c) the production of C₂H₄, CH₂, and CO; and also (d) the production of CH₂CO and C₂H₄. The structures for these pathways can be seen in the insets in Fig. 5.

Five key features can be observed in the gas-phase UED of cyclobutanone; these are highlighted by horizontal dashed lines in Fig. 4 at s = 1.1, 2.2, 4.5, 7.2, and 9.8 Å⁻¹. Matching these peaks with those observed in Fig. 5, one can see all four reaction products yield a negative feature at ∼1.1 Å⁻¹ and a positive feature at ∼2.2 Å⁻¹, likely making this peak arise from the carbon backbone of cyclobutanone, the breaking of which is common among all reaction pathways. Both features at ∼1.1 and ∼2.2 Å⁻¹ grow in intensity after ∼20 fs, matching the timescale shown in Sec. IV A. Both features continue to grow in intensity. Most notably, around 75 fs, Fig. 2 shows that we begin to observe α-CC bond breaking. This, coupled with the high intensity observed in the static signal in Fig. 5(c) at 2.2 Å⁻¹, suggests that the stepwise mechanism to form CO, C₂H₄, and CH₂ requires ∼75 fs to form these products. However, we must note that the features in the signal have contributions from all the pathways shown in Fig. 5.

An additional broad positive feature can be observed centered at 4.5 Å⁻¹, with a signal that decays after ∼115 fs. Figure 5(c) shows a broad negative feature between the values ∼3.75 and ∼5.5 Å⁻¹, indicating this reaction pathway causes the depletion of the broad signal. Once again, due to the structural similarities of other reaction pathways, the net signal is a compounded signal with contributions from all products. A similar depletion can be seen in Fig. 5(d). Due to the higher proportion of trajectories being classified as belonging to (d), this likely has a stronger effect on the signal. Depletion of signal in this region is observed in Figs. 5(b)–5(d), having a common feature of β-CC bond breaking, likely making this responsible for the signal depletion.

Further peaks can be seen at 7.5 and 9.8 Å⁻¹. A peak at 7.5 Å⁻¹ can be seen in all pathways with a similar intensity (∼10%); therefore, yielding little structural information other than the molecule has moved away from the equilibrium geometry. In contrast, Fig. 5(c) and, to a lesser extent, Fig. 5(d) both show a signal at 9.8 Å⁻¹. Thus, it is likely this feature arises from the breaking of the β-CC common to both reaction products.

This work was undertaken in response to the “Prediction Challenge: Cyclobutanone Photochemistry” and presents simulated ultrafast electron diffraction (UED) signals for gas-phase cyclobutanone upon photoexcitation into the S₂ electronic state. The main dissociation pathways of photoexcited cyclobutanone were identified with the help of AIMC non-adiabatic dynamics. Then, using these AIMC trajectories, the electronic diffraction was calculated using the IAM method. The calculated UED pattern was compared with static signals for representative structures for the different dissociation pathways observed in the AIMC dynamics. Overall, five key features in the UED pattern can be used to distinguish the reaction channels observed in the AIMC simulation. We find that there is a significant overlap between many features due to a high degree of structural similarity between the different photoproducts, combined with a significant degree of symmetry in cyclobutanone. However, ultimately, we found strong correlations between the timescales and products evident in the simulations and features in the overall UED signal (shown in Fig. 4).

The extent of the work was limited by the strict deadline inherent in the challenge, and it is, therefore, straightforward to identify avenues for further work. To begin with, the presented simulations include only the first 200 fs of the cyclobutanone dissociation dynamics, and the degree of sampling, i.e., the number of trajectories propagated, was also limited. Previously, we developed a technique that allowed us to take the pump pulse shape into account and account for the dynamics during the pulse.⁴⁸ We have not used this approach here and assumed instant excitation, but it can be straightforwardly done. With longer propagation times, it will also be easy to account for some coupling between coherent states using the so called train basis functions.³⁵ This approach does not require additional trajectories or electronic structure calculations. It is anticipated to improve the accuracy of the results, although we do not expect qualitative changes since the underlying trajectories remain the same. All of these improvements will be the subject of subsequent work. Longer simulation times will make it possible to make comparisons to the long-term dynamics observed in the experiment, while more trajectories should, in principle, allow us to move beyond the independent and semiclassical trajectory approximation used when calculating the UED signals. In addition, as discussed in the UED theory section, should the experimental data indicate that more subtle effects in the scattering were observed, then ab initio simulations of the scattering signal, going beyond the independent atom model, are clearly of interest.

In summary, the present work demonstrates the capability of AIMC to simulate photodynamics in a challenging molecule and that UED signals can be predicted straightforwardly from the simulations. We also note that the AIMC simulations should, in principle, provide a better basis for the prediction of experimental signals that reflect the degree of coherence in the molecule during the dynamics, such as nonlinear spectroscopies or coherent mixed scattering.^64,65

D.V.M. and D.V.S. acknowledge the support of EPSRC Grant No. EP/P021123/1. Furthermore, D.S., A.K., and D.M. acknowledge the support from EPSRC Programme Grant No. EP/X026973/1. L.H. and A.K. acknowledge the funding from the Leverhulme Trust (Grant No. RPG-2020-208), and A.K. additionally acknowledges EPSRC Grant Nos. EP/V049240/2, EP/V006819/2, and EP/X026698/1. A.K. also acknowledges the support from the Chemical Sciences, Geosciences, and Biosciences Division, Office of Basic Energy Sciences, Office of Science, US Department of Energy, Grant No. DE-SC0020276. This work was carried out using the computational facilities of the Advanced Computing Research Centre, University of Bristol.

The authors have no conflicts to disclose.

Dmitry V. Makhov: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Lewis Hutton: Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – review & editing (equal). Adam Kirrander: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Dmitrii V. Shalashilin: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.