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

We present the result of our calculations of ultrafast electron diffraction (UED) for cyclobutanone excited into $S_2$ electronic state, which are based on the non-adiabatic dynamics simulations with \textit{Ab Initio} Multiple Cloning (AIMC) method with the electronic structure calculated at the SA(3)-CASSCF(12,12)/aug-cc-pVDZ level of theory. The key features in the UED pattern were identified that can be used to distinguish between the reaction pathways observed in the AIMC dynamics, although there is a significant overlap between representative signals due to structural similarity of the products. The calculated UED pattern can be compared with experiment.


I. INTRODUCTION
Ultrafast electron diffraction (UED) has evolved into a powerful method for structural dynamics. 1,2Although 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 (timedependent) molecular model commensurate with the experimental data, [5][6][7][8][9] the gold standard for interpreting ultrafast experiments remains comparison to high-quality simulations of the photoexcited target molecule.However, such simulations remain challenging, and their veracity depends keenly on numerous methodological choices.As a much needed step towards surveying and evaluating good practice, the Journal of Chemical Physics recently announced the Prediction Challenge: Cyclobutanone Photochemistry to which this paper is a response.
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 a time-resolved UED signals are recorded.][16][17][18][19] In this work, the photodynamics of cyclobutanone is simulated using Ab initio Multiple cloning (AIMC) [20][21][22] 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 quantum dynamics of nuclear wave functions.AIMC was successfully applied before [23][24][25][26][27] to simulate the process of the photodissociation of a number of heterocyclic molecules.Based on AIMC dynamics results, the 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.

A. AIMC
As the AIMC methodology was extensively described before [20][21][22][23] , here we provide only a summary of the technique.The AIMC method represents the further development of the Multi Configurational Ehrenfest (MCE) [28][29][30][31] approach and makes use of the following wave-function ansatz: where R and r are electronic and nuclear coordinates respectively.The electronic part of each basis function is represented in a basis of are adiabatic electronic states |φ I ⟩, and the nuclear parts is a moving Gaussian Coherent State: which is a Gaussian-shaped de Broglie wave centred at R n with momentum P n and phase γ n .The motion of Gaussians |χ n ⟩ is guided by the Ehrenfest force: I , the equations of motion (3) must be solved simultaneously with the equations for a (n) where the right side is the electronic Hamiltonian for nth basis function.Finally, phase γ n is propagated semiclassically as γn = P n Ṙn /2.It is well known that Ehrenfest trajectories misguide basis set when several non-interacting electronic states have significant amplitudes.The cloning procedure is applied in AIMC approach in order to address this issue.In principle, the cloning can be viewed as straightforward way of spawning employed in the Multiple Spawning method 32 .The idea of cloning is to replace a basis function with two clones, each of which is guided in most by just one potential energy surface.In the simplest case of two electronic states: where The total contribution of two clones into wave function |Ψ(R, r,t)⟩ is exactly the same as that of the original basis function.However, cloning increases the size of basis set creating additional flexibility, as two clones can now move in different directions.The cloning is applied when the breaking force ∇V J 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 AIMC method, as it allows AIMC to reproduce the bifurcation of the wave function at conical intersections.
The trajectories in AIMC approach can be calculated independently using potential energies forces and non-adiabatic coupling vectors calculated "on the fly" by an electronic structure code.Then, time-depended Schrödinger equation for amplitudes c n is solved in post-processing in 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 train guided by the same trajectories are among those techniques 33 .It has been demonstrated that MCE can produced the results, which are well converged 34 and AIMC, its ab initio direct dynamics version, is more accurate than Surface Hopping or Ehrenfest Dynamics 35,36 .A technique which allows to take into account pulse shape and dynamics which occur during the excitation has been developed 37 .In its simplest form, which is used in the present paper, AIMC can yield qualitative or semiquantitative picture of the process, similar to that given by Surface Hoping, Ab initio Multiple Spawning (AIMS) 32 and many other popular techniques.

B. Ultrafast electron diffraction
For the modelling of the UED signals, we anticipate that the experimental data in this challenge can be modelled 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, [38][39][40][41] albeit at significantly higher computational cost.
The total (energy-integrated) scattering cross section into the solid angle dΩ at time t is given by, [42][43][44] dσ dΩ where s = k 0 − k 1 is the scattering vector expressed in terms of the wave vectors of the incoming and outgoing electrons.
The scattering is given in units of the Rutherford cross section (dσ /dΩ) Rh which includes the s −4 scaling factor. 45,46Note that the expression above does not account for the duration of the electron pulse, which may be accounted for via a temporal convolution of the predicted signal by the instrument response function for the experiment.General expressions that account for the full wavefunction in Eq. (1), including the non-local nature of the individual Gaussian coherent states, have been derived previously. 42iven the sparse basis used in the present simulations, we resort here to the diagonal bracket-averaged Taylor (BAT) expansion approximation 42 and assume that expansion coefficients are independent of time, c n ≈ c n (t), giving the total scattering intensity as, In this simplified form, sufficient for our present needs, the scattering from each trajectory is given by IAM as, 47 where S inel (s) is the inelastic scattering, which is independent of molecular geometry and isotropic, as underscored by its dependence only on the amplitude of the momentum transfer vector, s = |s|.It is given by an incoherent summation over the individual atomic contributions, with N at the number of atoms in the molecule.The corresponding elastic contribution is given by the form factor F(s, R n (t)), where f e A (s) are the atomic form factors and R nA (t) the position vector for atom A in trajectory n.The form factors for electron scattering are f e A = ( f x A − Z A ), where f x A the xray scattering form factor and Z A is the atomic number. 45,46oth f x A (s) and S A (s) are tabulated. 48For high energy electron scattering it is sometimes necessary to use form factors with relativistic corrections, 49,50 but this is not done presently.
When the target is a gas of anisotropic molecules, |F(s, R n (t))| 2 in Eq. ( 9) is replaced by its rotationally averaged counterpart, ⟨|F(s, where R nAB (t) = |R nA (t) − R nB (t)| is the distance between atoms A and B in trajectory n.

III. COMPUTATIONAL DETAILS
The trajectories were calculated using our own AIMC code, where potential energies, forces and non-adiabatic couplings given "on the fly" by MOLPRO 52 electronic structure package at 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 have been shown to be unimportant for dynamics after excitation into S 2 , therefore we do not include them in our simulations. 53We 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. 19The 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 to the second excited state within Condon approximation.As in our previous simulations 20,[23][24][25][26][27] the cloning thresholds were taken as 5 • 10 −6 a.u. and 2 • 10 −3 a.u.for the magnitude of breaking force and non-adiabatic coupling, respectively.
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 multidimensional wave-function is randomly sampled with a small number of Gaussians, these Gaussians will be located far away from each other with no coupling between them.Running Gaussians closer together would be an inefficient use of CPU 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 brunch.
We initially run an ensemble of 39 Efrenfest trajectories, which give rise to 121 branches in the process of cloning.All trajectories were propagated for 200 fs with 0.06 fs (2.5 a.u.) timestep.A relatively small number of trajectories and short duration of the dynamics is due to the strict deadline for this work.Nevertheless, despite not very good statistics, our calculations show clear UED patterns for the cyclobutanone photodynamics.

A. Dynamics
Figure 1 presents the dynamics of the populations for S 0 , S 1 and S 2 electronic states.For the first 25 fs of dynamics, the molecules stay in S 2 state, then S 2 → S 1 population transfer starts.The growing population of S 1 state immediately initiate the next step of population transfer, from S 1 into the ground state.Within next 10 fs, the S 1 state population reaches the equilibrium level of about 20%, when the rates of S 2 → S 1 and S 1 → S 0 transfers are about the same.For the rest of our dynamics, the S 1 state population is fluctuating around this level, while S 2 state exhibit the exponential decay into a ground state S 2 → S 1 → S 0 .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 time, simultaneously with the beginning of the non-adiabatic decay of S 2 state, by breaking β -CC bonds.Within next 50 fs, 30% of these bonds break, which corresponds to 60% of the opened rings.At this stage of the dynamics, α-CC bonds are starting to break; this happens mostly in already opened rings creating ethylene (CH 2 )(CH 2 ) and ethenone (CO)(CH 2 ) molecules.In some of these ethenone molecules, C=C bonds also later break, creating CO and CH2 radicals.
After about 100 fs time, 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 behaviour in the number of broken bonds.
By the end of the dynamics, the yield in the (CH 2 )(CH 2 ) + (CO)(CH 2 ) dissociation channel is 40.6% , the yield in the (CH 2 )(CH 2 ) + (CO) + (CH 2 ) dissociation channel is 3.5%, and the yield of ring opening is 31.0%;also 17.2% 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.

B. Ultrafast electron diffraction
The AIMC simulations presented in section IV A serve as a framework to calculate the total rotationally averaged UED pattern for cyclobutanone, using the methodology presented in Section II B. The UED signal thus obtained is given in Fig. 3, plotted as percent difference %I tot (s,t), where, I tot (s,t) is the signal at time t and I tot (s, 0) is the reference signal at t = 0, i.e. the pump-off signal.
To aid the interpretation of Fig. 3, we have also calculated the static signal for all reaction products observed in the AIMC simulations, shown in Fig. 4. Representative structures FIG. 3. Gas-phase UED pattern, %I tot (s,t) in Eq. ( 13), for cyclobutanone calculated using the IAM with all trajectories and branches from our dynamics simulations.Five key features in the UED pattern are highlighted with horizontal dashed lines.
were taken from the trajectories showing each reaction observed and the IAM was then used to calculate the UED signal for each structure, briefly comprising of: a) α-CC bond breaking, b) β -CC bond breaking, c) the production of C 2 H 4 , CH 2 and CO, and also d) production of CH 2 CO and C 2 H 4 .The structures for these pathways can be seen in the insets of Fig. 4.
Five key features can be observed in the gas-phase UED of cyclobutanone, these are highlighted by horizontal dashed lines in Fig. 3 at s = 1.1, 2.2, 4.5, 7.2 and 9.8 Å −1 .Matching these peaks with those observed in Fig. 4, one can see all four reaction products yield a negative feature at ∼1.1 Å −1 and a positive feature at ∼2.2 Å −1 .Both features at ∼1.1 and ∼2.2 Å −1 grow in intensity after approximately 20 fs, matching the timescale shown in Section 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 of 4c) at 2.2 Å −1 , suggest that the stepwise mechanism to form CO, C 2 H 4 and CH 2 requires approximately 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. 4.
An additional broad positive feature can be observed centered at 4.5 Å −1 , with a signal that decays after approximately 115 fs. Figure 4c) shows a broad negative feature between the values ∼3.75 Å −1 and ∼5.5 Å −1 , 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. 4d), due to the higher proportion of trajectories being classified as belonging to d) this likely has a stronger effect on the signal.
Further peaks can be seen at 7.5 and 9.8 Å −1 .A peak at 7.5 Å −1 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. 4c) and to a lesser extent Fig. 4d) both show a signal at 9.8 Å −1 .Thus, it is likely this feature arises from the breaking of the β -CC common to both reaction products.

V. CONCLUSIONS
This work was undertaken in response to the "Prediction Challenge: Cyclobutanone Photochemistry" and presents simulated ultrafast electron diffraction (UED) signals for gasphase cyclobutanone upon photoexcitation into the S 2 electronic state.The main dissociation pathways of photoexcited cyclobutanone were identified with the help of AIMC nonadiabatic 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 reac-tion channels observed in the AIMC simulation.We find that there is 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. 3).
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 which allows to include the pump pulse shape into account and allows to account for the dynamics during the pulse 37 .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 21 .This approach does not require additional trajectories and electronic structure cal-culations.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.Also, 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. 54,55KNOWLEDGEMENT DM and DS acknowledge the support of EPSRC grant EP/P021123/1.AK acknowledges funding from the Leverhulme Trust (RPG-2020-208), and EPSRC grant EP/V049240.

FIG. 1 .
FIG. 1.The dynamics of S 2 (red), S 1 (yellow) and S 0 (blue) electronic state populations for cyclobutanone molecule after its photoexcitation into S 2 state.

FIG. 2 .
FIG. 2. The share of broken α-(red) and β -(blue) CC bonds as a function of time for cyclobutanone molecule after its photoexcitation into S 2 state.

FIG. 4 .
FIG. 4. Static signals obtained for representative structures of pathways observed in AIMC simulations given in percent differences.Panel a) α-CC ring opening, b) β -CC ring opening, c) dissociation to form CO, C 2 H 4 and CH 2 , and finally d) CH 2 CO and C 2 H 4 .The structures shown as insets are the geometries from which the static signal is calculated.The five features shown in Figure 3 are highlighted with vertical dashed lines.
V I is potential energy surface of the Ith electronic state, d IJ a non-adiabatic coupling vector, and M is a diagonal matrix of atomic masses.As the force depends on Ehrenfest amplitudes a arXiv:2402.10349v1 [physics.chem-ph]15 Feb 2024 where