# Predicting the photodynamics of cyclobutanone triggered by a laser pulse at
200 nm and its MeV-UED signals – a trajectory surface hopping and XMS-CASPT2
perspective

Jiří Janoš  Department of Physical Chemistry, University of Chemistry and
Technology, Technická 5, Prague 6, 166 28, Czech Republic  Centre for
Computational Chemistry, School of Chemistry, University of Bristol, Bristol
BS8 1TS, United Kingdom  Joao Pedro Figueira Nunes  Diamond Light Source Ltd,
Didcot, UK.  Daniel Hollas  Centre for Computational Chemistry, School of
Chemistry, University of Bristol, Bristol BS8 1TS, United Kingdom  Petr
Slavíček<EMAIL_ADDRESS>Department of Physical Chemistry, University
of Chemistry and Technology, Technická 5, Prague 6, 166 28, Czech Republic
Basile F. E. Curchod<EMAIL_ADDRESS>Centre for Computational
Chemistry, School of Chemistry, University of Bristol, Bristol BS8 1TS, United
Kingdom

###### Abstract

This work is part of a prediction challenge that invited
theoretical/computational chemists to predict the photochemistry of
cyclobutanone in the gas phase, excited at 200 nm by a laser pulse, and the
expected signal that will be recorded during a time-resolved megaelectronvolt
ultrafast electron diffraction (MeV-UED). We present here our theoretical
predictions based on a combination of trajectory surface hopping with XMS-
CASPT2 (for the nonadiabatic molecular dynamics) and Born–Oppenheimer
molecular dynamics with MP2 (for the athermal ground-state dynamics following
internal conversion), coined (NA+BO)MD. The initial conditions were sampled
from Born–Oppenheimer molecular dynamics coupled to a quantum thermostat. Our
simulations indicate that the main photoproducts after 2 ps of dynamics are CO
+ cyclopropane (50%), CO + propene (10%), and ethene and ketene (34%). The
photoexcited cyclobutanone in its second excited electronic state S2 can
follow two pathways for its nonradiative decay: (i) a ring-opening in S2 and a
subsequent rapid decay to the ground electronic state, where the photoproducts
are formed, or (ii) a transfer through a closed-ring conical intersection to
S1, where cyclobutanone ring opens and then funnels to the ground state.
Lifetimes for the photoproduct and electronic populations were determined. We
calculated a stationary MeV-UED signal [difference pair distribution function
– $\Delta$PDF$(r)$] for each (interpolated) pathway as well as a time-resolved
signal [$\Delta$PDF$(r,t)$ and $\Delta I/I(s,t)$] for the full swarm of
(NA+BO)MD trajectories. Furthermore, our analysis provides time-independent
basis functions that can be used to fit the time-dependent experimental UED
signals [both $\Delta$PDF$(r,t)$ and $\Delta I/I(s,t)$] and potentially
recover the population of photoproducts. We also offer a detailed analysis of
the limitations of our model and their potential impact on the predicted
experimental signals.

##  I Introduction

The field of nonadiabatic molecular dynamics simulations, concerned with the
dynamics of molecular systems beyond the Born–Oppenheimer approximation, has
experienced a growing interest by the physical chemistry and chemical physics
community over the past three decades.Tully (2012); Persico and Granucci
(2014); Agostini and Curchod (2019); Mai and González (2020); González and
Lindh (2020) This interest was greatly stimulated by the development of
frameworks for on-the-fly nonadiabatic dynamics around the 90s,Tully (1990,
1998); Martínez, Ben-Nun, and Levine (1996); Martínez and Levine (1997);
Bittner and Rossky (1995); Kapral and Ciccotti (1999); Donoso and Martens
(1998); Coker and Xiao (1995) combined with advances in efficient electronic-
structure methods for excited states,Schirmer (1982); Finley  _et al._ (1998);
Koch  _et al._ (1995); Doltsinis and Marx (2002); Casida (1995); Petersilka,
Gossmann, and Gross (1996); González, Escudero, and Serrano-Andrés (2012);
Matsika (2021) which allowed for dynamics simulations of photophysical and
photochemical processes unraveled by the emergence of femtochemistry.Zewail
(2000) The field of nonadiabatic dynamics has continued its development over
the years and reached a level of maturity that permits its use to interpret
complex spectroscopic experiments or rationalize photophysical and
photochemical events.

Nonadiabatic molecular dynamics remains a challenging technique to operate:
the overall nuclear dynamics is sensitive to the photoexcitation processSuchan
_et al._ (2018); Pieroni  _et al._ (2023); Prlj, Hollas, and Curchod (2023)
and is out-of-equilibrium in nature, accounting for nuclear quantum effects is
often required,Crespo-Otero and Barbatti (2018); Curchod and Martínez (2018);
Persico and Granucci (2014); Subotnik  _et al._ (2016); Hammes-Schiffer
(2022); Ryabinkin, Joubert-Doriol, and Izmaylov (2017); Ghosh  _et al._
(2020); Agostini, Tavernelli, and Ciccotti (2018) and including the effect of
an environment challenges most frequently used methods for this
purpose.Santoro  _et al._ (2021); Bondanza  _et al._ (2022); Tavernelli,
Curchod, and Rothlisberger (2011) Perhaps most importantly, the nonadiabatic
dynamics can be very sensitive to the electronic-structure method employed to
describe the different potential energy surfaces and couplings between
them.Janoš and Slavíček (2023) While photophysical processes, where a
photoexcited molecule explores a limited portion of the nuclear configuration
space, are often less problematic to simulate for larger molecular systems
than photochemical processes, the formation of photoproducts stresses more the
approximations underlying nonadiabatic molecular dynamics simulations. Hence,
it appears fair to say that, in 2024, nonadiabatic molecular dynamics
simulations can be performed to study the photochemistry of small- to medium-
size molecules, in vacuo preferably, but at the cost of extensive tests prior
to any production runs.

But even in the best-case scenario described above, can nonadiabatic molecular
dynamics be predictive at all? This is the central question of this Special
Topic, to which this article contributes. The goal of this series of works is
to predict the outcome of an ultrafast pump-probe gas-phase experiment before
it is conducted. Experimentally, a molecule will be photoexcited by a short
pump laser pulse at 200 nm and subsequently probed at different time delays by
a probe, which in this case is a bundle of MeV electrons that will diffract on
the molecule, revealing its geometrical structure. The precise details of the
planned experiments are given in Sec. II.5. The molecule chosen for this
challenge is cyclobutanone, a cyclic ketone that can undergo various
photochemical reactions depending on the excitation wavelength. The current
knowledge of cyclobutanone photochemistry triggered by excitation in either
the S1 or S2 electronic state is described in the following paragraphs and
schematically summarized in Fig. 1.

Figure 1: Schematic representation of the photochemistry of cyclobutanone
triggered by photoexcitation to S1($n\pi^{*}$) or S2($3s$) as described by
literature (see main text). Red symbols indicate the pathways discussed and
questions addressed in the present work.

The (photo)reactivity of cyclobutanone upon photoexcitation in its first
excited electronic state (S1) is governed by the ring strain in the four-
member ring as for the other cyclic ketones.Xia  _et al._ (2015); Kao  _et
al._ (2020) The smallest analog of the cyclic-ketone family, the highly
strained cyclopropanone, undergoes an ultrafast deactivation upon excitation
in its lowest excited electronic state S1, releasing CO within a hundred
femtoseconds.Janoš and Slavíček (2023); Cui and Fang (2011) The larger and
less strained analogues of this series, cyclopentanone and cyclohexanone,
exhibit a slower deactivation involving intersystem crossings to the lowest
triplet electronic state T1 (with a $n\pi^{*}$ character).Xia  _et al._
(2015); Kao  _et al._ (2020) Cyclobutanone lies somewhere at the interface
thanks to its medium ring strain, with internal conversion (IC) and
intersystem crossing (ISC) processes competing when the molecule is excited in
its S1 electronic state. The reported photoproducts of cyclobutanone upon
excitation in S1 are ethene and ketene (denoted as C2 products), or CO and
cyclopropane/propene (coined C3 products), and both families of photoproducts
were observed in the gas phase and in solution.Denschlag and Lee (1967); Turro
and Southam (1967) The ratio of photoproducts was found to vary between a
photoexcitation wavelength of 340 nm (C3:C2 ratio 7:1) and 315 nm (C3:C2 ratio
1:2), but it remains constant below 315 nm.Campbell, Schlag, and Ristow
(1967); Tang and Lee (1976) This stark change in photoproduct quantum yields
was explained by an activation barrier to ring-opening in the first excited
state with a value predicted between 2 to 7 kcal/mol.Xia  _et al._ (2015);
Diau, Kötting, and Zewail (2001) When the wavelengths used produce
cyclobutanone with internal energy below the barrier (340-315 nm), the direct
dissociation is hampered, and ISC from S1 to triplet states becomes
competitive (the ISC proceeds on a nanosecond timescaleDiau, Kötting, and
Zewail (2001)). When photoexcited above the barrier ($<315$ nm), cyclobutanone
undergoes a Norrish type-I cleavage leading to the formation of both the C3
and C2 products. Time-dependent experiments in solution have revealed that the
dissociation happens on a subpicosecond timescale (0.65 ps),Kao  _et al._
(2020) although previous experiments suggested a slightly longer lifetime of 5
picoseconds with 307 nm pump.Diau, Kötting, and Zewail (2001) Theoretical
calculations using the ab initio multiple spawning (AIMS) method combined with
SA2-CAS(12,11)/6-31G* method predicted the S1 lifetime of $\approx 450$ fs and
the ratio of C3:C2 products 1:1.Liu and Fang (2016)

The photodynamics of cyclobutanone resulting from an excitation into its
second excited electronic state (S2, exhibiting a Rydberg 3s character in the
Franck-Condon region) has been studied both experimentallyTrentelman  _et al._
(1990); Kuhlman, Sølling, and Møller (2012); Kuhlman  _et al._ (2013) and
theoretically.Kuhlman  _et al._ (2012) Kuhlman and coworkers performed time-
resolved mass spectroscopy (TRMS) and time-resolved photoelectron spectroscopy
(TRPES) experiments on cyclobutanone, triggered by a 200 nm laser pumpKuhlman,
Sølling, and Møller (2012) and interpreted them using multi-configuration
time-dependent Hartree (MCTDH) quantum-dynamics simulations.Kuhlman  _et al._
(2012) The time dependence of the parent ion signal from TRMS was fitted to a
sequential biexponential decay with two-time constants $\tau_{1}=0.08\pm 0.01$
ps and $\tau_{2}=0.74\pm 0.01$ ps. The authors assumed that the longer
lifetime, $\tau_{2}$, must reflect internal conversion to the S1($n\pi^{*}$)
state because the S2 electronic state is of Rydberg character – such an
electronic-state character leads to generally bound states and a direct
dissociation or isomerization of the 3s state appeared unlikely to the
authors. An experimental signal corresponding to a fragment with 42 $m/z$
(cyclopropane, propylene, or ketene) yielded similar lifetimes. Based on the
correspondence between the parent and fragment lifetimes, the authors stated
that the processes following IC to the S1($n\pi^{*}$) should be ultrafast. The
shorter lifetime, $\tau_{2}$, was attributed to the formation of an
intermediate structure in the S2 state, probably the planarization of the
initially puckered ground-state geometry. The TRPES experiments corroborated
the TRMS results with similar lifetimes $\tau_{1}=0.31\pm 0.06$ ps and
$\tau_{2}=0.74\pm 0.02$ ps. Furthermore, oscillations in the TRPES signal were
recorded, exhibiting a frequency of 35 cm-1. The MCTDH simulations supporting
this experiment used a 5-dimensional vibronic coupling Hamiltonian fitted to
EOM- and LR-CCSD energies.Kuhlman  _et al._ (2012) The theoretical model
includes the following modes: the carbonyl stretch to describe the changing
C=O bond length, the ring-puckering and carbonyl out-of-plane deformation to
describe the changing structure of the ring, and the angle between the C=O
bond and the C-C-C plane combined with ring modes affecting the length of the
C-C bonds. The calculated decay of the S2-photoexcited cyclobutanone to the S1
electronic state was also fitted by a biexponential function, leading to time
constants of 0.95 ps and 6.32 ps. The lifetimes observed for the IC process
were attributed to (i) an immediate activation of coupling modes leading to a
fast population transfer from S2 towards S1 (direct mechanism), and (ii) an
internal vibrational energy redistribution followed by a slower population
transfer to S1 (indirect mechanism). We note that no population reached the
ground electronic state S0 or the S3(3p) electronic state in this dynamics.
The trapping in the S1 state is not surprising as these simulations are
restricted by the underlying harmonic approximation of the linear vibronic
coupling model and cannot describe direct bond dissociation. This restriction
also hampered a possible dissociation in the S2. Another experiment by
Trentelman and coworkersTrentelman  _et al._ (1990) measured rotational
spectra of CO fragments following photoexcitation of cyclobutanone at 193 nm,
revealing two groups of fragments with different rotational temperatures. The
high-energy fragments ($\approx 3000$ K) were attributed to direct
dissociation into CO and cyclopropane. The low-temperature fragments ($\approx
200$ K) were attributed to the formation of an excited ketene, which then
dissociates into CH2 and CO. The ratio of these channels was found to be
85:15. The authors estimated the ratio of C3:C2 products to be 57:43 based on
the assumption that cyclobutanone in S2 undergoes an ultrafast IC to the S1
state, which then drives the dissociation, and extrapolation of the quantum
yields in the S1 state.

Moving to the current prediction challenge, the time-resolved MeV-UED
experiment on cyclobutanone in vacuo will use a 200 nm pump pulse to promote
the molecule in its S2($3s$) state. Time-resolved MeV-UED experiments are
sensitive to the geometry of the molecule probed at a given time following
photoexcitationYang  _et al._ (2018); Wolf  _et al._ (2019); Centurion, Wolf,
and Yang (2022) (electronic information can be obtained too in some casesYang
_et al._ (2020)) and therefore offers a perfect tool to challenge the nuclear
wave packets produced by nonadiabatic molecular dynamics simulations and
unravel the actual mechanistic details of the photodissociation that have
eluded earlier experiments (question marks in Fig. 1). In the following, we
propose an in silico version of this experiment, where we simulate the
excited-state dynamics of cyclobutanone and the subsequent dynamics of its
photoproducts in the ground electronic state from the perspective of
nonadiabatic molecular dynamics. We also predict the MeV-UED signal stemming
from our simulations. The time constraints imposed by this prediction
challenge mean that we had to opt for a compromise between computational
efficiency and accuracy – XMS-CASPT2 for the electronic structure and
trajectory surface hopping (TSH) for the nonadiabatic molecular dynamics – and
that the number of tests we would usually conduct to validate the electronic-
structure methods had to be kept to a minimum. We therefore supplement our
results section with a discussion on the potential limitations of our model
(Sec. IV) and their possible influence on the results of our simulations.

##  II Methods

###  II.1 Electronic-structure methods and benchmark

The electronic structure of cyclobutanone in this work was mainly described by
the extended multistate complete active space second-order perturbation theory
(XMS-CASPT2)Park and Shiozaki (2017) based on a state-averaged complete active
space self-consistent field (SA-CASSCF) reference wavefunction as implemented
in the BAGEL package.Shiozaki (2018) The state-averaging procedure considered
three singlet states with equal weights. The active space comprised 8
electrons in 8 orbitals, namely two $\sigma_{\mathrm{CC}}$ adjacent to
carbonyl, $\pi_{\mathrm{CO}}$, $n_{\mathrm{O}}$, $\pi_{\mathrm{CO}}^{*}$,
$Ry(3s)$, and two corresponding $\sigma^{*}_{\mathrm{CC}}$ orbitals (see Fig.
S1 for a representation of the active-space orbitals). All calculations used
the aug-cc-pVDZ basis set combined with cc-pVTZ-jkfit auxiliary basis for the
density fitting. This overall level of electronic-structure theory is denoted
XMS-CASPT2(8/8) in the following and was used to optimize the electronic-state
minima, transition states (TSs), and minimum energy conical intersections
(CIs) of cyclobutanone, see the supplementary material (SM) for the
corresponding Cartesian coordinates. Linear interpolations in internal
coordinates (LIICs) between these critical geometries were used to investigate
the potential energy surfaces and further benchmark the level of electronic-
structure theory. In particular, the XMS-CASPT2(8/8) strategy described above
was thoroughly compared (excitation energies and LIICs) to equation-of-motion
coupled-cluster singles doubles (EOM-CCSD) and linear-response time-dependent
density functional theory (LR-TDDFT) employing the Tamm-Dancoff approximation
(TDA) with the B97XD exchange-correlation functional. Both the EOM-CCSD and
LR-TDDFT/TDA/B97XD methods were performed in Gaussian 09, Revision D.01.Frisch
_et al._ (2013) The excitation energies were also benchmarked against the EOM-
CC3 methodPaul, Myhre, and Koch (2021) implemented in the eT code.Folkestad
_et al._ (2020) We note that the presence of a carbonyl group in the studied
molecule and the importance of a low-lying electronic state with a $n\pi*$
character hampers the use of ADC(2).Marsili, Prlj, and Curchod (2021) See the
SM for additional details of this comparison not discussed in the main text.

Møller–Plesset second-order perturbation theory (MP2) as implemented in
Gaussian 09, Revision D.01Frisch  _et al._ (2013) was also used to model the
ground state and sample ground-state nuclear density of cyclobutanone and its
photoproducts. For optimization of the photoproducts (CO, cyclopropane,
propene, ketene, and ethene) and also ground-state continuation of the
nonadiabatic molecular dynamics simulations (see below), we used the MP2/aug-
cc-pVDZ level of theory. This strategy was benchmarked against XMS-CASPT2 and
unrestricted PBE0 (UPBE0) calculations (see SM). The UPBE0 calculations were
done in Gaussian 09, Revision D.01.Frisch  _et al._ (2013) with a mixed guess
for the initial wave function. For sampling the ground-state nuclear density
with molecular dynamics, we employed the MP2 method with a smaller basis set
cc-pVDZ.

Spin-orbit coupling (SOC) interactions were modeled at the SA-CASSCF(8/8)
level of theory, state-averaging over three singlet and three triplet states
with equal weights [SA(3S,3T)-CASSCF(8/8)] using the same active space as
described above and the Molpro2012 package.Werner  _et al._ (2012) The
energies of the triplet states (and their character) obtained at this level of
theory were benchmarked against the EOM-CCSD and LR-TDDFT/TDA/B97XD methods
(see SM). The SOC matrix elements were evaluated using the Breit–Pauli
Hamiltonian $\hat{H}_{\mathrm{SOC}}$ and the values reported here correspond
to the magnitude of the SOC between a given singlet electronic state Sn and
the MS sublevels of a triplet electronic state Tm:

$H_{\mathrm{SOC}}^{\mathrm{T_{m},S_{n}}}=\sqrt{\sum_{\mathrm{M_{S}}=-1}^{1}\left|\left\langle\psi_{\mathrm{T_{m},M_{S}}}\middle|\hat{H}_{\mathrm{SOC}}\middle|\psi_{\mathrm{S_{n},0}}\right\rangle\right|^{2}}$

###  II.2 Photoabsorption cross-section

Two approaches were used to simulate the photoabsorption cross-section of
cyclobutanone: (i) Franck-Condon Herzberg-Teller (FCHT) method to obtain an
accurate depiction of the vibronic structure of the
$\mathrm{S}_{2}\leftarrow\mathrm{S}_{0}$ band and (ii) the nuclear ensemble
approach (NEA) to reproduce the shape of this band and sample initial
conditions for the nonadiabatic molecular dynamics simulations.Crespo-Otero
and Barbatti (2012)

We calculated the vibrationally-resolved absorption spectrum of cyclobutanone
using the FCHT terms at the LR-TDDFT/TDA/B97XD/aug-cc-pVDZ level of theory, as
implemented in Gaussian 09, Revision D.01.Frisch  _et al._ (2013) This
calculation employed a harmonic approximation around the optimized minima of
the ground electronic state S0 and the second excited electronic state S2, and
considered a temperature of 298 K. We note that the benchmark described above
and presented in the SM indicates that LR-TDDFT/TDA/B97XD/aug-cc-pVDZ offers
an adequate description of S2 minimum.

The NEA was used to capture the envelope of the
$\mathrm{S}_{2}\leftarrow\mathrm{S}_{0}$ band of the photoabsorption cross-
section.Crespo-Otero and Barbatti (2012) The ground-state nuclear probability
density of cyclobutanone was represented by 5000 geometries that were sampled
from ab initio Born–Oppenheimer molecular dynamics (BOMD) at the MP2/cc-pVDZ
with our molecular dynamics code ABIN.Hollas  _et al._ (2019) Both thermal
effects (298 K) and zero-point energy effects were included by the quantum
thermostat (QT) based on the generalized Langevin equation.Ceriotti, Bussi,
and Parrinello (2009, 2010) The cyclobutanone molecule was first equilibrated
and then propagated for 96 ps with a time step of 20 a.u. ($\sim$ 0.5 fs).
Recent studies suggest that QT-BOMD offers a more reliable strategy to
approximate the ground-state nuclear probability of flexible molecules.Prlj
_et al._ (2022); Prlj, Hollas, and Curchod (2023) The photoabsorption cross-
section spectra were obtained using LR-TDDFT/TDA/B97XD/aug-cc-pVDZ excitation
energies ($\Delta E_{\mathrm{S}_{0}\leftarrow\mathrm{S}_{2}}$) and transition
dipole moments ($\boldsymbol{\mu}_{\mathrm{S}_{0}\leftarrow\mathrm{S}_{2}}$)
between S0 and S2 on the support of the 5000 geometries extracted from each
distribution, using

$\sigma_{\mathrm{S}_{2}\leftarrow\mathrm{S}_{0}}(E)=\frac{\pi}{3\hbar\varepsilon_{0}c}\frac{1}{N_{s}h\sqrt{2\pi}}\sum_{i=1}^{N_{s}}\Delta
E_{\mathrm{S}_{0}\leftarrow\mathrm{S}_{2},i}|\boldsymbol{\mu}_{\mathrm{S}_{0}\leftarrow\mathrm{S}_{2},i}|^{2}\exp{\left(-\frac{\left(E-\Delta
E_{\mathrm{S}_{0}\leftarrow\mathrm{S}_{2},i}\right)^{2}}{2h^{2}}\right)},$

where $N_{s}$ is the number of sampled geometries (5000 in the present case)
and $h$ is the the Gaussian broadening factor, which was set according to
Silverman’s rule of thumb

$h\approx\left(\frac{4s^{5}}{3N_{s}}\right)^{\frac{1}{5}}$

where $s$ is the standard deviation of the excitation energies.Sršeň, Hollas,
and Slavíček (2018); Silverman (1986)

###  II.3 Nonadiabatic molecular dynamics

The photodynamics of cyclobutanone upon photoexcitation at 200 nm was
simulated with two versions of the trajectory surface hopping (TSH) method
implemented in the ABIN codeHollas  _et al._ (2019): the fewest-switches
surface hopping (FSSH)Tully (1990) and Landau-Zener surface hopping
(LZSH).Suchan, Janoš, and Slavíček (2020) The initial conditions (nuclear
positions and momenta) were selected from the 5000 geometries sampled from the
MP2/cc-pVDZ QT-BOMD simulations described above. An energy filterMartínez-Mesa
and Saalfrank (2015) was applied to the 5000 sampled geometries: only the
geometries with a S2$\leftarrow$S0 electronic transition falling within an
energy window between 6.1764 and 6.2220 eV – 200.74 and 199.27 nm respectively
– were selected. The applied energy window reflects the energy width of the
pulse that will be used in the experiment (see below). The 119 initial
conditions selected by the energy windowing were then vertically promoted to
the S2 electronic state, where the trajectories were initiated. For the FSSH
dynamics, the energy-based decoherence scheme was applied using the
recommended 0.1 a.u. value.Granucci and Persico (2007) The nonadiabatic
transitions were mediated by the nonadiabatic coupling vectors (the latter
being calculated only when the energy difference between coupled electronic
states dropped below 1.5 eV). The nuclear velocities were rescaled along the
nonadiabatic coupling vector after a successful hop, and the time step was set
to 10 a.u. ($\sim$0.25 fs). The total energy conservation was carefully
monitored for all trajectories. The active state chosen for the XMS-
CASPT2(8/8) method was stable until the molecule opened its ring, cleaving the
C-C bond adjacent to the carbonyl group. When this process happened, the
Rydberg $3s$ state was pushed higher in energy, forcing the Rydberg orbital to
rotate out of the active space slowly. Nevertheless, the orbital rotations
causing jumps in total energy occurred mainly when the molecule was already
evolving in the ground electronic state. This means that total energy was well
conserved in the segment of the dynamics evolving in S2, but suffered moderate
jumps of tenths of eV in total energy when the trajectories were evolving in
S0 and, in some cases, in S1. We note that a larger active space,
XMS(4)-CASPT2(8/9), would be more robust and stable but computationally
intractable (in particular, given the time constraints of this prediction
challenge).

To collect 2 ps of simulation time for all trajectories, we switched from the
time-consuming FSSH/XMS-CASPT2 level of theory to ab initio BOMD with MP2/aug-
cc-pVDZ (time step of 10 a.u) once the molecule reached the ground electronic
state and left the nonadiabatic region. The criteria for switching from
FSSH/XMS-CASPT2 to BOMD/MP2 were i) a formation of stable C2 or C3 products
and ii) a distance between the nearest carbon atoms of the two dissociated
fragments got bigger than 8 Å. When a ground-state cyclobutanone was formed,
the FSSH/XMS-CASPT2 simulations were continued. We coin this strategy of
combining nonadiabatic and Born–Oppenheimer molecular dynamics (NA+BO)MD in
the following. We note that the strategy of switching between nonadiabatic
molecular dynamics to BOMD was successfully used (and benchmarked) in earlier
studies.Mignolet, Curchod, and Martínez (2016); Pathak  _et al._ (2020) In 8
cases, SA-CASSCF failed to converge before the dissociated products departed
far enough to fulfil criterion ii). In these cases, we switched to a BOMD at
the geometry where the SA-CASSCF convergence issue was experienced. We stress
that, for all trajectories, stable ground-state C2 or C3 products were always
formed by the time the switch to BOMD was operated.

###  II.4 Analysis of the nonadiabatic molecular dynamics simulations

The time-dependent electronic-state populations were evaluated as an average
of the FSSH active electronic states over all trajectories, i.e. the states on
which the molecule evolves. We estimated the error bars of the $i$-th
electronic-state population at a given time $t$ using a multinomial
distribution as
$\epsilon_{i}(t)=z\sqrt{\frac{1}{N}p_{i}(t)\left(1-p_{i}(t)\right)}$, where
$N$ is the number of trajectories, $p_{i}(t)$ is the population of $i$-th
electronic state at time $t$, and $z=1.96$ corresponds to a 95 % confidence
interval. To account for the temporal smearing of the results, the calculated
populations $p_{i}(t)$ were then convoluted with a normalized Gaussian
envelope $\chi_{\mathrm{pul.env.}}$ of the excitation pulse (described in the
experimental details) based on

$\tilde{p}_{i}(t)=\int_{-\infty}^{\infty}{p}_{i}(t-\tau)\chi_{\mathrm{pul.env.}}(\tau)\mathrm{d}\tau$

To fit the population traces and obtain the lifetimes, we considered a
consecutive reaction scheme
$\mathrm{S}_{2}\rightarrow\mathrm{S}_{1}\rightarrow\mathrm{S}_{0}$. For more
details on the fitting and unconvoluted population traces, see SM.

To calculate the time-dependent relative populations of the different
photoproducts, we developed a classification scheme that sorts molecular
geometries extracted from the trajectories into six possible photoproduct
categories (see Fig. 2). We distinguish the following (photo)products: CO +
cyclopropane, CO + propene, CO + ^.CH2-CH2-CH2^. biradical, ethene + ketene,
any transient/intermediate opened structures, and cyclobutanone, which is both
initial structure and one of the photoproducts when ring-closing happens. The
classification scheme used a series of geometrical parameters (see Fig. S5 in
the SM for the parameters and categories) and was validated by visual analysis
of the trajectories. Its sensitivity to the choice of parameters was tested
and this scheme could capture all the types of structures encountered during
the nonadiabatic and ground-state dynamics. The quantum yield for a given
(photo)product $k$, $\phi_{\mathrm{k}}(t)$, was evaluated as
$\phi_{\mathrm{k}}(t)=\frac{N_{k}(t)}{N}$, where $N$ is the total number of
trajectories and $N_{k}(t)$ the number of trajectories showing the
(photo)product $k$ at time $t$. The statistical error was derived from a
multinomial distribution as
$\epsilon_{k}=z\sqrt{\frac{1}{N}\phi_{\mathrm{k}}(t)\left(1-\phi_{\mathrm{k}}(t)\right)}$,
where $z=1.96$ is used to express a 95% confidence interval. The resulting
populations were again convoluted with the experimental pulse to account for
the temporal smearing of the pulse.

Figure 2: The different categories of (photo)products considered in our
analysis of nonadiabatic dynamics.

Finally, we note that we tested weighting the result by transition dipole
moments, but we did not observe any noticeable effects on the resulting
quantities.

###  II.5 Details of the planned experiment

We summarize here the main details of the time-resolved MeV-UED experiment
that will be conducted on cyclobutanone. A complete description of the
prediction challenge is available online.pre  The gas-phase sample of
cyclobutanone will be irradiated with 200 nm light ($\approx 80$ fs cross-
correlation), and electron-diffraction images will be obtained with an
expected 150-fs time resolution (FWHM) and 0.6 Å spatial resolution
(2$\pi$/Smax). The scattering vector S will range from 1-10 Å-1. Diffraction
images will be collected for time delays ranging from -1 ps to tens of ps in
variable step sizes. The region around time zero (-200 fs to 200 fs) will be
scanned with 30-fs step sizes. Longer positive delays will be scanned with
step sizes of up to several ps. The intensity of the 200 nm excitation light
will be kept as low as possible (5 μJ) to avoid multiphoton excitation,
exciting $\sim 10\%$ of the molecules.

###  II.6 Calculation of UED signal

The theoretical static pair-distribution function [$\text{PDF}(r)$] of all
contributing interatomic distances $r$ and time-dependent fractional change
signal ($\Delta I/I(s,t)$, with $I(s)$ the one-dimensional scattering
intensity and $s$ the momentum transfer vector) and difference pair
distribution functions [$\Delta\text{PDF}(r,t)$] signals reported in this work
were calculated within the independent atom modelMcMorrow and Als-Nielsen
(2011); Parrish and Martínez (2019) (IAM) using the nuclear geometries from
LIICs and nuclear geometries obtained along the entire swarm of FSSH
trajectories (see Ref. 76 and its SM for the detailed equations).

###  II.7 Carbon footprint of the calculations presented in this work

The simulations in this work used approximately 145’920 CPU hours (AMD EPYC
7513), drawing $\sim$2.26 MWh. Based in the Czech Republic, this consumption
leads to a carbon footprint of 1.12 T CO2e, which is equivalent to a carbon
sequestration of 101.86 tree-years or the emission of a return flight between
New York City and San Francisco (calculated using green-algorithms.org
v2.2Lannelongue, Grealey, and Inouye (2021)).

##  III Results and discussion

We start this Section by discussing a static picture of the photodynamics of
cyclobutanone (Sec. III.1), where we first confirm the nature of the excited
electronic state reached with the experimental 200 nm laser pulse. Then, we
introduce the possible photochemical pathways of the molecule in terms of
LIICs and discuss the mechanism. Finally, we present a steady-state picture of
the expected MeV-UED experimental signal produced along the two proposed
pathways. In the second part of this Section (Sec. III.2), we move to a time-
resolved picture and discuss the results of our FSSH simulations before
calculating the time-resolved signal that our nonadiabatic molecular dynamics
trajectories predict for the first 2 ps of dynamics following photoexcitation
of cyclobutanone at 200 nm.

###  III.1 Static picture of the photochemistry of cyclobutanone

####  III.1.1 Photoabsorption cross-section of cyclobutanone

The experimental 200 nm pump laser pulse targets the low-energy region of the
second absorption band in the absorption spectrum of cyclobutanone (Fig. 3a).
This absorption band is well separated from the first one, which is located at
around 280 nm, but it slightly overlaps with the third band visible in the
spectrum near 175 nm.

Let us begin by assigning the character of the electronic states responsible
for these bands and validate the selected electronic structure methods. We
calculated the first five excitation energies at the optimized geometry of the
ground-electronic state (MP2/cc-pVTZ) with various electronic-structure
methods (Table 1). The excitation to the first excited electronic state S1,
producing the low-energy absorption band, corresponds to a $n\pi^{*}$
transition located on the carbonyl moiety of cyclobutanone, with a
corresponding very low oscillator strength. The second absorption band belongs
solely to a transition towards the S2 electronic state, which exhibits a
$nRy(3s)$ character, i.e. transition from the non-bonding orbital to the
Ryberg 3s orbital. The electronic-state character will be denoted $3s$ in the
following. The electronic energy of the S2 electronic state is well separated
from that of the other states in the Franck-Condon region and appears to be
the target of our calculations – we will confirm this assumption below by
calculating the photoabsorption cross-section corresponding to this
transition. The following peak located at 175 nm is attributed to three
nearly-degenerate transitions to the excited electronic states S3-S5, all
three transitions showing a $nRy(3p)$ character. Comparing the transition
energies and oscillator strengths obtained with different electronic-structure
methods (Table 1), we observe a consistent performance between the coupled-
cluster-based techniques, XMS-CASPT2(8/8) and LR-TDDFT/TDA – XMS-CASPT2
underestimating the oscillator strength for the brighter transitions, as
observed for other molecular systems.Prlj, Hollas, and Curchod (2023) Given
our interest in further validating the nature of the second absorption band by
reproducing its vibronic structure and obtaining consistent initial conditions
with the NEA method, we chose to use LR-TDDFT/TDA/B97XD for the calculation of
photoabsorption cross-section calculations that will be reported later in this
Section. The SM contains a more extensive benchmark of the electronic-
structure methods, including also the role of triplet electronic states and a
larger basis set, and performance in the vicinity of the S${}_{2}(3s)$
minimum-energy geometry.

Transition |  Character | EOM-CC3 | EOM-CCSD |  LR-TDDFT/TDA B97XD | XMS-CASPT2(8/8)
---|---|---|---|---|---
$\Delta E$ (eV) |  $f$ |  $\Delta E$ (eV) |  $f$ |  $\Delta E$ (eV) |  $f$ |  $\Delta E$ (eV) |  $f$
S1$\leftarrow$S0 |  $n\pi^{*}$ |  4.39 |  $4\cdot 10^{-7}$ |  4.41 |  0.0000 |  4.38 |  0.0000 |  4.37 |  0.0000
S2$\leftarrow$S0 |  $nRy(3s)$ |  6.36 |  0.0380 |  6.43 |  0.0378 |  6.58 |  0.0376 |  6.25 |  0.0167
S3$\leftarrow$S0 |  $nRy(3p)$ |  6.96 |  0.0016 |  7.05 |  0.0022 |  7.15 |  0.0019 |  - |  -
S4$\leftarrow$S0 |  $nRy(3p^{\prime})$ |  7.12 |  0.0000 |  7.19 |  0.0001 |  7.33 |  0.0001 |  - |  -
S5$\leftarrow$S0 |  $nRy(3p^{\prime\prime})$ |  7.15 |  0.0021 |  7.23 |  0.0016 |  7.33 |  0.0003 |  - |  -
Table 1: Excitation energies ($\Delta E$) and oscillator strengths ($f$)
obtained at different levels of electronic-structure theory on the support of
the ground-state geometry of cyclobutanone (MP2/cc-pVTZ). All excited-state
methods employed the aug-cc-pVDZ basis set.

Before discussing the photoabsorption cross-section of cyclobutanone, it is
crucial to briefly mention its ground-state geometry. Cyclobutanone is not
planar in its ground-electronic state and, as mentioned in the Introduction,
the far-infrared and microwave spectra of cyclobutanone revealed a low
vibrational mode with a frequency of 35 cm-1 that corresponds to a ring-
puckering motion.Borgers and Strauss (1966); Scharpen and Laurie (1968);
Alonso  _et al._ (1992) A fit of the spectral lines indicates a double-well
ground-state potential for the inversion of the ring with a barrier of 4.8-7.6
cm-1. The zero-point energy of the double-well potential is approximately 9.2
cm-1 above the barrier.Scharpen and Laurie (1968) This observation is
particularly important for the sampling of the ground-state nuclear density
using a harmonic approximation. We tested different electronic-structure
methods and found a wide range of S0-minimum geometries and transition-state
energies, indicating that the ring-puckering appears to be highly sensitive to
the description of dispersion interactions offered by a given method (see SM
for a summary of these tests). Although small, these effects strongly
influence the shape of the ground-state potential energy surface and, as such,
the sampling of geometries representative of the ground-state nuclear
probability density. Our benchmark shows that MP2 was in agreement with CCSD
and XMS-CASPT2, as was DFT with the B97XD functional, which gives the best
result among the tested exchange-correlation functionals. We note that the
second lowest vibrational mode with a frequency of 395 cm-1 appears to be
harmonic.

Figure 3: Experimental and theoretical photoabsorption cross-sections of
cyclobutanone. a) Experimental photoabsorption cross-section measured by
Udvarhazi and El-SayedUdvarhazi and El-Sayed (1965) and digitized by the MPI-
Mainz UV/VIS Spectral Atlas of Gaseous Molecules of Atmospheric
Interest.Keller-Rudek  _et al._ (2013) The energy distribution of the
experimental pump laser pulse (200 nm, $\approx 80$ fs cross-correlation) is
superimposed to highlight the region of photoexcitation (see Sec. II). b)
Calculated S2$\leftarrow$S0 photoabsorption cross-section for cyclobutanone
using the FCHT and the NEA strategies. In both cases, excitation energies and
all required electronic quantities were calculated with LR-TDDFT/TDA/B97XD.
For the NEA, 5000 geometries were sampled from QT-BOMD. The highlighted area
(ICs) represents the energy window employed to select the initial conditions
for the subsequent nonadiabatic molecular dynamics, mimicking the role of the
experimental laser pulse. No shifts were applied to the calculated
photoabsorption cross-sections.

We now move to the simulation of the photoabsorption cross-section for the
second band of cyclobutanone. To reproduce the experimentally-observed
vibronic structure of this band (Fig. 3a), we first used the FCHT method based
on DFT/B97XD (ground-state optimized geometry) and LR-TDDFT/TDA/B97XD (the
minimum of the S2 electronic state, transition energy, transition dipole
moment and its nuclear derivative). While the S0 ring-puckering mode is not
adequately described within the harmonic approximation, the target S2
electronic state is harmonic. Therefore, the vibrational splitting in the
photoabsorption cross-section should not be affected as it depends only on the
vibrational energy levels of the S2 electronic state. The ground-state
anharmonicity of the ring puckering mode may, however, impact the absolute
shift of the photoabsorption cross-section and its intensity. The FCHT
photoabsorption cross-section is presented in Fig. 3b and, despite the
potential anharmonicity issue, shows an excellent agreement with the
experimental cross-section, capturing the first few peaks of the vibronic
progression with accuracy. This result further confirms the electronic
character of this band and validates the electronic-structure methods employed
in this work. We stress that no shifts were applied to the calculated
photoabsorption cross-sections presented here. The modes responsible for the
peak in the photoabsorption cross-section at around 200 nm – excited by the
pump pulse – correspond to a ring puckering and C=O stretching.

We also calculated the photoabsorption cross-section using the NEA, which
proposes to use the reflection principle based on a set of sampled geometries
representative of the ground-state nuclear probability density. The NEA is
faster to deploy than the FCHT but does not give access to the vibronic
structure, as the nuclear wavefunctions for the excited electronic state of
interest are not calculated. The sampling of nuclear coordinates required by
the NEA means that this technique also offers a simple strategy to select
initial conditions for any subsequent nonadiabatic molecular dynamics (nuclear
momenta are then also required). The ground-state probability density was
sampled using QT-BOMD (see Sec. II), which is able to describe anharmonic
potential energy surfaces such as the ring puckering mode of cyclobutanone.
QT-BOMD appears to outperform the other common sampling strategies like the
harmonic Wigner distribution for flexible molecules.Prlj, Hollas, and Curchod
(2023) The QT-BOMD simulation used MP2/cc-pVDZ for the ground-state
electronic-structure quantities, while the excitation energy and transition
dipole moment required for the NEA was obtained with LR-TDDFT/TDA/B97XD for
each sampled geometry. The resulting NEA photoabsorption cross-section was
constructed from 5000 geometries and recovers the low-energy spectral envelope
of the FCHT photoabsorption cross-section (Fig. 3b). To select the initial
conditions (nuclear positions and momenta) for the nonadiabatic molecular
dynamics simulations, we used an energy window (box in Fig. 3b) mimicking the
excitation laser pulse to filter the geometries in the photoexcitation region
of interest.

####  III.1.2 Possible relaxation pathways of cyclobutanone upon
photoexcitation at 200 nm

Upon photoexcitation with a 200 nm laser pulse, our earlier results indicate
that cyclobutanone will be promoted into its S${}_{2}(3s)$ electronic state.
After an initial relaxation away from the Franck-Condon region and towards the
planar S2 minimum, two different relaxation pathways appear to be possible to
reach S1 and S0 (these pathways will be further confirmed by our nonadiabatic
molecular dynamics simulations presented below). To map the energy profile
along these pathways, we optimized the critical geometries characterizing each
pathway – energy minima (min) and conical intersections (CIs) – at the XMS-
CASPT2(8/8)/aug-cc-pVDZ and connected these geometries by interpolated
structures using LIICs. We then used XMS-CASPT2(8/8) to obtain the electronic
energies for S0, S1, and S2 on the support of the LIICs geometries – see Fig.
4 for the final result (the electronic-structure benchmark validating the use
of XMS-CASPT2(8/8) can be found in the SM). The first pathway – coined the ’CI
pathway’ (left part of Fig. 4) – shows how cyclobutanone can funnel from
S2(min) to S1 preserving its closed ring via the CI-closed(S2/S1) conical
intersection (0.58 eV above S2(min)), after which the molecule can relax
towards the S1 minimum. From the S1 minimum, cyclobutanone can open the ring,
overcoming a barrier of approximately 0.18 eV in S1, and hit the CI(S1/S0)
conical intersection, which takes the molecule back to the ground electronic
state, leading to the formation of photoproducts or possibly the reformation
of cyclobutanone via ring closure.

Figure 4: Two possible relaxation pathways upon photoexcitation of
cyclobutanone in S2 (vertical arrow) from its S0 minimum. Electronic energies
(XMS-CASPT2(8/8)/aug-cc-pVDZ) are represented along the LIICs characterizing
two relaxation pathways – the CI pathway [left of S0(min)] and the adiabatic
pathway [right of S0(min)]. As LIICs do not accurately describe barrier
energies, we denote the energies of the true transition state with a star. A
shaded area highlights the electronic state likely to drive the dynamics of
the molecule from S2 to S0. The molecular structure representative of each
located critical point is superimposed on the pathways (the spatial
orientation of each structure was chosen to highlight the specific distortions
of the molecular framework).

The second nonradiative pathway involves a change of electronic character
along the S2 electronic state, leading to the direct ring opening of
cyclobutanone after overcoming a barrier of about 0.15 eV – this pathway is
coined the ’adiabatic pathway’ (right part of Fig. 4). We note that the S3
state stays energetically separated from S2 in the proximity of the S2
transition state, limiting the possibilities for a nonadiabatic transition
towards S3 (see SM). Since the activation barrier for this process is smaller
than the activation barrier for the CI pathway, one may expect that this
pathway dominates the nonradiative decay of cyclobutanone. Once the ring-
opened molecule crosses the S2 transition state, it decays toward the CI-
open(S2/S1) and can then funnel through CI(S1/S0) (the same conical
intersection as for the CI pathway). Interestingly, the ring-opening process
in S2 means that the three electronic states considered become close in energy
in the direct vicinity of CI-open(S2/S1) (see electronic energies near CI-
open(S2/S1) in Fig. 4). This proximity between S2, S1, and S0 may lead to
efficient transfers through the three states in this region of configuration
space. We will see that this prediction is corroborated by the nonadiabatic
molecular dynamics simulations described below.

At this stage, we should comment on the possibility of ISC processes in the
photodynamics of cyclobutanone following photoexcitation at 200 nm. We
calculated the SOC magnitude between the driving singlet excited state and
nearby triplet electronic states with the SA3-CASSCF(8/8)/aug-cc-pVDZ method
along our LIICs (see SM) and observed that the SOC magnitude remains smaller
than 20 cm-1 in the vicinity of the S2 minimum. Such a value of SOC would
require the photoexcited molecule to exhibit an extended lifetime in the S2
electronic state, i.e., picoseconds,Favero, Granucci, and Persico (2013) to
undergo an efficient deactivation via ISC. By analogy with the photodynamics
of cyclobutanone excited in the S1 electronic state, a cyclobutanone with a
long-lived S2 state can be prepared by photoexcitation below the dissociation
barrier. In the S1 state, excitation above 315 nm (red edge of the first band
in the photoabsorption cross-section) creates a photoexcited molecule with
internal energy below the dissociation barrier, opening the ISC channel for
deactivation. We note that the 200 nm laser pump that will be used to
photoexcite cyclobutanone will also target the red edge of the photoabsorption
cross-section (see Fig. 3a). Therefore, a similar scenario could be plausible
for the S2 ring-opening barrier, as the S2 barrier height is 0.145 eV (both
with and without zero-point energy, ZPE) and comparable to that of the S1
state, 0.176 eV without ZPE and 0.146 with ZPE. Thus, an ISC channel might be
operational on a long timescale during the planned MeV-UED experiment,
hampering the IC process. However, earlier experimental studiesKuhlman,
Sølling, and Møller (2012) at 200 nm (and our simulations presented below)
reveal an ultrafast internal conversion of cyclobutanone from the S2
electronic state and cast doubt on the ISC hypothesis. We finally note that
for the ring-opened S1 structure, the SOC magnitude between the S1 state and
neighboring triplet states rises up to 40 cm-1 yet we will show below that the
lifetime of this species in S1 is again too short for ISC to occur.

####  III.1.3 Steady-state UED signals of cyclobutanone

In the following, we offer a first prediction of a quantity that can be
retrieved from the MeV-UED experiment, the pair-distribution function (PDF),
on the support of the two reaction pathways we discussed earlier.

Let us begin by inspecting the PDF for the ground-state cyclobutanone (Fig. 5)
resulting from the ground-state optimized geometry and that obtained as an
average over the PDFs for each initial condition (used to initiate our
nonadiabatic molecular dynamics). We note that the PDF acquired for the 119
initial conditions is in close agreement with the PDF obtained from the 5000
geometries sampled from the QT-BOMD trajectory and, therefore, is believed to
offer a reliable depiction of a PDF for the S0 nuclear probability density of
cyclobutanone, see SM. First, we decompose the average PDF obtained for the
initial conditions into contributions of different atom pairs observed in the
molecule (see Fig. 5 and the molecular structure for the color code). The
first peak of the overall PDF for the initial conditions (centered around 1.5
Å) corresponds to bonded pairs of atoms – the shortest distance comes from the
carbonyl bond (red curve in Fig. 5), while the main contribution (lime curve)
originates from all the C-C bonds of the molecular framework. The second peak
of the main PDF for the initial conditions, centered at around 2.5 Å, is a
combination of non-bonded atom pairs, mainly composed of distances between the
carbon or oxygen of the carbonyl with the three opposite carbons in the
cyclobutanone ring. Interestingly, while the second peak of the PDF for the
ground-state optimized geometry exhibits a lower intensity than the first
peak, this order is reversed in the PDF obtained from the whole set of initial
conditions. As stressed before, the ground-state minimum is not fully
representative of the ground-state nuclear probability density, and the non-
bonded distances appear to play a more significant role when the nuclear
density is considered. Thus, care needs to be taken when using optimized
geometries to analyze an experimental UED signal, and adequate nuclear
probability density should be considered instead (reinforcing the message of
Ref. 76).

Figure 5: Steady-state PDF of cyclobutanone. The PDF was calculated for the
ground-state minimum (XMS-CASPT2(8/8)/aug-cc-pVDZ) and as an average of the
PDF calculated for each initial condition (119) used for the nonadiabatic
dynamics (sampled from the QT-BOMD trajectory, MP2/cc-pVDZ). The PDF signal
obtained from the initial conditions is decomposed into different atom-pair
distances relevant to cyclobutanone, highlighted in the molecular structure
given as inset. Note that the PDF signal depicted by the lime curve
corresponds to all the C-C bonds of cyclobutanone.

Let us now move to a more conceptual analysis of our data by predicting the
PDF signal along the two LIIC pathways discussed in Fig. 4. Needless to stress
that such an analysis needs to be taken with a grain of salt based on the
previous paragraph as it solely relies on critical geometries and interpolated
geometries between them – most likely not representative of a time-resolved
signal that we will discuss later based on nonadiabatic molecular dynamics
simulations. The difference PDFs ($\Delta$PDFs, obtained as the PDF for a
given geometry on the LIIC minus the PDF of the reference S0 minimum geometry)
for both LIICs introduced earlier are presented in Fig. 6. The initial step of
both pathways is the relaxation of cyclobutanone from its S0 minimum geometry
to the S2 minimum geometry, which manifests itself as a fringe of positive and
negative $\Delta$PDF signals caused by the elongation of all distances between
the carbonyl group and the other carbon atoms due to planarization. This
elongation leads to negative $\Delta$PDF signals for each of the original S0
minimum atom-pair distances (shown in Fig. 5) and positive $\Delta$PDF signals
for the now increased atom-pair distances, hence the appearance of fringe on
the overall $\Delta$PDF plot for this part of both pathways.

Figure 6: Simulated $\Delta$PDF signal (PDF for a given geometry along the
LIIC minus the PDF of the reference S0 minimum geometry) along the two
deactivation pathways – a) the CI pathway and b) the adiabatic pathway – of
cyclobutanone excited at 200 nm. The geometries used to predict $\Delta$PDF
signals are the same as those used for the electronic energies presented in
Fig. 4, with the only exception of the transition-state geometries (TS) along
the CI pathway that was added for a more detailed analysis of the ring-opening
process.

We shall now focus on the CI pathway. The fringe pattern persists during the
passage through CI-closed(S2/S1) and until cyclobutanone reaches the S1
minimum. The only difference in this part of the path appears around 1 Å. This
change is caused by the C=O bond that first shortens from 1.17 Å to 1.12 Å
when going from the S2 minimum to the CI-closed(S2/S1) – explaining the
positive $\Delta$PDF signal observed near 1 Å around CI-closed(S2/S1). The C=O
bond then stretches to 1.30 Å when reaching the S1 minimum, consistent with
the $n\pi^{*}$ character of this state and leading to a depletion of the
$\Delta$PDF signal at the original S0 C=O bond distance and the appearance of
a negative $\Delta$PDF feature at around 1.20 Å. The following ring-opening
process over the transition state leads to a substantial depletion of the
$\Delta$PDF signal around 2 Å and a connected increase around 3 Å when the
molecule reaches the CI(S1/S0), clearly attesting from the loss of a bonded
interaction and shorter non-bonded atom-pair distances, and the appearance of
further (non-bonded) distances from the ring-opened structure. Focusing on the
adiabatic pathway, we observe a very similar pattern to that of the CI pathway
when the cyclobutanone ring opens, simply arriving earlier in the deactivation
pathway as the S2 TS already results in a ring opening. The overall
$\Delta$PDF pattern barely changes between CI-open(S2/S1) and CI(S1/S0) as
both geometries are open and rather similar geometrically. We note that the
evolution of different bond distances depicted in the inset of Fig. 5 is
presented in the SM for more insight on the $\Delta$PDFs along the LIICs.

The differences between the pathways do not appear to be very pronounced, and
distinguishing between them in the experiment might be difficult. The
evolution of the photoexcited cyclobutanone on S2 is likely to be visible, and
an interesting question will be the time the molecule spends near the S2
minimum before finding its way to the CI-closed(S2/S1) or the S2 TS. The
vibronic structure observed in the photoabsorption cross-section (Fig. 3) and
the broadening of its peaks seem to indicate that cyclobutanone should at
least remain in this region of configuration space for a few hundred
femtoseconds, which is corroborated by earlier time-resolved photoelectron
spectroscopy and time-resolved mass spectroscopy experiments.Kuhlman, Sølling,
and Møller (2012) While the ring opening of cyclobutanone will be clearly
visible when it happens, both pathways – ring opening in S1 or S2 – lead to
very similar features in the $\Delta$PDF. Probably the most significant
difference between the pathways lies in the variation of the $\Delta$PDF
signal for molecular structures near the CI-closed(S2/S1) and S1 minimum.
While the CI-closed(S2/S1) will only act as a gateway (the nuclear wavepacket
will only visit this part of configuration space for a short amount of time),
the S1 minimum could be sufficiently long-lived to be observed. The S1 minimum
is characterized by an extended C=O bond and a shorter distance between the
carbonyl group and the opposite carbon. As a result, the fringe is shifted to
shorter distances compared to the S2 minimum (see Fig. 6a).

###  III.2 Time-resolved photochemistry of cyclobutanone

####  III.2.1 Analysis of the nonadiabatic molecular dynamics simulations

Figure 7: Relative populations of the photoproducts (left) and the electronic
states (right) obtained from the (NA+BO)MD for a swarm of 119 trajectories.
The population traces were convoluted with the envelope of the laser pulse.
The values report the fitted lifetimes for the population of cyclobutanone in
an excited electronic state ($\tau_{\mathrm{CB}}$) and the open structure
($\tau_{\mathrm{open}}$), as well as the electronic population of the S2 state
($\tau_{\mathrm{S_{2}}}$) and the S1 state ($\tau_{\mathrm{S_{1}}}$). The
details on the fitting procedure and resulting curves are described in the SM.

We begin this discussion of our (NA+BO)MD simulations of cyclobutanone,
photoexcited in its S2 electronic state, by inspecting the resulting
populations predicted by the swarm of trajectories. Let us first focus on the
time-dependent relative populations of the cyclobutanone and photoproducts
(left panel of Fig. 7). The population of the parent cyclobutanone molecule
rapidly decays (lifetime of $\sim$359 fs) following photoexcitation due to a
ring-opening induced by the cleavage of one of the C-C bonds adjacent to the
carbonyl group. This open structure acts as an intermediate to the formation
of either CO + C3 (cyclopropane or propene) or ethene + ketene (C2), or the
reformation of a hot ground-state cyclobutanone with a ratio of 71:42:4. The
C3 products observed are dominated by cyclopropane (83%), followed by propene
(17%). The formation of cyclopropane or propene is mediated by the appearance
of an intermediate ^.CH2-CH2-CH2^. biradical structure. We did not observe the
conversion of cyclopropane to propene (or vice versa) in our simulation, i.e.,
the propene was always formed directly from the photochemical reaction and not
later in the ground state via the isomerization from cyclopropane, as
suggested in the literature. For the C2 channel, we observed the decomposition
of ketene into CH2 and CO in the ground electronic state as a minor channel
(only 2 trajectories). A small population of hot ground-state cyclobutanone
(S${}_{0}^{\\#}$) is reformed during the photochemical process, yet it was
observed to be unstable and decomposed into one of the other photoproducts
after several hundreds of femtoseconds. The relative populations of the
photoproducts after 2 ps are summarized in Table 2.

photoproducts | relative population in %
---|---
CO + cyclopropane | $49.6\pm 9.0$
CO + propene | $10.1\pm 5.4$
CO + biradical | $0.0\pm 0.0$
ethene + ketene | $33.6\pm 8.5$
ethene + CO + CH2 | $1.7\pm 2.3$
cyclobutanone | $1.7\pm 2.3$
cyclobutanone(S${}_{0}^{\\#}$) | $3.4\pm 3.2$
Table 2: Relative populations of the photoproducts of cyclobutanone following
irradiation at 200 nm as predicted after 2 ps of (NA+BO)MD. We note that these
populations vary over time and have to be considered as a snapshot of the
photoproduct populations in the earlier times following the photoexcitation.

The decay of the electronic-state populations exhibits similar time scales,
with the lifetime of the S2 population predicted to be $\sim$334 fs and that
of the S1 state $\sim$66 fs. While the majority of the trajectories proceeded
via the expected cascade
$\mathrm{S}_{2}\rightarrow\mathrm{S}_{1}\rightarrow\mathrm{S}_{0}$, 12 out of
the 119 trajectories hopped directly from S2. to S0 These hops always happened
during the adiabatic pathway close to the CI-open geometry (see Fig. 4), where
all three electronic states are coupled, and their respective FSSH electronic
coefficients are non-zero.

From a mechanistic perspective, our (NA+BO)MD trajectories indicate that both
the adiabatic and CI pathways are operative, with the adiabatic pathway being
the dominant channel. A detailed analysis of the hopping geometries between
the S2 and S1 states using a multidimensional scaling technique indeed reveals
that the hops through the CI-open (adiabatic pathway) versus those through the
CI-closed (CI pathway) follow a 3:1 ratio (see SM for additional details).
Apart from these two pathways, we observed a few additional minor channels:
the molecule can dissociate into ethene and ketene in the S1 electronic state
and then transfer to the ground state, or the molecule can remain open in the
ground state for a few hundreds of femtoseconds before dissociating and
forming one of the products.

Finally, we note that we also ran 45 trajectories using an LZSH procedure
(using the same level of electronic-structure theory) to probe the sensitivity
of the dynamics to the nonadiabatic dynamics algorithm. We did not observe any
qualitative discrepancy for the 45 trajectories within 700 fs of dynamics,
reinforcing our confidence in the FSSH trajectories. We also ran a set of 10
LZSH/SA(3S,3T)-CASSCF(8/8) trajectories, including the triplet states, to
confirm that the ISC pathway should be minor for photoexcitation of
cyclobutanone to its S2 electronic state, see SM for more details.

####  III.2.2 Time-resolved MeV-UED signal

We present here the prediction of time-resolved and time-averaged MeV-UED
signals based on the swarm of (NA+BO)MD trajectories. The $\Delta I/I(s,t)$
signal (Fig. 8) was evaluated from the molecular structures along the
(NA+BO)MD trajectories using the IAM and convoluted with the instrument
response function (150 fs FWHMShen  _et al._ (2019)). We note that this
convolution smears out some fine details of the signal, like the early-time
oscillations produced by the photoexcited cyclobutanone in the S2 electronic
state. The total $\Delta I/I(s,t)$ signal (left panel of Fig. 8) only exhibits
weak variations over the course of time, in particular after 500 fs, and
displays two strong positive ($\sim$ 2.5 and 7.5 $\mathrm{\AA}^{-1}$) and
negative (between 0 and 1.5 $\mathrm{\AA}^{-1}$) features. The rest of the $s$
range considered is mostly composed of weak negative contributions. As
introduced in Ref. 76, we exploited one advantage of theory over experiment:
the total $\Delta I/I(s,t)$ can be decomposed into the contributions of each
different species (right panels of Fig. 8), providing an alternative
interpretation tool for the features observed in the total signal. From this
decomposition, we can directly observe that the contribution of the
photoproducts (ethene + ketene and C3 \+ CO) to the $\Delta I/I(s,t)$ signal
is about 5 times stronger than that of cyclobutanone or its open structure and
explain the strong features highlighted before.

Figure 8: Time-resolved $\Delta I/I(s,t)$ MeV-UED signal obtained from the
swarm of (NA+BO)MD trajectories and convoluted with the instrument response
function. The total $\Delta I/I(s,t)$ signal is presented on the left panel
and then decomposed into the contribution of cyclobutanone and each category
of photoproducts.

Given the rather time-independent nature of the decomposed $\Delta I/I(s,t)$
UED signals in the region between 1 and 2 ps, we average them to produce
representative $\Delta I/I(s)$ signals for each category of photoproducts and
cyclobutanone (left panel of Fig. 9). These signals can be used as time-
independent basis functions to fit the corresponding experimental signal and
decompose it into the contributions of the main different molecular structures
(as performed in Ref. 76). The features (amplitudes and locations of the
peaks) of each time-independent $\Delta I/I(s)$ basis function produced are
different enough to be used in a fitting procedure. We note that
distinguishing between the different C3 products (cyclopropane and propene)
would most likely not be possible experimentally due to experimental
resolution – this observation motivated the categories proposed for our
general analysis of the (NA+BO)MD results. We also deployed the same procedure
to produce averaged $\Delta$PDF$(r)$ for the photoproducts and cyclobutanone –
these signals being easier to interpret in terms of molecular parameters as
discussed later.

Figure 9: Representative $\Delta I/I(s)$ (left) and $\Delta$PDF$(r)$ (right)
signals for cyclobutanone and the different families of photoproducts,
obtained from the averaging of the theoretical $\Delta I/I(s,t)$ and
$\Delta$PDF$(r,t)$ signals for each category of molecules in the time interval
1 - 2 ps.

Let us finally move to the interpretation of the predicted time-resolved
$\Delta$PDF$(r,t)$ signal calculated from our (NA+BO)MD trajectories (Fig. 10)
by using the basis functions introduced above (right panel of Fig. 9).
Focusing on the short distances, we observe a constant weak positive signal at
around 1 $\mathrm{\AA}$ from the early times of the dynamics, corresponding to
a loss of single carbon bonds in cyclobutanone and formation of double bonds
in ethene and ketene (see the basis function for ethene + ketene in Fig. 9).
Centered at around 2 $\mathrm{\AA}$, we observe two strong negative features
caused by a loss of chemical bonds during dissociation. This loss can be
equally attributed to the dissociation into C2 and C3 photoproducts or the
open structure, see Fig. 9. While the described signal below 2.5
$\mathrm{\AA}$ is constant throughout the entire dynamics, the region above
2.5 $\mathrm{\AA}$ shows time-dependent features. Within the first 250 fs of
dynamics, two positive signals between 3 and 4 $\mathrm{\AA}$ appear due to
the formation of the open structure, which is consistent with its relative
population calculated in Fig. 7. This positive signal appears then to spread
to larger distances $r$ as the fragments dissociate, and the signal in the
region between 3 and 4 $\mathrm{\AA}$ becomes negative due to the formation of
the photoproducts. Since the positive signal characteristic of the
dissociation is low, we expect that a direct observation of the dissociation
features on a long timescale might be challenging. With that said, Figure 9
shows that dissociation produces a very strong feature at low $s$, which can
be clearly identified in the UED signal, even at later times. To distinguish
between dissociation and ring-opening, it would be important that the low $s$
signal ($0.5>s>1$) is recorded.

We note that the fringe pattern representative of the S2 minimum predicted by
our static calculations along the LIICs can be visible only in the
unconvoluted signal during the first 40 fs. Nevertheless, it’s quickly
overlayed by the ring opening and the convolution with the instrument response
function smears it out. Still, the fringe pattern might be visible in the
experiment considering the potential limitations of our simulations (see Sec.
IV below), as the ring-opening of cyclobutanone in the excited state may be
too rapid in our simulations (due to the accuracy of the electronic-structure
method and the internal energy of the photoexcited cyclobutanone).

Figure 10: Time-resolved $\Delta$PDF$(r,t)$ signal obtained from the swarm of
(NA+BO)MD trajectories and convoluted with the instrument response function.
Contour lines around the signal map were added to improve the visibility of
weak features.

##  IV Limitations of the current model and their impact on the predicted
experimental signals

In this Section, we summarize the limitations of the model used for our
predictions, arising either from the approximations underlying our methodology
or the time constraint imposed by this prediction challenge. We grouped the
limitations into four categories corresponding to the various components of
our modeling: level of electronic-structure theory, the influence of initial
conditions, approximations of the nonadiabatic molecular dynamics, and
generation of the UED signal. We also discuss the implications of these
limitations on the predictions of the experimental signals and the
corresponding “confidence intervals”.

###  IV.1 Electronic-structure theory

The electronic structure has been recognized as one of the most quality-
determining components of nonadiabatic molecular dynamics simulations. We
selected the XMS-CASPT2 method for our simulations due to i) the homolytic
bond breaking that requires a multiconfigurational method and ii) the
sensitivity of the activation barriers for dissociation to dynamic correlation
observed in our LIIC. All the other methods at hand (LR-TDDFT, ADC(2), SI-SA-
REKS, SA-CASSCF, or FOMO-CASCI) were either not compatible with a homolytic
dissociation or would not provide an accurate depiction of the energy barriers
along the dissociation profile. Nevertheless, the performance of XMS-CASPT2
strongly depends on the choice of its parameters, such as the active space. We
offered a detailed benchmark to justify our methodology, but with more time at
hand, we would have investigated further the following potential improvements
of our XMS-CASPT2 calculations.

(a) Inclusion of the two $\sigma_{\mathrm{CC}}$ and corresponding two
$\sigma^{*}_{\mathrm{CC}}$ orbitals opposite to the carbonyl group, which
facilitate the formation of the ketene and ethene photoproducts. Our active
space does not include these orbitals from the beginning, and our simulations
rely on (slowly) rotating these orbitals in the active space once the adjacent
C-C bond is open and its corresponding $\sigma$ orbitals become irrelevant.
Therefore, the formation of ketene and ethene depends on uncontrolled orbital
rotations and can be hampered compared to the release of CO, for which all the
relevant orbitals are present in the active space from the beginning of the
dynamics. This extended (12/12) active space would require more computational
time than this challenge allowed and was, as such, deemed out of reach.

(b) Inclusion of higher excited electronic states in the state averaging to
stabilize the active space and prevent jumps in total energy. The issues with
energy discontinuities in our simulations partially stemmed from the S2
electronic state increasing in energy upon dissociation (due to the character
of this electronic state), leading eventually to a root flipping. However, as
for point (a) above, this increased state-averaging protocol would require the
inclusion of $Ry(3p)$ orbitals in the active space, meaning an even larger
active space and prohibitive computational costs for this challenge.

(c) Inclusion of the triplet electronic states in the state averaging. This
inclusion would allow us to monitor the influence of SOC along the
trajectories (more rigorously than our tests with LZSH). However, the (8/8)
active space used in this work is unstable with respect to the T2
($\pi\pi^{*}$) electronic state, which quickly rises in energy and swaps with
other triplet electronic states, causing a jump in total energy.

###  IV.2 Influence of initial conditions

The selection of initial conditions in photodynamics simulations is a topic of
intense research in the community and the modeled experiment represents a
challenging (and exciting) case in this context. Our modeling relies on a
sudden vertical excitation of the molecule, projecting the ground state on the
S2 excited electronic state, combined with an energy-windowing scheme, that
accounts for the energy spread of the laser pulse. The main limitations of
this approach are:

(a) The vertical excitation approximation is unable to correctly describe
excitation to well-defined vibrational states, which is expected with the
experimental pump laser pulse. Although we used the windowing scheme, our
simulations still remain within the frame of a sudden excitation of the
molecule, which is justified only for a $\delta$ laser pulse. As a result, the
spread of total energy for our simulated nuclear wave packet at the beginning
of the nonadiabatic dynamics is about 0.84 eV, while it should be around 0.05
eV for the experimental pulse. This observation is crucial if one considers
that the photoexcited cyclobutanone has to cross a barrier following
photoexcitation to the S2 electronic state. The excess of internal energy in
our initial nuclear wave packet may result in a shorter lifetime of the
excited molecule in S2 and hamper the ISC channel.

(b) The sudden approximation used in this work leads to the formation of a
perfect nuclear wavepacket in S2, meaning that the experimental signal
simulated is most likely too coherent. A more suitable generation of initial
conditions can be achieved by a technique developed by Saalfrank et
al,Martínez-Mesa and Saalfrank (2015) which includes the laser pulse
explicitly to provide time-dependent initial conditions. However, an
implementation of this strategy for molecular systems is currently not
available.

(c) The ground-state nuclear density was sampled at 298 K, as the true
temperature after a nozzle expansion was not clearly defined.

###  IV.3 Approximations of the nonadiabatic molecular dynamics strategy

Our results are based on FSSH trajectories that were benchmarked against the
simpler LZSH technique. Although both methods yielded similar population
traces, several limiting factors of this trajectory-based approach need to be
highlighted.

(a) FSSH may suffer from its independent trajectory approximation. We observed
for some of our trajectories that the FSSH electronic coefficients were non-
zero for tens of femtoseconds following nonadiabatic transitions. The
decoherence is accounted for via the ad hoc energy-based decoherence
correction in the present study, yet this correction appears to compete with
the transfer of electronic amplitude during the dynamics through the
nonadiabatic regions. This observation calls for additional tests on the
reliability of decoherence-corrected FSSH, for example by comparing its result
with that obtained with a multiple-spawning based approach like AIMSWISS.

(b) The three-state hops observed in our dynamics (12 out of 119 trajectories)
might also pose a challenge for the FSSH algorithm, and further benchmarking
with a multiple-spawning-based approach like AIMSWISS would be desirable.

(c) The simulations do not incorporate ISC processes, rationalized by our
static calculations and tests using LZSH/CASSCF. It would be safer to monitor
the magnitude of SOC along the trajectories and preserve the possibility for
ISC channels. As stated above, this approximation is mostly motivated by the
availability and cost of the electronic structure and the potential need for a
large number of trajectories to describe ISC adequately.Persico and Granucci
(2014)

(d) Zero-point energy leakage can affect our semiclassical dynamics on the 2
ps timescale.Pieroni  _et al._ (2023) As a result, our lifetimes are likely to
be shorter than predicted by quantum dynamics.

###  IV.4 Generation of the UED signal

Our prediction can also be affected by the IAM applied to approximate the UED
signal from our FSSH trajectories. Although the effect of the electronic
charge density is often negligible, it was shown that it could play a
significant role in the evaluation of signal.Champenois  _et al._ (2023)

###  IV.5 Direct implications of these limitations on the predicted MeV-UED
signals

Based on the limitations highlighted below, we would expect to see a longer
lifetime of cyclobutanone in its photoexcited S2 electronic state. The time-
dependent relative populations of the photoproducts obtained should be
considered with larger error bars than the statistical ones provided. We
expect that the active space used favours the CO release over the formation of
ethene and ketene, although other effects also might play a role. The ratio of
photoproducts and lifetimes predicted might also differ because the pump laser
pulse will most likely form the excited-state wave packet with lower internal
energy than our sudden excitation estimate, which may appear at the edge of
the dissociation barrier and CI-closed. If this happens, we would then expect
the MeV-UED signal representative of the S2(min) to last for longer and,
possibly, to observe an alternative ISC channel on a scale of a few to tens of
picoseconds. The ISC channel would dominantly form CO and
cyclopropanone/propene based on the experiments with 340 nm pump pulseTang and
Lee (1976). We note, however, that a previous experiment using the same
excitation wavelength observed subpicosecond deactivation processes.Kuhlman,
Sølling, and Møller (2012)

##  V Conclusion

In this work, we investigated the first two picoseconds of the cyclobutanone
photodynamics following an excitation with a laser pulse centered at 200 nm
and predicted the expected time-resolved MeV-UED signal for the processes
observed. Our (NA+BO)MD trajectories show an almost complete return of the
nuclear wavepacket to the ground state within the timescale of our
simulations. The main photoproduct observed from this photochemical reaction
after 2 ps is CO (60%) accompanied by either cyclopropane (50%) or propene
(10%), while the second photoproduct was ethene and ketene (34%). An
additional minor product was detected: ethene with CO and CH2, which is formed
by a dissociation of the ground-state ketene and a hot ground-state
cyclobutanone. The formation of photoproducts follows two main pathways. In
the adiabatic pathway, the photoexcited cyclobutanone ring opens next to the
carbonyl group in the S2 electronic state to form an open biradical structure,
which then decays rapidly via internal conversion to the ground electronic
state S0, where the photoproducts are formed. In the CI pathway, the
photoexcited cyclobutanone undergoes an internal conversion to S1 through a
conical intersection characterized by a closed cyclobutanone structure. The
C-C bond next to the carbonyl group is then cleaved and photoproducts are
formed upon internal conversion to the S0 state. Our (NA+BO)MD trajectories
predict that the ratio of the adiabatic to CI pathways is approximately 3:1.
The decomposition of cyclobutanone is characterized by a lifetime estimated to
be $359\pm 13$ fs, which correlates with the S2-state population lifetime of
$334\pm 12$ fs. Similarly, the lifetime of the ring-open structure ($114\pm
11$ fs) correlates with the lifetime of the population in the S1 state ($66\pm
10$ fs).

The main outcome of this work is the prediction of a series of MeV-UED signals
for the photochemistry of cyclobutanone excited at 200 nm. A stationary MeV-
UED signal ($\Delta$PDF) was predicted along interpolated paths in internal
coordinates, characterizing the critical structures relevant for the adiabatic
and CI pathways. Our swarm of (NA+BO)MD trajectories was also used to predict
time-resolved MEV-UED signals ($\Delta$PDF and $\Delta I/I$) that would be
observed for the 2 ps nuclear wavepacket dynamics simulated in our work. We
decomposed the $\Delta I/I$ signal into its molecular components
(cyclobutanone and the different photoproducts observed) and suggested time-
independent basis functions – based on the average, athermal structures
observed in the dynamics – that could be used to distinguish between the
following molecules when they are evolving in the hot ground state: hot
cyclobutanone, C3 products + CO, ethene + ketene, and the ring-open structure.

We further stress the challenges one faces when predicting such experiments
and the expected limitations of our approaches, which are difficult to
overcome either due the current limitations of the methodologies or the time
constraints imposed by this prediction challenge. As emphasized in recent
work,Janoš and Slavíček (2023) the level of electronic-structure theory
appears to play a decisive role in the quality of the predictions, both
stationary and time-resolved ones, given in particular the sensitivity of the
excited-state dissociation barrier in the S2 electronic state to the quantum-
chemical method used and its importance in controlling the efficiency of the
ring-opening process. This dissociation process appears to require both a
multiconfigurational treatment and a proper inclusion of dynamic correlation,
leading us to employ the XMS-CASPT2 method, which, however, is highly
sensitive to the precise selection of its underlying active space and the
number of electronic states considered, leading to very computationally-
demanding simulations. The price to be paid for using XMS-CASPT2 forced us to
compromise and combine it with the mixed quantum/classical method FSSH, which
may not offer the most accurate depiction of nonadiabatic transitions.

Another challenge we experienced was the adequate description of the initial
conditions for the (NA+BO)MD simulations. Since the ring-puckering mode
exhibits a double-well potential and is highly anharmonic, the harmonic
approximation underlying the usual (harmonic) Wigner sampling may lead to
artifacts (as observed recently for another anharmonic moleculePrlj, Hollas,
and Curchod (2023)). More involved approaches to sampling the ground-state
density are necessary, like ab initio Born–Oppenheimer molecular dynamics
coupled to a quantum thermostat used in the present work. The effects of the
laser pulse are also nontrivial to include when dealing with a narrow spectral
bandwidth of about 0.05 eV. The usual vertical approximation, even with an
energy windowing, remains a crude approximation to the real excitation process
and may lead to a photoexcited molecule with higher internal energy –
problematic when the photochemical processes involve a barrier in the excited
electronic state. We believe that further methodological work on sampling the
ground-state density, including the effect of a laser pulse explicitly, and
electronic-structure theory are necessary to address these challenges.

While uncomfortable, the time constraint imposed by this prediction challenge
to perform our simulations is actually a common scenario when collaborating
with experimentalists, who may require rapid theoretical support to analyze
their data. Thus, this situation overall calls for the development of
diagnostic tools that would at least allow us to test the sensitivity of fast
predictions for more reliable interpretations of experimental results.

## Supplementary Material

The Supplementary Material contains a depiction of the active-space orbitals,
a detailed electronic-structure benchmark for the ground- and excited-
electronic states, an analysis of the higher excited electronic states of
cyclobutanone and their potential influence on the nonadiabatic dynamics, the
detailed classification of photoproducts for the analysis of the nonadiabatic
molecular dynamics, an analysis of the influence of spin-orbit coupling and
intersystem crossing, a benchmark of the electronic structure for the ground-
state of C3 products, a correlation between interatomic distances and UED
signal for the LIIC pathways, a comparison between the UED PDF signal for a
large number of ground-state geometries and the subset of these geometries
used as initial conditions for the dynamics, an analysis of the hopping
geometries using multidimensional scaling, details about the fitting procedure
for the populations, and a comparison between LZSH and FSSH dynamics. A zip
file is provided with all the critical geometries located (xyz format) as well
as the time-independent basis functions for the $\Delta$PDF and $\Delta I/I$
signal of each (ground-state) photoproduct (represented in Fig. 9).

###### Acknowledgements.

The authors would like to thank CECAM for supporting the organization of the
CECAM workshop ’Triggering out-of-equilibrium dynamics in molecular systems’
in March 2023, which led to the idea of this prediction challenge. This
project has received funding from the European Research Council (ERC) under
the European Union’s Horizon 2020 research and innovation programme (Grant
agreement No. 803718, project SINDAM) and the EPSRC Grants EP/V026690/1. J.J.
and P.S. were supported by the Czech Science Foundation, project number
23-07066S. This work was supported by the grant of Specific university
research – grant No A2_FCHI_2023_048. This article is based upon work from
COST Action CA18212 - Molecular Dynamics in the GAS phase (MD-GAS), supported
by COST (European Cooperation in Science and Technology).

## Data Availability Statement

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

## References

  * Tully (2012) J. C. Tully, J. Chem. Phys. 137, 22A301 (2012).
  * Persico and Granucci (2014) M. Persico and G. Granucci, Theor. Chem. Acc. 133, 1 (2014).
  * Agostini and Curchod (2019) F. Agostini and B. F. E. Curchod, Wiley Interdiscip. Rev. Comput. Mol. Sci. , e1417 (2019).
  * Mai and González (2020) S. Mai and L. González, Angew. Chem. Intl. Ed. 59, 16832 (2020).
  * González and Lindh (2020) L. González and R. Lindh,  _Quantum Chemistry and Dynamics of Excited States: Methods and Applications_ (John Wiley & Sons, 2020).
  * Tully (1990) J. C. Tully, J. Chem. Phys. 93, 1061 (1990).
  * Tully (1998) J. C. Tully, Faraday Discuss. 110, 407 (1998).
  * Martínez, Ben-Nun, and Levine (1996) T. J. Martínez, M. Ben-Nun, and R. D. Levine, J. Phys. Chem. 100, 7884 (1996).
  * Martínez and Levine (1997) T. J. Martínez and R. D. Levine, J. Chem. Soc., Faraday Trans. 93, 941 (1997).
  * Bittner and Rossky (1995) E. R. Bittner and P. J. Rossky, J. Chem. Phys. 103, 8130 (1995).
  * Kapral and Ciccotti (1999) R. Kapral and G. Ciccotti, J. Chem. Phys. 110, 8919 (1999).
  * Donoso and Martens (1998) A. Donoso and C. C. Martens, J. Phys. Chem. A 102, 4291 (1998).
  * Coker and Xiao (1995) D. Coker and L. Xiao, J. Chem. Phys. 102, 496 (1995).
  * Schirmer (1982) J. Schirmer, Phys. Rev. A 26, 2395 (1982).
  * Finley  _et al._ (1998) J. Finley, P.-Å. Malmqvist, B. O. Roos, and L. Serrano-Andrés, Chem. Phys. Lett. 288, 299 (1998).
  * Koch  _et al._ (1995) H. Koch, O. Christiansen, P. Jørgensen, and J. Olsen, Chem. Phys. Lett. 244, 75 (1995).
  * Doltsinis and Marx (2002) N. L. Doltsinis and D. Marx, Phys. Rev. Lett. 88, 166402 (2002).
  * Casida (1995) M. E. Casida, in  _Recent Advances in Density Functional Methods_ , edited by D. P. Chong (Singapore, World Scientific, 1995) p. 155.
  * Petersilka, Gossmann, and Gross (1996) M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys. Rev. Lett. 76, 1212 (1996).
  * González, Escudero, and Serrano-Andrés (2012) L. González, D. Escudero, and L. Serrano-Andrés, ChemPhysChem 13, 28 (2012).
  * Matsika (2021) S. Matsika, Chem. Rev. 121, 9407 (2021).
  * Zewail (2000) A. H. Zewail, J. Phys. Chem. A 104, 5660 (2000).
  * Suchan  _et al._ (2018) J. Suchan, D. Hollas, B. F. E. Curchod, and P. Slavíček, Faraday Discuss. (2018).
  * Pieroni  _et al._ (2023) C. Pieroni, F. Becuzzi, L. Creatini, G. Granucci, and M. Persico, J. Chem. Theory Comput. 19, 2430 (2023).
  * Prlj, Hollas, and Curchod (2023) A. Prlj, D. Hollas, and B. F. E. Curchod, J. Phys. Chem. A 127, 7400 (2023).
  * Crespo-Otero and Barbatti (2018) R. Crespo-Otero and M. Barbatti, Chem. Rev. 118, 7026 (2018).
  * Curchod and Martínez (2018) B. F. E. Curchod and T. J. Martínez, Chem. Rev. 118, 3305 (2018).
  * Subotnik  _et al._ (2016) J. E. Subotnik, A. Jain, B. Landry, A. Petit, W. Ouyang, and N. Bellonzi, Annu. Rev. Phys. Chem. 67, 387 (2016).
  * Hammes-Schiffer (2022) S. Hammes-Schiffer, Philos. Trans. Royal Soc. A 380, 20200377 (2022).
  * Ryabinkin, Joubert-Doriol, and Izmaylov (2017) I. G. Ryabinkin, L. Joubert-Doriol, and A. F. Izmaylov, Acc. Chem. Res. 50, 1785 (2017).
  * Ghosh  _et al._ (2020) S. Ghosh, S. Giannini, K. Lively, and J. Blumberger, Faraday Discuss. 221, 501 (2020).
  * Agostini, Tavernelli, and Ciccotti (2018) F. Agostini, I. Tavernelli, and G. Ciccotti, Eur. Phys. J. B 91, 139 (2018).
  * Santoro  _et al._ (2021) F. Santoro, J. A. Green, L. Martinez-Fernandez, J. Cerezo, and R. Improta, Phys. Chem. Chem. Phys. 23, 8181 (2021).
  * Bondanza  _et al._ (2022) M. Bondanza, B. Demoulin, F. Lipparini, M. Barbatti, and B. Mennucci, J. Phys. Chem. A 126, 6780 (2022).
  * Tavernelli, Curchod, and Rothlisberger (2011) I. Tavernelli, B. F. E. Curchod, and U. Rothlisberger, Chem. Phys. 391, 101 (2011).
  * Janoš and Slavíček (2023) J. Janoš and P. Slavíček, J. Chem. Theory Comput. 19, 8273 (2023).
  * Xia  _et al._ (2015) S.-H. H. Xia, X.-Y. Y. Liu, Q. Fang, and G. Cui, J. Phys. Chem. A 119, 3569 (2015).
  * Kao  _et al._ (2020) M.-H. H. Kao, R. K. Venkatraman, M. N. R. Ashfold, and A. J. Orr-Ewing, Chem. Sci. 11, 1991 (2020).
  * Cui and Fang (2011) G. Cui and W. Fang, J. Phys. Chem. A 115, 1547 (2011).
  * Denschlag and Lee (1967) H. O. Denschlag and E. K. C. Lee, J. Am. Chem. Soc. 89, 4795 (1967).
  * Turro and Southam (1967) N. J. Turro and R. M. Southam, Tetrahedron Lett. 8, 545 (1967).
  * 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. Lee, J. Phys. Chem 80, 1833 (1976).
  * Diau, Kötting, and Zewail (2001) E. W. Diau, C. Kötting, and A. H. Zewail, ChemPhysChem 2, 294 (2001).
  * Liu and Fang (2016) L. Liu and W. H. Fang, J. Chem. Phys. 144 (2016), 10.1063/1.4945782.
  * Trentelman  _et al._ (1990) K. A. Trentelman, D. B. Moss, S. H. Kable, and P. L. Houston, J. Phys. Chem 94, 3031 (1990).
  * 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._ (2013) T. S. Kuhlman, M. Pittelkow, T. I. Sølling, and K. B. Møller, Angew. Chem. Intl. Ed. 52, 2247 (2013).
  * Kuhlman  _et al._ (2012) T. S. Kuhlman, S. P. Sauer, T. I. Sølling, and K. B. Møller, J. Chem. Phys. 137 (2012), 10.1063/1.4742313.
  * Yang  _et al._ (2018) J. Yang, X. Zhu, T. J. Wolf, Z. Li, J. P. F. Nunes, R. Coffee, J. P. Cryan, M. Gühr, K. Hegazy, T. F. Heinz, _et al._ , Science 361, 64 (2018).
  * 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, _et al._ , Nat. Chem. , 1 (2019).
  * Centurion, Wolf, and Yang (2022) M. Centurion, T. J. Wolf, and J. Yang, Annu. Rev. Phys. Chem. 73, 21 (2022).
  * Yang  _et al._ (2020) J. Yang, X. Zhu, J. P. F. Nunes, J. K. Yu, R. M. Parrish, T. J. A. Wolf, M. Centurion, M. Gühr, R. Li, Y. Liu, B. Moore, M. Niebuhr, S. Park, X. Shen, S. Weathersby, T. Weinacht, T. J. Martinez, and X. Wang, Science 368, 885 (2020).
  * Park and Shiozaki (2017) J. W. Park and T. Shiozaki, J. Chem. Theory Comput. 13, 2561 (2017).
  * Shiozaki (2018) T. Shiozaki, Wiley Interdiscip. Rev. Comput. Mol. Sci. 8, e1331 (2018), arXiv:1707.03771 .
  * Frisch  _et al._ (2013) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski, and D. J. Fox, “Gaussian 09 Revision D.01,” (2013).
  * Paul, Myhre, and Koch (2021) A. C. Paul, R. H. Myhre, and H. Koch, J. Chem. Theory Comput. 17, 117 (2021).
  * Folkestad  _et al._ (2020) S. D. Folkestad, E. F. Kjønstad, R. H. Myhre, J. H. Andersen, A. Balbi, S. Coriani, T. Giovannini, L. Goletto, T. S. Haugland, A. Hutcheson, I.-M. Høyvik, T. Moitra, A. C. Paul, M. Scavino, A. S. Skeidsvoll, Å. H. Tveten, and H. Koch, J. Chem. Phys. 152 (2020).
  * Marsili, Prlj, and Curchod (2021) E. Marsili, A. Prlj, and B. F. E. Curchod, Phys. Chem. Chem. Phys. 23, 12945 (2021).
  * Werner  _et al._ (2012) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, Wiley Interdiscip. Rev. Comput. Mol. Sci. 2, 242 (2012).
  * Crespo-Otero and Barbatti (2012) R. Crespo-Otero and M. Barbatti, Theor. Chem. Acc. 131, 1 (2012).
  * Hollas  _et al._ (2019) D. Hollas, J. Suchan, M. Ončák, and P. Slavíček, “PHOTOX/ABIN: Pre-release of version 1.1, https://github.com/ PHOTOX/ABIN,” (2019).
  * Ceriotti, Bussi, and Parrinello (2009) M. Ceriotti, G. Bussi, and M. Parrinello, Phys. Rev. Lett. 103, 030603 (2009).
  * Ceriotti, Bussi, and Parrinello (2010) M. Ceriotti, G. Bussi, and M. Parrinello, J. Chem. Theory Comput. 6, 1170 (2010), arXiv:1204.0822 .
  * Prlj  _et al._ (2022) A. Prlj, E. Marsili, L. Hutton, D. Hollas, D. Shchepanovska, D. R. Glowacki, P. Slavíček, and B. F. E. Curchod, ACS Earth Space Chem. 6, 207 (2022).
  * Sršeň, Hollas, and Slavíček (2018) Š. Sršeň, D. Hollas, and P. Slavíček, Phys. Chem. Chem. Phys. 20, 6421 (2018).
  * Silverman (1986) B. Silverman,  _Density Estimation for Statistics and Data Analysis_ (Taylor & Francis, 1986).
  * Suchan, Janoš, and Slavíček (2020) J. Suchan, J. Janoš, and P. Slavíček, J. Chem. Theory Comput. 16, 5809 (2020).
  * Martínez-Mesa and Saalfrank (2015) A. Martínez-Mesa and P. Saalfrank, J. Chem. Phys. 142, 194107 (2015).
  * Granucci and Persico (2007) G. Granucci and M. Persico, J. Chem. Phys. 126, 134114 (2007).
  * Mignolet, Curchod, and Martínez (2016) B. Mignolet, B. F. E. Curchod, and T. J. Martínez, Angew. Chem. 128, 15217 (2016).
  * Pathak  _et al._ (2020) S. Pathak, L. M. Ibele, R. Boll, C. Callegari, A. Demidovich, B. Erk, R. Feifel, R. Forbes, M. Di Fraia, L. Giannessi, _et al._ , Nat. Chem. 12, 795 (2020).
  * (73) “Prediction challenge: Cyclobutanone photochemistry,” https://publishing.aip.org/publications/journals/special-topics/jcp/prediction-challenge-cyclobutanone-photochemistry/, accessed: 2023-11-28.
  * McMorrow and Als-Nielsen (2011) D. McMorrow and J. Als-Nielsen,  _Elements of modern X-ray physics_ (John Wiley & Sons, 2011).
  * Parrish and Martínez (2019) R. M. Parrish and T. J. Martínez, J. Chem. Theory Comput. 15, 1523 (2019).
  * Figueira Nunes  _et al._ (2024) J. P. Figueira Nunes, L. M. Ibele, S. Pathak, A. R. Attar, S. Bhattacharyya, R. Boll, K. Borne, M. Centurion, B. Erk, M.-F. Lin, R. J. G. Forbes, N. Goff, C. S. Hansen, M. Hoffmann, D. M. P. Holland, R. A. Ingle, D. Luo, S. B. Muvva, A. H. Reid, A. Rouzée, A. Rudenko, S. K. Saha, X. Shen, A. S. Venkatachalam, X. Wang, M. R. Ware, S. P. Weathersby, K. Wilkin, T. J. A. Wolf, Y. Xiong, J. Yang, M. N. R. Ashfold, D. Rolles, and B. F. E. Curchod, J. Am. Chem. Soc. 0, 0 (2024).
  * Lannelongue, Grealey, and Inouye (2021) L. Lannelongue, J. Grealey, and M. Inouye, Adv. Sci. 8, 2100707 (2021).
  * Borgers and Strauss (1966) T. R. Borgers and H. L. Strauss, J. Chem. Phys. 45, 947 (1966).
  * Scharpen and Laurie (1968) L. R. H. Scharpen and V. W. Laurie, J. Chem. Phys. 49, 326 (1968).
  * Alonso  _et al._ (1992) J. L. Alonso, R. Spiehl, A. Guarnieri, J. C. López, and A. G. Lesarri, J. Mol. Spectrosc. 156, 341 (1992).
  * Udvarhazi and El-Sayed (1965) A. Udvarhazi and M. A. El-Sayed, J. Chem. Phys. 42, 3335 (1965).
  * Keller-Rudek  _et al._ (2013) H. Keller-Rudek, G. K. Moortgat, R. Sander, and R. Sörensen, Earth Syst. Sci. Data 5, 365 (2013).
  * Favero, Granucci, and Persico (2013) L. Favero, G. Granucci, and M. Persico, Phys. Chem. Chem. Phys. 15, 20651 (2013).
  * Shen  _et al._ (2019) X. Shen, J. P. F. Nunes, J. Yang, R. K. Jobe, R. K. Li, M.-F. Lin, B. Moore, M. Niebuhr, S. P. Weathersby, T. J. A. Wolf, C. Yoneda, M. Guehr, M. Centurion, and X. J. Wang, Struct. Dyn. 6, 054305 (2019).
  * Champenois  _et al._ (2023) E. G. Champenois, N. H. List, M. Ware, M. Britton, P. H. Bucksbaum, X. Cheng, M. Centurion, J. P. Cryan, R. Forbes, I. Gabalski, K. Hegazy, M. C. Hoffmann, A. J. Howard, F. Ji, M.-F. Lin, J. P. F. Nunes, X. Shen, J. Yang, X. Wang, T. J. Martinez, and T. J. A. Wolf, Phys. Rev. Lett. 131, 143001 (2023).