# Prediction Through Quantum Dynamics Simulations: Photo-excited Cyclobutanone

Olivia Bennett  Dept. of Chemistry, University College London, 20 Gordon St.,
London WC1H 0AJ, U.K.  Antonia Freibert  Department of Physics, University of
Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany  K. Eryn Spinlove  Dept.
of Chemistry, University College London, 20 Gordon St., London WC1H 0AJ, U.K.
Graham A. Worth<EMAIL_ADDRESS>Dept. of Chemistry, University College
London, 20 Gordon St., London WC1H 0AJ, U.K.

###### Abstract

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

##  I Introduction

Molecular Quantum dynamics (QD) has become a powerful computational tool for
understanding fundamental reactivity. By solving the time-dependent
Schrödinger equation these simulations can follow the time evolution of nuclei
after a molecular system is prepared in a particular way, for example in a
molecular beam or pump-probe experiment. The key feature is that the quantum
nature of the nuclei is taken into account, which is essential for systems
involving either tunnelling or non-adiabatic effects.

A present question, that is to be addressed in this paper, is whether QD
simulations have predictive power, which is necessary for them to become a
mature and trusted tool. This work addresses a challenge issued last year and
aims to predict outcomes of experiments on cyclobutanone conducted at the SLAC
Megaelectronvolt Ultrafast Electron Diffraction facility. These experiments
commenced in late January 2024. The molecules will be photo-excited at 200 nm
with a short (approximately 100 fs wide) pulse and diffraction images will be
recorded over many picoseconds. Up to 200 fs, images will be recorded every 30
fs, and then the delay time will be increased. The challenge to simulation is
to produce signals that can be compared to the experimental outputs and aid in
their interpretation. Ultrafast electron diffraction is proving its value by
providing data on structural dynamics in photochemistry Liu  _et al._ (2020);
Champenois  _et al._ (2021). Examples are the dissociation of water Cao, Wang,
and Zhong (2021), the ring-opening of cyclohexadiene Wolf  _et al._ (2019) and
imaging the passage of CH3I as it passes through a conical intersection Yang
_et al._ (2018).

Cyclobutanones are an interesting target as they exhibit distinct chemical
reactivity compared to cyclic ketones with larger rings, primarily owing to
their inherent ring strain of around 105-120 kJmol-1. Belluš and Ernst (1988);
Bach and Dmitrenko (2006) 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. Trost (1986) Since the mid-1970‘s, 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. Lumbroso  _et al._
(2015); Izquierdo  _et al._ (2005); Lee-Ruff and Wells (2008); Mihovilovic
_et al._ (2005) Reviews of relatively standard unsubstituted cyclobutanone
reactions can be found in Trost (1986); Belluš and Ernst (1988), and, for
example, in a more recent publication it has been reported that an isonitrile-
based, four-component Ugi reaction,Dömling and Ugi (2000) utilising
cyclobutanone, effectively produced an aspartame analogue. Pirrung and Wang
(2009)

The earliest synthetic protocol for the production of cyclobutanone, from
cyclobutanecarboxylic acid, authored by N. Kirscher Kirscher (1905), suffered
from inefficiency, yielded low quantities and featured several sequential
reaction steps. Since then, several high-yield synthetic methods have been
devised. One synthetic scheme involves the epoxidation of
methylidenecyclopropane to 1-oxaspiro[2.2]pentane, leading to a subsequent
lithium-catalyzed rearrangement. Salaun, Champion, and Conia (1988) An
alternative synthetic scheme entails a dialkylation process (utilizing
1-bromo-3-chloropropane) of 1,3-dithiane. This is followed by deprotection of
the ketone through treatment with a mercury salt and cadmium carbonate.
Seebach and Beck (1988)

There has been extensive research into the photo-induced behaviours of
cyclobutanone, and it’s derivatives, both theoretically and experimentally. An
early study of cyclobutanone showed two channels to photolytic decay; the C2
channel, whereby ketene and ethylene are produced, and the C3 channel, whereby
cyclopropane/propene and CO are produced, with a ratio of 2:3. Benson and
Kistiakowsky (1942)

Further theoretical and experimental studies have shown that the ratio of
photoproducts is dependant upon the wavelength of excitation. Upon excitation
between 340-240nm the (n$\pi^{*}$) S1 state is accessed, with the absorption
maximum at around 280nm. Moule (1976); Hemminger, Carless, and Lee (1973) Upon
excitation to this state, three vibrational modes are activated corresponding
to CO stretching, CO out-of-plane wagging and ring puckering. Diau, Kötting,
and Zewail (2001) Studies have indicated that photoexcitation at shorter
wavelengths ($<$315nm range) an IC process via $\alpha$-cleavage occurs
resulting in the production of photofragments with a ratio (C2:C3) of
approximately 2:1. Campbell, Schlag, and Ristow (1967); Tang and Lee (1976) In
the 345-315nm range, an ISC process occurs via a triplet state, followed by
the production of the photofragments, with a ratio (C2:C3) of approximately
1:2 at 326.3 nm Tang and Lee (1976); Hemminger and Lee (1972) and
approximately 1:7 at 343.7 nm. Tang and Lee (1976) 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 kJmol-1 Diau, Kötting, and Zewail
(2001) and 29 kJmol-1. Kao  _et al._ (2020); Xia  _et al._ (2015) By
comparison, the barrier to ring opening of cyclopentanone and cyclohexanone is
in excess of 63 kJmol-1. Diau, Kötting, and Zewail (2001) Experimentally, the
formation of hot ketene fragments from cyclobutanone, in a cyclohexane
solution occurred, at the time delays of 0.25 ps, indicating an ultrafast
ring-opening pathway. The subsequent relaxation of the vibrationally hot
photoproducts and the growth of the ketene fundamental band were fitted to
biexponential functions with shared time constants of 7 ps and 550 ps. Kao
_et al._ (2020) A theoretical study using the Ab-Initio Multiple Spawning
method estimated the S1 lifetime to be around 484 fs with an average time
taken for the $\alpha$-cleavage occurring in 176.6 fs. Liu and Fang (2016) It
should be noted that theoretical findings are in the gas phase.

Experimental and theoretical studies have also explored the second excited
(S2) state of cyclobutanone, and other cyclic ketones. In the Franck-Condon
region the S2 state exhibits Rydberg 3s character and undergoes an IC process,
facilitated by the low frequency ring puckering mode, to the S1 (n$\pi^{*}$)
state. Kuhlman, Sølling, and Møller (2012) 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 . Kuhlman  _et al._ (2012);
Kuhlman, Sølling, and Møller (2013) 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$\rightarrow$n$\pi^{*}$ process in cyclobutanone, cyclopentanone and
cyclohexanone were determined to be 13:2:1. Kuhlman, Sølling, and Møller
(2012)

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
simulation of the experimental signal and relate it to the evolving molecular
geometry.

##  II Methodology and Computational Details

All quantum chemistry calculations were performed using the Molpro 2022
package Werner  _et al._ (2022, 2012, 2020). The MCTDH and DD-vMCG simulations
used a development version of the Quantics Suite Worth  _et al._ (2023); Worth
(2020), with the Vibronic Coupling Hamiltonians built using the VCHam package
Cattarius, Markmann, and Worth (2007)

###  II.1 Quantum Chemistry Calculations

The starting point is to choose the level of electronic structure theory able
to describe the system of interest, here cyclobutanone, balancing accuracy
with cost. We start by finding and characterising the ground-state equilibrium
structure. For this the coupled-cluster singles and doubles (CCSD) method was
used with a 6-311++G** basis sets.

The next step is to choose the level of theory needed to describe the excited-
states. To describe potential long range motion, a multi-configurational
method is best and as cyclobutanone is quite small with only five heavy atoms,
a complete active space self-consistent filed (CASSF) wavefunction is
possible. After trying a number of different complete active space (CAS) sets
of orbitals, it was decided that a CAS with 8 electrons in 10 orbitals would
be suitable. The orbitals include the $\pi$ orbitals and low lying Rydberg
orbitals. The CAS orbitals are shown in Fig. 1. After further tests
calculating excitation energies at the $\mathrm{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.

Figure 1: The 6-311++G** CAS(8,10) orbitals of cyclobutanone where orbitals 19
and 20 are the HOMO and LUMO respectively.

###  II.2 Vibronic Coupling Model Hamiltonian

Assuming a diabatic electronic basis, the Hamiltonian can be written in matrix
form as

$\mathbf{H}=\hat{T}\mathbf{1}+\mathbf{W}$ (1)

where $\hat{T}$ is the kinetic energy operator, $\mathbf{1}$ a unit matrix and
$\mathbf{W}$ the diabatic potential matrix. Thus in this representation the
kinetic energy operator is diagonal in the electronic states and the nuclear-
electronic, non-adiabatic, coupling is represented by potential-like functions
as the off-diagonal elements of the potential matrix. The key feature is that,
unlike the usual adiabatic representation, the non-adiabatic couplings do not
contain singularities Köppel, Domcke, and Cederbaum (1984); Worth and
Cederbaum (2004).

The Vibronic Coupling Model represents the diabatic potentials using Taylor
expansions around the Franck-Condon (FC) point and truncated to a low order.

$\mathbf{H}=\mathbf{H}_{0}+\mathbf{W}^{(0)}+\mathbf{W}^{(1)}+\mathbf{W}^{(2)}\ldots$
(2)

In the standard Linear Coupling (LVC) scheme emplayed here, the zero-order
Hamiltonian is taken to be the ground-state Hamiltonian in the Harmonic
approximation.

$\mathbf{H}_{0}=\sum_{\alpha}\frac{\omega_{\alpha}}{2}\left(-\frac{\partial^{2}}{\partial
Q_{\alpha}^{2}}+Q_{\alpha}^{2}\right)\mathbf{1}$ (3)

Here, the coordinates used are mass-frequency scaled normal modes with
$\omega_{\alpha}$ representing the ground-state frequencies. The zero-order
diabatic matrix

$W_{ij}^{(0)}=E_{i}\delta_{ij}$ (4)

is a diagonal matrix with the energies of the states at the Franck-Condon
point. The first-order expansion matrices

$\displaystyle W^{(1)}_{ii}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\kappa^{(i)}_{\alpha}Q_{\alpha}$ (5) $\displaystyle
W^{(1)}_{ij}$ $\displaystyle=$
$\displaystyle\sum_{\alpha}\lambda^{(ij)}_{\alpha}Q_{\alpha}$ (6)

contain the linear on-diagonal coupling parameters, $\kappa^{(i)}_{\alpha}$
along with the off-diagonal parameters, $\lambda^{(ij)}_{\alpha}$. The former
are the gradients of the potential surfaces at the Franck-Condon point, while
the latter are the non-adiabatic coupling between states $i$ and $j$. Higher-
order matrices contain only on-diagonal terms as required to model the
diabatic surfaces, hence the name linear vibronic coupling. Terms up to
fourth-order were required for some modes due to the anharmonicity of the
surfaces.

The parameters are only non-zero if they obey the symmetry rule that the
product of normal modes and states involved in the expansion term must contain
the totally symmetric representation. For example, the linear on-diagonal
terms in terms of the related Hamiltonian matrix element is

$\kappa^{(i)}_{\alpha}=\langle\psi_{i}|\frac{dH_{\mathrm{el}}}{dQ_{\alpha}}|\psi_{i}\rangle_{Q_{\alpha}}$
(7)

where $H_{\mathrm{el}}$ is the usual electronic Hamiltonian with $\psi_{i}$,
the electronic wavefunction at the Franck-Condon geometry, and the derivative
is also evaluated at this point. This matrix element is only non-zero if
$\Gamma_{i}\times\Gamma_{i}\times\Gamma_{\alpha}\supset A_{1}$, where
$\Gamma_{i}$, $\Gamma_{\alpha}$ are the irreps of the electronic states and
normal mode, respectively, and $A_{1}$ the totally symmetric irrep. Thus only
totally symmetric vibrations have non-zero $\kappa$ parameters. In a similar
way, $\lambda$ parameters are only non-zero if
$\Gamma_{i}\times\Gamma_{j}\times\Gamma_{\alpha}\supset A_{1}$ and thus only
vibrations with the correct symmetry can couple states. These rules are easily
extended to higher order terms.

###  II.3 Multi-Configuration Time-Dependent Hartree (MCTDH)

The MCTDH method provides a variational solution to the time dependent
Schrödinger equation (TDSE) using the wavefunction ansatz

$\Psi(\mathbf{Q},t)=\sum_{j_{1}\ldots j_{p}s}A_{j_{1}\ldots
j_{p}s}\varphi^{(1)}_{j_{1}}(\mathbf{q}_{1},t)\ldots\varphi^{(p)}_{j_{p}}(\mathbf{q}_{p},t)|s\rangle$
(8)

i.e. a full direct product expansion of the wavefunction in terms of $p$ sets
of low-dimensional basis functions, $\varphi_{i}(\mathbf{q_{i}},t)$, known as
single-particle functions (SPFs). These functions are time-dependent and
depend on a set of physical coordinates
$\mathbf{q}_{i}=(Q_{\alpha},Q_{\beta},\ldots)$. They are represented in a
time-independent primitive basis, $\chi_{r}$

$\varphi^{(\kappa)}_{i}(\mathbf{q}_{\kappa},t)=\sum_{r}c_{ri}(t)\chi_{r}(\mathbf{q}_{\kappa})$
(9)

which provides an underlying grid. The functions $|s\rangle$ are time-
independent vectors that define the electronic state populated. This is
usually referred to as the single-set formulation of MCTDH as a single set of
SPFs are used to describe all states.

The method is thus a contraction scheme from the full basis set to a
variational one. The resulting equations of motion for the time evolution of
the SPFs and the expansion coefficients $A_{j_{1}\ldots j_{p}s}$ are well
described in the literature Beck  _et al._ (2000); Meyer, Gatti, and Worth
(2009) and will be not be given here.

For large systems, the Multi-Layer Multi-Configuration Time-Dependent Hartree
(ML-MCTDH) variant must be used Wang and Thoss (2003); Manthe (2008); Vendrell
and Meyer (2011). In this, the SPFs are expanded in the form of the MCTDH
ansatz Eq. (8). These new basis functions (SPFs of the first layer) can in
turn be expanded in this form. This procedure is continued to provide layers
of functions with a set of primitive grid functions forming the lowest layer.
In this way, a full tensor contraction scheme is set up and used to
variationally solve the TDSE.

The ML-MCTDH method was used in this work to allow simulations including all
degrees of freedom. The Hamiltonian is provided by the vibronic model
described above. The primitive basis functions used for all coordinates were
harmonic oscillator discrete variable representations (DVRs) Beck  _et al._
(2000). Calculations need to be converged with respect to the basis sets. In
order to achieve this, the size of the primitive basis is checked to ensure
the wavepacket does not significantly populate the end grid points (population
less than $10^{-6}$). The lowest natural populations of the SPFs are also kept
below $10^{-3}$. The latter was ensured throughout the calculation by
dynamically growing the basis in any layer whenever the population of the
least important SPF approaches this limit Mendive-Tapia  _et al._ (2017).

The MCTDH calculations using the vibronic coupling Hamiltonian are able to
simulate the absorption spectrum from the Fourier Transform of the
autocorrelation function

$I(\omega)\sim\int_{0}^{T}dt\,g(t)\langle\Psi(0)|\Psi(t)\rangle\exp(i\omega
t)$ (10)

The function $g(t)$ is used to remove artefacts from the Fourier Transform.
The chosen form is

$g(t)=\cos\left(\frac{\pi t}{T}\right)\times\exp\left(-\frac{t}{\tau}\right)$
(11)

where the first term is to ensure that the autocorrelation function goes to
zero at the end of the simulation, time $T$, and the second term a damping
function. The damping time used was 150 fs.

In addition, diabatic state populations are directly obtained from the
wavefunction which give a rate of electronic relaxation as the molecule
changes electronic configuration. The short-time dynamics may also be analysed
for the evolution of the molecular geometry in terms of excitation of the
vibrational coordinates.

###  II.4 Direct Dynamics variational Multi-Configuration Gaussian (DD-vMCG)

To create a more flexible wavefunction that has a direct connection to
molecular structures, the wavefunction ansatz can be written as a
superposition of Gaussian functions

$\Psi(\mathbf{q},t)=\sum_{js}A_{js}G_{j}(\mathbf{q},t)|s\rangle$ (12)

with each Gaussian function a separable product of one-dimensional Gaussian
functions

$G_{j}(\mathbf{q},t)=g^{(1)}_{j}(q_{1},t)\ldots g^{(f)}_{j}(q_{f},t)$ (13)

having the form

$g^{(\kappa)}_{j}(q_{\kappa},t)=\exp\left(\zeta_{j\kappa}q_{\kappa}^{2}+\xi_{j\kappa}q_{\kappa}+\eta_{j}\right)$
(14)

where $\zeta_{j\kappa},\xi_{j\kappa}$ and $\eta_{j}$ are quadratic, linear and
scalar parameters respectively.

Solving the TDSE using this ansatz, and the Dirac-Frenkel variational
principle, leads to the variational Multiconfigurational Gaussian (vMCG)
method Burghardt, Meyer, and Cederbaum (1999); Richings  _et al._ (2015). The
expansion coefficients evolve with equations similar to those for MCTDH.
However, the non-orthogonality of the GWP basis must be taken into account

$i\dot{A}_{is}=\sum_{jk}S^{-1}_{ij}\left(H_{jk}-i\tau_{jk}\right)A_{k}$ (15)

with the overlap matrix, $\mathbf{S}$, Hamiltonian matrix, $\mathbf{H}$ and
overlap time-derivative, $\boldsymbol{\tau}$

$\displaystyle S_{ij}$ $\displaystyle=$ $\displaystyle\langle
G_{i}|G_{j}\rangle$ (16) $\displaystyle H_{ij}$ $\displaystyle=$
$\displaystyle\langle G_{i}|H|G_{j}\rangle$ (17) $\displaystyle\tau_{ij}$
$\displaystyle=$ $\displaystyle\langle G_{i}|\dot{G}_{j}\rangle$ (18)

In the standard approach, the widths of the Gaussians, $\zeta$, are kept fixed
(frozen Gaussians) and the scalar parameter, $\eta$ is fixed by the
requirement for the Gaussians to be normalised and for the phase to be kept
zero. The time-dependence is then carried by the linear parameters, $\xi$.
Collecting the set of parameters for a multi-dimensional Gaussian into a
vector, $\Lambda_{j\alpha}=\xi_{j\alpha}$ , the equations of motion (EOM) for
the Gaussian functions can be written as

$i\dot{\boldmath\Lambda}=\boldmath{X}+\boldmath{C}^{-1}\boldmath{Y}_{R}$ (19)

Using the relationship between the general form of Eq. 14 and the Gaussian
wavepackets of Heller Heller (1975), the linear parameters can be written in
terms of the coordinate and momentum at the centre of the Gaussian

$\dot{\Lambda}_{i\alpha}=\dot{\xi}_{i\alpha}=-2\zeta_{i}\dot{q}_{\alpha}+i\dot{p}_{\alpha}$
(20)

and the vector $\boldmath{X}$ in Eq. (19) is related to the classical
equations of motion

$X_{i\alpha}=-2\zeta^{(\kappa)}_{i}\frac{{p}_{\alpha}}{m_{\alpha}}-i\left.\frac{\partial
V}{\partial q_{\alpha}}\right|_{\mathbf{q}_{i}}$ (21)

Thus the Gaussians move along trajectories that have a classical component
with an additional variational coupling between leading to faster convergence
of the wavefunction.

Due to the localised nature of a Gaussian function, it is reasonable to
calculate integrals of the potential energy using a Local Harmonic
Approximation (LHA) in which the potential is expanded to second order around
the centre of a Gaussian

$V(\mathbf{q})=V(\mathbf{q}_{i})+\sum_{\alpha}\left.\frac{\partial V}{\partial
q_{\alpha}}\right|_{\mathbf{q}_{i}}(q_{\alpha}-q_{i\alpha})+\frac{1}{2}\sum_{\alpha\beta}\left.\frac{\partial^{2}V}{\partial
q_{\alpha}\partial
q_{\beta}}\right|_{\mathbf{q}_{i}}(q_{\alpha}-q_{i\alpha})(q_{\beta}-q_{i\beta})$
(22)

Solutions to the TDSE no longer converge on the exact result, but the
integrals can all be performed analytically and the algorithm can be used for
direct dynamics simulations.

Direct dynamics simulations use this formalism to calculate the potential from
quantum chemistry calculations on-the-fly. The LHA requires the energies,
gradients and Hessians at the Gaussian central coordinate, and these can all
be provided by a quantum chemistry calculation. In the standard DD-vMCG
protocol, these quantities are calculated and stored in a database. New points
are added to the database only when a GWP has coordinates that are
significantly different from any point stored in the database. The surfaces
experienced by the evolving GWPs are provided by Shepard interpolation between
the points in the database. Calculations are run in the diabatic picture, and
the potentials are diabatised using the propagation diabatisation scheme. For
full details see Ref. Christopoulou, Freibert, and Worth (2021).

The major effort in the vMCG method is the inversion of the $\mathbf{C}$
matrix in Eq. (19) which has the size of $(N_{f}\times n)^{2}$, where $N_{f}$
is the number of degrees of freedom and $n$ the number of basis functions.
This effort can be reduced by partitioning the system and making the
wavefunction a multi-configurational product from two (or more) sets of
functions.

$\Psi(\mathbf{Q},t)=\sum_{j_{1}\ldots j_{p}s}A_{j_{1}\ldots
j_{p}}G^{(1)}_{j_{1}}(\mathbf{q}_{1},t)\ldots
G^{(p)}_{j_{p}}(\mathbf{q}_{p},t)|s\rangle$ (23)

This is the G-MCTDH ansatz Burghardt, Meyer, and Cederbaum (1999). The EOM are
the same as Eqs. (15, 19) but now the function $G_{i}$ is a configuration,
i.e. a product of basis functions,
$G_{i}=G^{(1)}_{i_{1}}G^{(2)}_{i_{2}}\ldots$ and the propagation of the
Gaussian parameters include mean-field operators connecting the different
partitions.

###  II.5 Gas Phase Ultrafast Electron Diffraction (GUED) Signal

GUED uses high energy pulses of electrons to provide structural information in
the form of a diffraction pattern. Electron scattering by a potential field (a
molecule) can be described by measuring the momentum transfer from an electron
to the target on collision. For a set of atoms under the independent atom
model, this can be formulated to give an intensity as a function of scattered
momentum transfer $s$ that has two parts: an “atomic” due to scattering
directly off an atom and a “molecular” due to interference from scattering off
neighbouring atoms Centurion, Wolf, and Yang (2022). The atomic part is

$I_{\mathrm{at}}(s)=\sum_{i=1}^{N}|f_{i}(s)|^{2}$ (24)

and the molecular part is

$I_{\mathrm{mol}}(s)=\sum_{i=1}^{N}\sum_{j\neq
i}^{N}f_{i}(s)f_{j}(s)\frac{\sin(sr_{ij})}{sr_{ij}}$ (25)

where $N$ is the number of atoms, $f_{i}$ are the “atomic form factors”
accounting for the structure of an atom in terms of nuclei and electrons and
$r_{ij}$ denotes the inter-atomic distance between the $i^{\mathrm{th}}$ and
$j^{\mathrm{th}}$ atom. The total intensity is simply the sum of these two
terms. Actual measurements record the scattering signal as a modified
scattering intensity

$sM(s,t)=s\frac{I_{\mathrm{mol}}(s,t)}{I_{\mathrm{at}}(s)}$ (26)

where the time dependence has been added to indicate that the signal is
obtained by scattering off the sample at time $t$. It is more useful to record
the signal as the difference between the signal at $t$ and an initial
(reference) signal

$\Delta sM(s,t)=s\frac{\Delta I_{\mathrm{mol}}(s,t)}{I_{\mathrm{at}}(s)}$ (27)

as the change in signal gives the structural changes in the molecule without
needing to extract $I_{\mathrm{mol}}$. This signal can be related to the
change in “pair distribution function” (PDF) by integrating over $s$ to give
the signal as a function of inter-atomic distances

$\Delta P(r,t)\approx r\int_{s_{\mathrm{min}}}^{s_{\mathrm{max}}}\Delta
sM(s,t)\sin(sr)e^{-\alpha s^{2}}ds.$ (28)

where $\alpha$ is a smoothing factor. Changes in the molecular geometry can
thus be seen in the PDF as losses and rises as atoms move apart. It is,
however, summed over all the pairs in a molecule so each inter-nuclear
distance may be related to a number of different atom pairs, making
interpretation with a molecular simulation difficult. From a simulation, the
PDF can be simulated by using the atomic coordinates along a trajectory to get
first the molecular intensity, Eq. (25), and then obtain the modified
scattering intensity before calculating the difference PDF of Eq. (28).

Finally, it needs to be taken into account that a molecule is not a classical
object moving along a single trajectory. The scattering takes place from the
molecular density given by the evolving wavepacket. This can be taken into
account by integrating over the molecular configurations

$\bar{I}_{\mathrm{mol}}(s)=\int
d\mathbf{q}\,\Psi^{2}(\mathbf{q})I_{\mathrm{mol}}(s)$ (29)

where $\Psi^{2}(\mathbf{q})$ is the probability of the molecule being in
configuration $\mathbf{q}$.

In vMCG, a useful approximation for calculating the expectation values of
operators that depend only on the atomic coordinates is

$\displaystyle\Psi^{2}(\mathbf{q})$ $\displaystyle=$
$\displaystyle\sum_{ij}0.5(A^{\ast}_{i}\langle
G_{i}\mathbf{(}\mathbf{q}_{i})|G_{j}(\mathbf{q_{j}})\rangle
A_{j}O(\mathbf{q}_{j})+c.c.)$ (30) $\displaystyle=$
$\displaystyle\sum_{j}\mathrm{GGP}_{j}O(\mathbf{q}_{j})$ (31)

where $\mathbf{r_{j}}$ is the molecular configuration defined by the centre
coordinate of the GWP, and $O(\mathbf{q}_{j})$ the value of the operator at
that point. The expectation value is thus approximated by a weighted sum along
the GWP trajectories, with the weighting given by the Gross Gaussian
Populations

$\mathrm{GGP}_{j}=Re\sum_{i}A^{\ast}_{i}S_{ij}A_{j}$ (32)

where $S_{ij}$ is the GWP overlap matrix. The final difference PDF from the
signal calculated over the molecular density can thus be approximated in this
way as a weighted sum

$\bar{I}_{\mathrm{mol}}\approx\sum_{i}\mathrm{GGP}_{i}I_{\mathrm{mol,i}}$ (33)

where $I_{\mathrm{mol,i}}$ is the molecular scattering intensity from vMCG
trajectory $i$. In fact, this can be carried through the procedure for
calculating the difference PDF and one can provide simply a weighted sum of
difference PDFs.

The modified scattering intensity for each trajectory at each time step was
calculated using a code to simulate electron diffraction signals provided by
Wolf et al.Wolf and Martinez (2021) with subsequent weighting according to Eq.
(33) to get the final time-dependent signal. The time-dependent difference
pair distribution function $\Delta P(r,t)$ was then obtained using
$s_{\mathrm{min}}=0~{}\mathring{A}^{-1}$,
$s_{\mathrm{max}}=10~{}\mathring{A}^{-1}$ and
$\alpha=0.036~{}\mathring{A}^{2}$ matching the experimental spatial resolution
of $0.6~{}\mathring{A}$.

##  III Results

###  III.1 Ground-State Equilibrium Structure

Using CCSD/6-311++G**, two stable structures of cyclobutanone were found. The
first is planar with $\mathrm{C_{2v}}$ symmetry. The second is bent and has
$\mathrm{C_{s}}$ symmetry. These are shown in Fig. 2 and the coordinates,
along with bond lengths and angles, are given in the Supplementary
Information. The ring puckering angle between the C1C2C3 and C4C2C3 planes was
found to be 171.7∘, and the angle of the C–O bond to the C1C2C3 plane was
12.3∘. Frequency calculations showed that the $\mathrm{C_{s}}$ structure is a
minimum and the $\mathrm{C_{2v}}$ structure a transition state. The
frequencies of the vibrational modes of both structures are listed in the
Supplementary Information. With one exception, the transition mode, changes in
the normal modes and frequencies are minimal between the two structures.

The transition mode with an imaginary frequency of 80 cm-1 is the ring
puckering motion, $\nu_{1}$ and at $\mathrm{C_{s}}$ it has a real frequency of
100 cm-1. The difference in energy between the structures, $\Delta
E=E_{\mathrm{C_{2v}}}-E_{\mathrm{C_{s}}}$ is the barrier height to the out-of-
plane bend and has a value of 0.0087 eV. The frequency of the transition mode,
$\nu_{1}$ means it has a zero point energy level at 0.0062 eV. Thus even
though the barrier is very small, the low frequency of the transition mode
means that it is a stable structure and at 0 K the equilibrium structure will
be $\mathrm{C_{s}}$. However, at room temperature, a significant population of
molecules will be in excited vibrational levels and the ground-state structure
will be closer to $\mathrm{C_{2v}}$. It was also found that at the CASSCF and
CASPT2 levels of theory used in the dynamics calculations the C2v structure is
the minimum energy on the ground-state. We therefore will need to examine the
dynamics starting from both structures.

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

###  III.2 Excitation Energies and Absorption Spectra

The Excitation energies and characters from both structures to the lowest two
singlet and three triplet states are listed in Table 1. Calculations used a
CAS with 8 electrons in 10 orbitals and a 6-311++G** basis. The orbitals are
shown in Fig. 1.

Little difference is found for the excitation from the two structures. In
agreement with previous work Diau, Kötting, and Zewail (2001); Kuhlman,
Sølling, and Møller (2013); Xia  _et al._ (2015), the lowest singlet state is
an excitation from the in-plane $\pi_{y}^{\ast}$ to the $\pi^{\ast}$ 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 T1 and T2 states are the same
configurations as S1 and S2, but lie slightly lower in energy. The T3 state
lies above S2 at 6.7 eV and is a $\pi$ to $\pi^{\ast}$ excitation.

Table 1: Excitation Energies (eV) of cyclobutanone calculation at the CASSCF and CASPT2 levels of theory with a CAS(8,10)/6-311++G** wavefunction. The symmetry labels refer to the structure which has been optimised using CCSD/6-311++G** |  | $\mathrm{C_{2v}}$ | $\mathrm{C_{s}}$ |
---|---|---|---|---
State | Character | CASSCF | CASPT2 | CASSCF | CASPT2 | Expt
S${}_{0}(A_{1})$ |  | 0.0 | -0.0087 111Energy calculated CCSD/6-311++G** |
S${}_{1}(A_{2})$ | $\pi_{y}^{\ast}\rightarrow\pi_{x}^{\ast}$ | 4.93 | 4.16 | 5.19 | 4.21 | 3.61 222Adiabatic 0-0 transition Ref. Tang and Lee (1976)
S${}_{2}(B_{2})$ |  $\pi_{y}^{\ast}\rightarrow$ sRy | 5.81 | 6.32 | 5.86 | 6.23 | 6.11 333Adiabatic 0-0 transition Ref. O’Toole  _et al._ (1991)
T${}_{1}(A_{2})$ | $\pi_{y}^{\ast}\rightarrow\pi_{x}^{\ast}$ | 4.59 | 3.77 | 4.58 | 3.84 |
T${}_{2}(B_{2})$ |  $\pi_{y}^{\ast}\rightarrow$ sRy | 5.73 | 6.27 | 5.85 | 6.18 |
T${}_{3}(A_{1})$ | $\pi_{x}\rightarrow\pi_{x}^{\ast}$ | 6.44 | 6.69 | 6.48 | 6.24 |

###  III.3 Vibronic Model Hamiltonian

The vibronic coupling model Hamiltonian of Sec. II.2 was obtained using the
mass-frequency scaled normal modes of the ground-state C2v structure, with the
imaginary frequency of $\nu_{1}$ taken as real. The ground-state and
excitation energies were calculated at the CASPT2 corrected
SA6-CAS(8,10)/6-311++G** level of theory at twenty one points along each
vibrational mode, treating the singlet and triplet manifolds separately. The
parameters of the model were then chosen to optimise the fit between the
adiabatic potential surfaces of the model and the calculated adiabatic
energies Cattarius  _et al._ (2001). 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 added as constant
values between the singlet and triplet states, using the magnitude of the
complex components to give an effective coupling treating each triplet state
as a single state. The excitation energies are given in Table 1. The spin-
orbit couplings, given as a normalised average over the three components of
the triplet state are given in Table 2.

Table 2: Averaged Spin-Orbit Coupling (cm-1) for cyclobutanone at the CASSCF level of theory with a SA6-CAS(8,10)/6-311++G** wavefunction State | T${}_{1}(A_{2})$ | T${}_{2}(B_{2})$ | T${}_{3}(A_{1})$
---|---|---|---
S${}_{0}(A_{1})$ | 33.94 | 1.90 | 0.00
S${}_{1}(A_{2})$ | 0.00 | 0.76 | 15.42
S${}_{2}(B_{2})$ | 0.64 | 0.00 | 1.45

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 3 and 4 list the linear $\kappa$ and
$\lambda$ parameters. Coupling strengths are defined as
$\kappa^{(i)}_{\alpha}/\omega_{\alpha}$ and
$\lambda^{(ij)}_{\alpha}/\omega_{\alpha}$, where $\omega_{\alpha}$ is the
frequency of the mode. Using these strengths as a criteria, a subset of 12 key
modes can be identified. $\nu_{5},\nu_{16},\nu_{19},\nu_{20}$ and $\nu_{21}$
are the totally symmetric vibrations with the highest $\kappa$ coupling
strengths, i.e. they have significant gradients on the potential surfaces at
the Franck-Condon point. Modes $\nu_{4}$ and $\nu_{9}$, with A2 symmetry carry
the coupling between S0 and S1. Modes $\nu_{3}$ and $\nu_{17}$, with B2
symmetry couple between S0 and S2, while modes $\nu_{2}$ and $\nu_{25}$, with
B1 symmetry carry the coupling between S1 and S2. Finally, $\nu_{1}$ has
symmetry B1 and dominates the coupling between T1 and T2. The coupling between
T1 and T3 is dominated by $\nu_{4}$, and between T2 and T3 by $\nu_{11}$. This
latter coupling though is quite weak and $\nu_{11}$ is not included in the key
modes as the high energy T3 state is unlikely to be important for the short
term dynamics. The key vibrations are plotted in Fig. 3.

The parameters for the full final Hamiltonian are given as the Quantics
operator file in the supplementary datasets. Cuts through the potential
surfaces for the full singlet and triple manifold are also given in the SI. In
the triplet manifold, the major features are due to the near degenerate T2 and
T3 states, which are also near degenerate with the S2 state. The spin-orbit
coupling is however weak. In this manifold, $\nu_{2},\nu_{3}$ and $\nu_{16}$
show double well structures in T2, while in $\nu_{21}$ the T2 and T3 states
cross at the Franck-Condon point and T2 looks to be dissociative. $\nu_{21}$
is the C-O stretch vibration and is also the most interesting mode in the sing
manifold. The S2 / S1 crossing is seen to negative values along this modes and
in the downhill direction from the FC point,. The S1 state is also long and
fairly flat, but not dissociative, out to positive values. The mode $\nu_{5}$,
which is the ring stretching along the C-C-O axis is also potentially
interesting as it has a strong gradient at the FC point in all excited states.

Table 3: On-diagonal linear parameters from the Vibronic Coupling Model of cyclobutatone calculated at the CASPT2 level with a SA6-CAS(10,8)/6-311++G** wavefunction. Mode | Symmetry |  $\kappa^{(2)}$ (eV) |  $\kappa^{(3)}$ (eV) |  $\kappa^{(4)}$ (eV) |  $\kappa^{(5)}$ (eV) |  $\kappa^{(6)}$ (eV)
---|---|---|---|---|---|---
$\nu_{5}$ | A1 | 0.0727 | 0.0272 | 0.0435 | 0.0423 | 0.0571
$\nu_{7}$ | A1 | 0.0509 | -0.0347 | 0.0186 | -0.0118 | -0.0118
$\nu_{10}$ | A1 | 0.0135 | 0.0331 | 0.0445 | 0.0116 | 0.0116
$\nu_{16}$ | A1 | -0.0577 | -0.1895 | -0.0197 | -0.0728 | 0.0073
$\nu_{19}$ | A1 | -0.0419 | -0.1246 | -0.0370 | -0.0970 | -0.0155
$\nu_{20}$ | A1 | 0.0508 | 0.1094 | 0.0377 | 0.1110 | 0.0475
$\nu_{21}$ | A1 | -0.3530 | 0.1541 | -0.3070 | -0.7209 | 0.0881
$\nu_{23}$ | A1 | -0.0008 | -0.0047 | -0.0098 | -0.0884 | -0.0884
$\nu_{24}$ | A1 | -0.0174 | -0.0251 | -0.0005 | 0.0007 | 0.0007
Table 4: Off-diagonal linear parameters from the Vibronic Coupling Model of cyclobutatone calculated at the CASPT2 level with a SA6-CAS(10,8)/6-311++G** wavefunction. Modes with all parameters less than $10^{-2}$ are omitted. Mode | Symmetry |  $\lambda^{(1,2)}/\omega$ (eV) |  $\lambda^{(1,3)}$ (eV) |  $\lambda^{(2,3)}$ (eV) |  $\lambda^{(4,5)}$ (eV) |  $\lambda^{(4,6)}$ (eV) |  $\lambda^{(5,6)}$ (eV)
---|---|---|---|---|---|---|---
$\nu_{1}$ | B1 | 0.0000 | 0.0000 | 0.0006 | 0.0569 | 0.0000 | 0.0000
$\nu_{2}$ | B1 | 0.0000 | 0.0000 | -0.0277 | -0.0035 | 0.0000 | 0.0000
$\nu_{3}$ | B2 | 0.0000 | -0.2513 | 0.0000 | 0.0000 | 0.0000 | -0.0000
$\nu_{4}$ | A2 | -0.1699 | 0.0000 | 0.0000 | 0.0000 | 0.0489 | 0.0000
$\nu_{8}$ | B2 | 0.0000 | 0.0149 | 0.0000 | 0.0000 | 0.0000 | 0.0000
$\nu_{9}$ | A2 | 0.1913 | 0.0000 | 0.0000 | 0.0000 | -0.0476 | 0.0000
$\nu_{11}$ | B2 | 0.0000 | -0.0307 | 0.0000 | 0.0000 | 0.0000 | 0.0364
$\nu_{12}$ | B1 | 0.0000 | 0.0000 | -0.0019 | -0.0006 | 0.0000 | 0.0000
$\nu_{14}$ | B1 | 0.0000 | 0.0000 | 0.0006 | 0.0527 | 0.0000 | 0.0000
$\nu_{17}$ | B2 | 0.0000 | 0.1293 | 0.0000 | 0.0000 | 0.0000 | -0.0172
$\nu_{18}$ | B2 | 0.0000 | 0.0739 | 0.0000 | 0.0000 | 0.0000 | 0.0319
$\nu_{22}$ | B2 | 0.0000 | -0.0701 | 0.0000 | 0.0000 | 0.0000 | -0.0000
$\nu_{25}$ | B1 | 0.0000 | 0.0000 | 0.1351 | -0.0761 | 0.0000 | 0.0000
$\nu_{26}$ | A2 | 0.0069 | 0.0000 | 0.0000 | 0.0000 | -0.0285 | 0.0000
$\nu_{27}$ | B1 | 0.0000 | 0.0000 | -0.0002 | -0.0543 | 0.0000 | 0.0000
Figure 3: The 12 key normal modes of cyclobutanone, calculated at
CCSD/6-311++G** level on the Ground-state C2v structure. Arrows indicate the
motion of the atoms during each vibration.

###  III.4 Importance of Triplet States

ML-MCTDH simulations were run for 200 fs on the full singlet and triplet
vibronic coupling operator of Sec. III.3. The tree defining the layering along
with the number of functions is given in the SI. Simulations started in both
the C2v and Cs structures. At the end of this, less than 0.1 % population
flows into the T1 state, and the populations of the T2 and T3 are even
smaller. This is not surprising given the small size of the coupling - for
both the S1 and S2 states there is zero SOC with the close lying triplet (T1
and T2, respectively).

Of course the vibronic coupling model treats the SOC with only a constant
value taken at the Frank-Condon point. To see if more ISC occurs if the
nuclear motion is taken into account, DD-vMCG simulations were run including
all 6 states and calculating the SOC at each point added to the DB. All
properties were calculated at the CAS(8,10) level. The wavepacket was
partitioned into 2 parts with the 12 key modes in one partition and treated
with 16 GWPs and the remaining modes in the second partition treated with 8
GWPs. This is a small basis set but with 128 trajectories allows a reasonable
searching of the main configuration space to see if the are regions where ISC
occurs.

At the end of the simulation, all three triplet states had accumulated
approximately 0.1 % population. Thus, including the nuclear motion has indeed
increased the ISC, but it is still not likely to be significant for the short-
term dynamics (over the first 0.5 ps). As a result, the triplet states will be
ignored in the final analysis.

###  III.5 Model Hamiltonian Dynamics

ML-MCTDH simulations were run for 200fs on the singlet manifold vibronic
coupling model of cyclobutanone. These started in the C2v and Cs structures
and included either all modes or just the key 12 modes:
$\nu_{1},\nu_{2},\nu_{3},\nu_{4},\nu_{5},\nu_{9},\nu_{16},\nu_{17},\nu_{19},\nu_{20},\nu_{21}$
and $\nu_{25}$. The tree structures and basis set sizes for converged
calculations are given in the SI.

Diabatic state populations after a vertical excitation to the S2 state are
shown in Fig. 4. Including all modes and starting from the planar C2v
structure, a fast step-wise decay of population is seen to S1 over the first
100 fs during which time 90 % of the population is lost from S2. Only a small
population rise is seen in S0 during this time, with only 5 % having reached
the ground-state after 200 fs. A similar picture is seen starting from the Cs
minimum energy structure, but now there is more population transfer to the
ground-state, with approximately 15 % transfer to S0 by 200 fs. Simulations
including only the key 12 modes show very similar population dynamics,
indicating that these are indeed the modes that carry the main 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 Kuhlman, Sølling, and Møller
(2013). The fact that the system retains a significant population in S1 on
this time-scale is supported by the fact that cyclobutanone is known to
fluoresce Shortridge Jr., Rusbult, and Lee (1971).

The minimum energy conical intersections (MECI) between the S2 / S1 and S1 /
S0 states from the vibronic coupling model were obtained by first minimising
the energy gap between the states of interest and then minimising the energy
along the intersection seam using a Lagrange constraint. The geometries in
terms of the normal modes are given in the SI. The S2 / S1 intersection was
found to be at 5.4 away from the FC point. Large displacements are made along
the $\nu_{5},\nu_{16},\nu_{19},\nu_{20}and\nu_{21}$ modes. The energy, at 5.97
eV is below the FC point at 6.32 eV. In the recently introduced classification
of conical intersection types Gómez, Spinlove, and Worth (2024), this is a
direct sloped intersection, which provides the step wise fast transfer
observed. The S1/ S0 intersection is further away from the FC point, at 11.5
units and an energy of 4.58 eV. It is thus energetically accessible, but the
displacements required to reach it involve large-scale motions along $\nu_{1}$
and all the B2 vibrations, hence it is not directly accessed and the transfer
to the ground-state is limited.

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

Absorption spectra calculated from the Fourier Transform of the
autocorrelation function are shown in Fig. 5. Similar results are obtained
with the full or 12D simulations. To remove artefacts 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 S1 were also calculated from
simulations started with a vertical excitation to S1.

The experimental spectrum for the S1 band is a broad, featureless band from
3.8 eV to 5 eV with the maximum at around 4.2 eV Diau, Kötting, and Zewail
(2001). This is a dark state and the intensity is low. The calculated spectrum
from the C2v structure is found to be in the correct energy region, but the
maximum is too low in energy by 0.4 eV, the band is too narrow and too
structured. Starting from the Cs structure, the structure is much reduced, but
the band is too broad and high in energy.

The experimental S2 band is again broad, but with a weak progression of peaks.
It starts at 6 eV and runs to 6.8 eV, with a maximum at 6.5 eV Diau, Kötting,
and Zewail (2001). The calculated spectrum from the C2v structure lies in the
correct energy region, perhaps 0.2 eV too low, and has too much structure. The
spectrum from the Cs structure is again too broad and high in energy.

The spectra indicate that the spectra are better reproduced by excitation from
the C2v structure even though this is a transition state. This is, however,
not a surprise due to the low barrier that means that at room temperature the
ground state wavefunction will not be purely in the ground-state Cs minimum
and will spread across the equivalent structures to be more C2v like.

Figure 5: Spectra from ML-MCTDH calculation on the singlet manifold vibronic
coupling model of cyclobutanone. (a), (c) The S1 spectrum starting from either
the (a) C2v or (c) Cs structure. (b), (d) The S2 spectrum starting from either
the (b) C2v or (d) Cs structure.

Finally, the expectation values of the normal modes can be examined to see
which modes are excited. Starting from the C2v structure, the modes
$\nu_{21},\nu_{5},\nu_{16}$ and $\nu_{19}$ all show large oscillations,
particularly $\nu_{21}$ and $\nu_{5}$. All are totally symmetric vibrations
and there is no loss of symmetry, or out-of-plane motion, evident over the 200
fs. Starting from the non-planar Cs structure, the expectation values for
modes $\nu_{1},\nu_{2}$ and $\nu_{20}$ now also undergo significant dynamics.
$\nu_{1}$ starts with a value around 2 and this decreases, showing the ring
becomes more planar over time. $\nu_{2}$, moves to a negative value, showing
the oxygen moves further out of plane.

###  III.6 Direct Dynamics

The final step to give a more complete molecular picture of the dynamics and
obtain the GUED signal, direct dynamics simulations were run 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.5. The
surfaces were calculated only at the CASSCF level 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 C2v transition state, and points were added to the database
keeping the C2v symmetry of the surface by generating appropriate symmetry
replicas of each point. An initial calculation of 100 fs was also run starting
with the C2v geometry and 8 GWPs to ensure the symmetry of the initial points
in the database.

As the sample in the experiments will be cooled below room temperature, the
main simulations were started from the Cs 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 22391 structures. A production run was then made for 500 fs, again with a
32 / 16 GWP basis set but not collecting more points. This final simulation
had population dynamics similar to the model Hamiltonian results of Fig. 4,
but with less crossing out of S2 and ending with 46 % in S1, 42 % in S2 and 12
% in S0.

The minimum energy conical intersections between the S2 / S1 and S1 / S0
surfaces provided the direct dynamics quantum chemistry database were obtained
by an optimisation procedure in the same way as those in the vibronic coupling
model. The MECI geometries are also provided in the SI both normal modes and
Cartesian coordinates are given in the SI. In both cases the MECI has moved a
bit closer to the FC point, with the S2 / S1 MECI now 3.2 units away and the
S1/S0 MECI 8.3 units away. The energies have also changed, with the S2 MECI
moving down to 5.2 eV and the S1 / S0 MECI up to 4.90 eV.

The simulation provides 512 trajectories (the centres of the GWP
configurations). Using the procedure in Sec. II.5, these can then be used to
simulate a pair distribution function (PDF), as would be produced by a GUED
experiment. The modified scattering intensity, the simulation of the original
signal is shown in the SI. The difference PDF as a function of time is shown
in Fig. 6 (b), with the reference PDF in Fig. 6(b). The reference is the pair
distribution function of the initial, Cs minimum energy structure. The peak at
1.5 Å is due to all the nearest neighbours. The peak at 2.3 Å is due to the
distances C1 \- C4, C2 \- C3 (both across the ring) as well as C2 \- O and C3
\- O. The peak at 3.2 Å is due to C4 \- O and the weak peak at 4.3 Å is due to
the distance from the protons attached to C4 to O. Atom numbering is as in
Fig. 2.

In Fig. 6, loss of intensity is observed at the neighbouring atom peak at 1.5
Å with 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
are the two new distances? 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 C2-C3 is seen to
undergo a small amplitude vibration. The across ring distance C1-C4, however,
undergoes a large amplitude vibration, stretching out to 2.5 Å which, along
with the C2 \- O and C3-O distances that extend out to 3 A Å is responsible
for the gain just below 3 Å. The gain at 3.5 Å is due to the C4-O vibration,
the loss of the original signal masked by the gain at 2.5 Å.

The overall picture then is that the molecule stretches along the C-C-O axis.
Looking at the bond lengths, it is clear that the C2-C4 distance barely
changes, and the stretch is dominated by the C-O group moving away from the
other carbon atoms. The dynamics also seems to divide into two regimes. In the
first 300 fs, the vibrations are coherent and stay in a tight grouping,
whereas at later times the distances spread and in most distances the
vibrational motion dies away. An exception is the C-O bond that seems to be
picking up amplitude. A final note on the dynamics is that the ring pucker
does not change much, 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
S1 Liu and Fang (2016) and static calculations of the S1 surface Xia  _et al._
(2015), there are no signs of ring opening over the 500 fs.

Figure 6: Pair Distribution Function (PDF) from 512 DD-vMCG trajectories
calculated on CASSCF potential surface. (a) Reference PDF for the initial
structure. (b) Difference PDF as a function of time. The blue peaks denote
loss with respect to the reference spectrum of (a), while red peaks denote
gain Figure 7: Distances between atoms in cyclobutanone from 512 DD-vMCG
trajectories calculated on CASSCF potential surface, weighted by the Gross
Gaussian Population of each trajectory. Some distances are not included due to
symmetry.

##  IV Prediction and Conclusion

We have performed quantum dynamics simulations to predict the dynamics of
cyclobutanone as it will be observed by a gas-phase ultrafast electron
diffraction experiment over 500 fs after initiating the dynamics by excitation
to the S2 electronic state. Cyclobutanone is a non-trivial molecule for
simulations. It been recently shown that the population dynamics for the
related cyclopropanone molecule are very sensitive to the electronic structure
method used Janoš and Slavíček (2023). The fact that the ground-state
equilibrium structure is a shallow well, and the molecule easily converts
between puckered structures through a planar transition state means that
choosing the initial conditions is also not straightforward as the geometry
will be temperature dependent. The potential involvement of triplet states
brings a further challenge.

The final results from the simulations indicate that triplet states are not
significantly populated over the first 500 fs. However, it is also found that
while the initial relaxation from S2 is fast, taking place in around 50 fs,
subsequent relaxation from S1 is slow. Thus at later time-scales there will be
crossing to the triplet manifold and this can result in photo-fragmentation
Chen and Ye (2003). 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 picking up energy towards the end of this time and it
is to be expected that the molecule will fragment to C3H6 \+ C-O at later
times. No ring puckering dynamics is observed. The pair distribution function
was obtained that would be measured by the electron scattering. It is seen
that this plot is not easy to interpret without the simulations to give a
molecular picture to the pair distances as each peak is due to multiple pairs,
and they overlap. Thus 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 basis set was not very large and so
this lack of convergence may mean some features are not seen. Also, the
potential surfaces used in the final simulations were only at the CASSCF
level, which may not be sufficient for a good description. And finally, the
approximations used in the dynamics, e.g. use of the LHA for integrals, and
the propagation diabatisation procedure may also bring errors. A comparison
with the experimental result will be a great test and show what is still in
need of improvement.

## Supplementary Information (SI)

Information supporting the work presented is available at DOI XXX.XXXX. This
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 are also given of datasets
containing the files from the simulations presented that are also available.

###### Acknowledgements.

This work has been funded by the EPSRC under the programme grant COSMOS
(EP/X026973/1). OB also thanks 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 ID 390715994,
from the International Max Planck Graduate School for Ultrafast imaging $\&$
Structural Dynamics (IMPRS-UFAST) and from the Christiane-Nüsslein-Vollhard-
Foundation.

## References

  * Liu  _et al._ (2020) Y. Liu, S. L. Horton, J. Yang, J. P. F. Nunes, X. Shen, T. J. A. Wolf, R. Forbes, C. Cheng, B. Moore, M. Centurion, K. Hegazy, R. Li, M.-F. Lin, A. Stolow, P. Hockett, T. Rozgonyi, P. Marquetand, X. J. Wang, and T. Weinacht, Phys. Rev. X 10, 021016 (2020).
  * Champenois  _et al._ (2021) E. G. Champenois, D. M. Sanchez, J. Yang, J. P. F. Nunes, A. Attar, M. Centurion, R. Forbes, M. Guhr, K. Hegazy, F. Ji, S. K. Saha, Y. Liu, M.-F. Lin, D. Luo, B. Moore, X. Shen, M. R. Ware, X. J. Wang, T. J. Martinez, and T. J. A. Wolf, Science 374, 178 (2021).
  * Cao, Wang, and Zhong (2021) J. Cao, X. Wang, and D. Zhong, Science 374, 34 (2021).
  * Wolf  _et al._ (2019) T. J. A. Wolf, D. M. Sanchez, J. Yang, R. M. Parrish, J. P. F. Nunes, M. Centurion, R. Coffee, J. P. Cryan, M. Gühr, K. Hegazy, A. Kirrander, R. K. Li, J. Ruddock, X. Shen, T. Vecchione, S. P. Weathersby, P. M. Weber, K. Wilkin, H. Yong, Q. Zheng, X. J. Wang, M. P. Minitti, and T. J. Martínez, Nat. Chem. 11, 504 (2019).
  * Yang  _et al._ (2018) J. Yang, X. Zhu, T. J. A. Wolf, Z. Li, J. P. F. Nunes, R. Coffee, J. P. Cryan, M. Gühr, K. Hegazy, T. F. Heinz, K. Jobe, R. Li, X. Shen, T. Veccione, S. Weathersby, K. J. Wilkin, C. Yoneda, Q. Zheng, T. J. Martinez, M. Centurion, and X. J. Wang, Science 361, 64 (2018).
  * Belluš and Ernst (1988) D. Belluš and B. Ernst, Angew. Chem., Int. Ed. Engl. 27, 797 (1988).
  * Bach and Dmitrenko (2006) R. D. Bach and O. Dmitrenko, J. Am. Chem. Soc. 128, 4598 (2006).
  * Trost (1986) B. M. Trost, in  _Small Ring Compounds in Organic Synthesis I_ , edited by A. de Meijere (Springer Berlin Heidelberg, Berlin, Heidelberg, 1986) pp. 3–82.
  * Lumbroso  _et al._ (2015) A. Lumbroso, S. Catak, S. Sulzer-Mossé, and A. De Mesmaeker, Tetrahedron Lett. 56, 2397 (2015).
  * Izquierdo  _et al._ (2005) S. Izquierdo, F. Rúa, A. Sbai, T. Parella, A. Álvarez-Larena, V. Branchadell, and R. Ortuño, J. Org. Chem. 70, 7963 (2005).
  * Lee-Ruff and Wells (2008) E. Lee-Ruff and D. Wells, Nucleosides Nucleotides Nucleic Acids 27, 484 (2008).
  * Mihovilovic  _et al._ (2005) M. Mihovilovic, P. Kapitan, J. Rydz, F. Rudroff, F. Ogink, and M. Fraaije, J. Mol. Catal., B Enzym. 32, 135 (2005).
  * Dömling and Ugi (2000) A. Dömling and I. Ugi, Angew. Chem. Int. Ed. 39, 3168 (2000).
  * Pirrung and Wang (2009) M. C. Pirrung and J. Wang, J. Org. Chem. 74, 2958 (2009).
  * Kirscher (1905) N. Kirscher, Zh. Russ. Fiz.-Khim. O-va. 37, 106 (1905).
  * Salaun, Champion, and Conia (1988) J. Salaun, J. Champion, and J. Conia, Org. Synth. 50-9, 320 (1988).
  * Seebach and Beck (1988) D. Seebach and A. Beck, Org. Synth. 50-9, 316 (1988).
  * Benson and Kistiakowsky (1942) S. W. Benson and G. B. Kistiakowsky, J. Am. Chem. Soc. 64, 80 (1942).
  * Moule (1976) D. C. Moule, J. Chem. Phys. 64, 3161 (1976).
  * Hemminger, Carless, and Lee (1973) J. C. Hemminger, H. A. J. Carless, and E. K. C. Lee, J. Am. Chem. Soc. 95, 682 (1973).
  * Diau, Kötting, and Zewail (2001) E. W.-G. Diau, C. Kötting, and A. H. Zewail, ChemPhysChem 2, 294 (2001).
  * Campbell, Schlag, and Ristow (1967) R. J. Campbell, E. W. Schlag, and B. W. Ristow, J. Am. Chem. Soc. 89, 5098 (1967).
  * Tang and Lee (1976) K. Y. Tang and E. K. C. Lee, J. Phys. Chem. 80, 1833 (1976).
  * Hemminger and Lee (1972) J. C. Hemminger and E. K. C. Lee, J. Chem. Phys. 56, 5284 (1972).
  * Kao  _et al._ (2020) M.-H. Kao, R. K. Venkatraman, M. N. R. Ashfold, and A. J. Orr-Ewing, Chem. Sci. 11, 1991 (2020).
  * Xia  _et al._ (2015) S.-H. Xia, X.-Y. Liu, Q. Fang, and G. Cui, J. Phys. Chem. A 119, 3569 (2015).
  * Liu and Fang (2016) L. Liu and W.-H. Fang, J. Chem. Phys. 144 (2016).
  * Kuhlman, Sølling, and Møller (2012) T. S. Kuhlman, T. I. Sølling, and K. B. Møller, ChemPhysChem 13, 820 (2012).
  * Kuhlman  _et al._ (2012) T. S. Kuhlman, S. P. A. Sauer, T. I. Sølling, and K. B. Møller, J. Chem. Phys. 137, 22A522 (2012).
  * Kuhlman, Sølling, and Møller (2013) T. S. Kuhlman, T. I. Sølling, and K. B. Møller,  _EPJ Web of Conferences 41, 02033 (2013) DOI: 10.1051/epjconf/20134102033 ©C Owned by the authors, published by EDP Sciences, 2013_ , EPJ Web of Conferences 41, 02033 (2013).
  * Werner  _et al._ (2022) H.-J. Werner, P. J. Knowles, P. Celani, W. Györffy, A. Hesselmann, D. Kats, G. Knizia, A. Köhn, T. Korona, D. Kreplin, R. Lindh, Q. Ma, F. R. Manby, A. Mitrushenkov, G. Rauhut, M. Schütz, K. R. Shamasundar, T. B. Adler, R. D. Amos, S. J. Bennie, A. Bernhardsson, A. Berning, J. A. Black, P. J. Bygrave, R. Cimiraglia, D. L. Cooper, D. Coughtrie, M. J. O. Deegan, A. J. Dobbyn, K. Doll, M. Dornbach, F. Eckert, S. Erfort, E. Goll, C. Hampel, G. Hetzer, J. G. Hill, M. Hodges, T. Hrenar, G. Jansen, C. Köppl, C. Kollmar, S. J. R. Lee, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, B. Mussard, S. J. McNicholas, W. Meyer, T. F. Miller III, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. A. Peterson, K. Pflüger, R. Pitzer, I. Polyak, M. Reiher, J. O. Richardson, J. B. Robinson, B. Schröder, M. Schwilk, T. Shiozaki, M. Sibaev, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, J. Toulouse, M. Wang, M. Welborn, and B. Ziegler, “Molpro, version 2022.1 , a package of ab initio programs,” (2022), see https://www.molpro.net.
  * Werner  _et al._ (2012) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, WIREs Comput. Mol. Sci. 2, 242 (2012).
  * Werner  _et al._ (2020) H.-J. Werner, P. J. Knowles, F. R. Manby, J. A. Black, K. Doll, A. Heßelmann, D. Kats, A. Köhn, T. Korona, D. A. Kreplin, Q. Ma, T. F. Miller, A. Mitrushchenkov, K. A. Peterson, I. Polyak, G. Rauhut, and M. Sibaev, J. Chem. Phys. 152 (2020).
  * Worth  _et al._ (2023) G. A. Worth, K. Giri, G. W. Richings, I. Burghardt, M. H. Beck, A. Jäckle, and H.-D. Meyer, “Quantics package, version 2.0 development. see http://www.chem.ucl.ac.uk/quantics,” (2023).
  * Worth (2020) G. A. Worth, Comp. Phys. Comm. 248, 107040 (2020).
  * Cattarius, Markmann, and Worth (2007) C. Cattarius, A. Markmann, and G. A. Worth, “The VCHAM program, see http://www.pci.uni-heidelberg.de/tc/usr/mctdh/,” (2007).
  * Köppel, Domcke, and Cederbaum (1984) H. Köppel, W. Domcke, and L. S. Cederbaum, Adv. Chem. Phys. 57, 59 (1984).
  * Worth and Cederbaum (2004) G. A. Worth and L. S. Cederbaum, Ann. Rev. Phys. Chem. 55, 127 (2004).
  * Beck  _et al._ (2000) M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000).
  * Meyer, Gatti, and Worth (2009) H.-D. Meyer, F. Gatti, and G. A. Worth,  _Multidimens. Quantum Dyn. MCTDH Theory Appl._ (Wiley-VCH, Weinheim, Germany, 2009).
  * Wang and Thoss (2003) H. Wang and M. Thoss, J. Chem. Phys. 119, 1289 (2003).
  * Manthe (2008) U. Manthe, J. Chem. Phys. 128, 164116 (2008).
  * Vendrell and Meyer (2011) O. Vendrell and H.-D. Meyer, J. Chem. Phys. 134, 44135 (2011).
  * Mendive-Tapia  _et al._ (2017) D. Mendive-Tapia, T. Firmino, H.-D. Meyer, and F. Gatti, Chem. Phys. 482, 113 (2017).
  * Burghardt, Meyer, and Cederbaum (1999) I. Burghardt, H.-D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 111, 2927 (1999).
  * Richings  _et al._ (2015) G. W. Richings, I. Polyak, K. E. Spinlove, G. A. Worth, I. Burghardt, and B. Lasorne, Int. Rev. Phys. Chem. 34, 269 (2015).
  * Heller (1975) E. J. Heller, J. Chem. Phys. 62, 1544 (1975).
  * Christopoulou, Freibert, and Worth (2021) G. Christopoulou, A. Freibert, and G. A. Worth, J. Chem. Phys. 154, 124127 (2021).
  * Centurion, Wolf, and Yang (2022) M. Centurion, T. J. Wolf, and J. Yang, Ann. Rev. Phys. Chem. 73, 21 (2022).
  * Wolf and Martinez (2021) T. J. Wolf and T. J. Martinez, “Diffraction_simulation, see https://github.com/thomasjawolf/diffraction_simulation,” (2021).
  * O’Toole  _et al._ (1991) L. O’Toole, P. Brint, C. Kosmidis, G. Boulakis, and P. Tsekeris, J. Chem. Soc., Faraday Trans. 87, 3343 (1991).
  * Cattarius  _et al._ (2001) C. Cattarius, G. A. Worth, H.-D. Meyer, and L. S. Cederbaum, J. Chem. Phys. 115, 2088 (2001).
  * Shortridge Jr., Rusbult, and Lee (1971) R. G. Shortridge Jr., C. F. Rusbult, and E. K. C. Lee, J. Am. Chem. Soc. 93, 1863 (1971).
  * Gómez, Spinlove, and Worth (2024) S. Gómez, K. E. Spinlove, and G. A. Worth, Phys. Chem. Chem. Phys. 26, 1829 (2024).
  * Janoš and Slavíček (2023) J. Janoš and P. Slavíček, J. Chem.Theory Comput. 19, 8273 (2023).
  * Chen and Ye (2003) Y. Chen and S. Ye, Int. J. Quantum Chem. 97, 725 (2003).