The photochemistry of Rydberg-excited cyclobutanone: Photoinduced processes and ground state dynamics Open Access

Ketones, i.e., organic compounds containing a carbonyl group, are among the simplest chromophores owing to their small size and low density of valence electronically excited states. Consequently, their photochemistry has been studied for several decades.^1,2 Cyclic ketones form an important subclass of these systems: Despite their apparently equivalent simplicity, their photochemistry is often comparatively complex and features activity across multiple competing photochemical channels most commonly characterized as Norrish Type I processes. In such a process, and upon excitation of an electron from the nonbonding orbital (n) on the carbonyl oxygen atom to the antibonding (π*) molecular orbital of the carbonyl group, a carbon–carbon bond adjacent to the carbonyl group (i.e., C_α; Fig. 1) cleaves, opening up the possibility for formation of a variety of fragmentory products on the electronic ground state surface after non-radiative relaxation. As an exemplar cyclic ketone, cyclobutanone has particular photochemistry owing to the high degree of ring strain in the cyclobutane ring (arising as a consequence of the small ring size).³ The particular photochemistry of cyclobutanone has attracted significant interest from both an experimental^3–16 and theoretical^17–24 perspective.

Many of these previous studies on cyclobutanone have aimed at elucidating the photochemistry and, in particular, the excited state dynamics arising from excitation at the weak S₁ ←S₀ absorption band (∼330–240 nm) associated with the symmetry-forbidden π* ← n transition. Early ultrafast spectroscopy of cyclobutanone in the gas phase carried out by the authors of the work of Diau et al.¹⁰ revealed that, upon excitation at 307 nm (i.e., slightly above energy of the π* ← n transition by ∼2 kcal mol⁻¹) α-cleavage occurs on a timescale of ∼5 ps, driven initially by the (i) cyclobutane ring puckering, (ii) carbonyl out-of-plane wagging, and (iii) carbonyl stretching modes. Following along the α-cleavage channel, an S₁/S₀ conical intersection (CI) is encountered via which non-radiative relaxation to the electronic ground state can proceed; this channel ultimately leads to the production of either (i) a vibrationally hot cyclobutanone on the S₀ state, (ii) diradical intermediates that fragment to yield C₂H₄ + CH₂CO (C2 products), or (iii) C₃H₆ + CO (C3 products). The C3:C2 fragmentation quantum yields are known to exhibit a strong wavelength dependence;^4,7,8,13,25 the C3:C2 product ratio is reported to be 0.5 on excitation at 313 nm and to increase to 0.8 at 248 nm and 1.0 at 200 nm.^13,26 At wavelengths longer than 315 nm, there is a marked increase in the C3:C2 product ratio; it is reported to be 2.0 at 326 nm and as high as 7.0 at 344 nm.^7,26 The increase is indicative of an alternative dormant mechanism in the S₁ state, which becomes operational at longer wavelengths; this alternative mechanism is proposed to involve photochemistry on the triplet manifold²⁶ as there is insufficient energy available at longer wavelengths to overcome the barrier encountered along the α-cleavage channel and access the S₁/S₀ CI. This mechanistic picture of cyclobutanone in the gas phase is largely maintained in solution as evidenced in the work of Kao et al.³ Norrish Type I α-cleavage is known to occur on the sub-picosecond timescale in solution. The involvement of photochemistry on the triplet manifold has been invoked in Norrish-type mechanisms,¹³ and the authors of the work of Kao et al.³ concede that—even with sufficient excess energy available to overcome the barrier encountered along the α-cleavage channel—triplet states may still play a mechanistic role, although this is likely limited to indirect channels active on timescales longer than 500 ps.³ 

The second absorption band of cyclobutanone is located between ∼206 and 182 nm and is assigned to the 3s ← n transition into the second electronically excited singlet state (S₂): a state that exhibits Rydberg character.¹⁴ The authors of the work of Trentelman et al.¹⁵ have studied the 193 nm photolysis of cyclobutanone, reporting that 57% of the electronically excited state species form the C3 products, 30% form the C2 products (on the electronically excited state), and 13% form the C2 products (on the hot electronic ground state). The observations presented in the work of Trentelman et al.¹⁵ suggest that the formation of C3 products is potentially a slow process that requires intersystem crossing (ISC) onto the triplet manifold. However, the proposed photochemical mechanism assumes an S₁ intermediate that survives long enough to undergo thermalization and, therefore, produce a statistical partitioning of products; this assumption does not necessarily hold given the speed at which S₀ ← S₁ relaxation is expected to take place. The authors of the work of Kuhlman et al.^11,27 have studied the 200 nm photolysis of cyclobutanone, comparatively, associating S₁ ← S₂ relaxation with two time constants of ∼350 and ∼750 fs. Using photoelectron spectroscopy and time-resolved mass spectrometry (TR-MS), the authors were able to resolve two species: the parent ion and a fragment ion with a mass-to-charge ratio m/z = 42 (corresponding to H₂C–C=O).^11,27 Both the species exhibited similar electronically excited state decay constants.^11,27 The authors also reported a blue shift in the photoelectron spectrum arising from a coherent oscillatory motion assigned to a low-frequency (35 cm⁻¹) cyclobutane ring puckering mode activated by removal of an electron from the carbonyl oxygen n orbital;^11,27 this process leads to a relaxation of the sp² hybridization of the carbonyl carbon due to increased mixing with components of the bonding orbitals in the carbonyl group as well as a relaxation of the adjacent carbon–carbon bonds.

From the perspective of theory, the authors of the work of Xia et al.¹⁷ carried out a quantum chemical investigation in order to locate key minima, transition states, and CIs for cyclobutanone, as well as to develop (relaxed) two-dimensional singlet (S₁) and triplet (T₁) excited state potential energy surfaces using high-level quantum chemical calculations. On the basis of these calculations, the authors of the work of Xia et al. proposed that ring opening occurs on the S₁ state through an accessible S₁/S₀ CI.¹⁷ The authors of the work of Liu et al.¹⁸ extended this work by performing excited state dynamics simulations under the ab initio multiple spawning (AIMS)^28,29 framework following simulated S₁ ← S₀ photoexcitation. The authors corroborated that S₀ ← S₁ relaxation occurs primarily through the S₁/S₀ CI located along the α-cleavage channel; relaxation was observed at geometries where the C_β ≈ 1.6 Å and C_α ≈ 2.5–3.5 Å.¹⁸ Comparatively, relaxation through CI associated with the formation of C2 products was observed at geometries where C_β ≈ 4.0 Å, although these geometries comprised as little as ∼15% of geometries at which relaxation was observed.¹⁸ The authors reported an excited state lifetime on the S₁ state of ∼500 fs and, on this basis, proposed that non-ergodic behavior was the driver of the change in C3:C2 product ratio as a function of excitation wavelength.¹⁸ In adopting an alternative approach, the authors of the work of Kuhlman et al.²⁴ developed a four-state, five-dimensional model Hamiltonian and carried out quantum dynamics to study the S₁ ← S₂ relaxation in cyclobutanone. The authors concluded that S₁ ← S₂ relaxation is driven by specific nuclear modes (including the cyclobutane ring puckering, carbonyl out-of-plane deformation, symmetric and asymmetric C–CO–C stretches, and carbonyl stretching modes) that couple the S₂ and S₁ states and promote population transfer on the sub-picosecond timescale.²⁴ While the model presented in the work of Kuhlman et al.²⁴ is informative, their potentials were built within a harmonic normal mode representation and are consequently unable to address satisfactorily large amplitude nuclear motions associated with the formation of fragmentary products.

Despite the concerted efforts of experiment and theory, there remain a significant number of open questions regarding the photochemistry and photochemical dynamics of cyclobutanone that evolve post-photoexcitation to the S₂ Rydberg (3s ← n) state. Nonetheless, the emergence of modern light sources^30,31 is enabling the study of ultrafast excited state dynamics using both x rays and electrons^32–35 with ever-increasing spatial and temporal resolution; experiments like this are providing the potential for novel and increasingly detailed insights into complex photochemical processes such as these. These developments are driving commensurate progress in theory and are starting to bring into focus a crucial question (which is that posed in the present Special Issue to which this article contributes): How accurate are modern excited state dynamics simulations really?

In this article, we combine quantum multiconfigurational time-dependent Hartree (MCTDH) and excited state trajectory surface-hopping dynamics (TSHD) simulations at the density functional [linear-response time-dependent density functional theory (LR-TDDFT)(PBE0)] and algebraic diagrammatic constructor [ADC(2)] levels of theory to explore the relaxation dynamics of cyclobutanone post-photoexcitation to the S₂ Rydberg (3s ← n) state. Our MCTDH and TSHD simulations are subsequently used to predict the scattering signals in an ultrafast electron diffraction (UED) experiment in reciprocal and real space to establish how the complex relaxation dynamics of cyclobutanone might manifest as an experimental observable. In such a UED experiment, free cyclobutanone molecules introduced in vacuo will be photoexcited at 200 nm [i.e., into the S₂ Rydberg (3s ← n) state] and probed via electron scattering at a sequence of temporal delays to image directly the evolving excited state structural dynamics.

All (LR-)TDDFT and ADC(2) quantum chemical calculations were carried out using TURBOMOLE (v7.4).^36,37 The density functional calculations on the electronic ground and excited states used density functional theory (DFT) and linear-response time-dependent DFT (LR-TDDFT), respectively, with the PBE0³⁸ density functional and the aug-cc-pVDZ^39,40 basis set. The Tamm–Dancoff⁴¹ approximation (TDA) was used throughout. Benchmarks of the level of electronic structure theory and the basis set are available in the supplementary material. Vibrational analysis were performed to confirm the absence of imaginary frequencies at the ground state minimum. All minimum-energy conical intersections (MECIs) were optimized using Turbomole (v7.4) coupled with an external (penalty-function-based) optimizer.⁴² Linear interpolation in internal coordinates (LIICs) plots of the potential energy surfaces between critical geometries and these MECI are available in the supplementary material.

Additional quantum chemical calculations were carried out at critical geometries using the second-order algebraic diagrammatic construction scheme [ADC(2)] and n-electron valence state perturbation theory (NEVPT2). All ADC(2) calculations were carried out using Turbomole (v7.4);^36,37 all NEVPT2 calculations were carried out using ORCA.^43,44

All expansion coefficients were obtained using the in-house-developed VCMaker^50,51 software, available at Ref. 52.

The quantum dynamics simulations over the multidimensional potential energy surface(s) were carried out using the MCTDH^53,54 approach as implemented in Quantics.⁵⁵ The initial wavefunction for the electronic ground state was built using one-dimensional harmonic oscillator functions with an expectation value of zero momentum along all coordinates; this wavefunction was then projected vertically from the electronic ground state (S₀) into the S₂ Rydberg (3s ← n) state at the Franck–Condon geometry (Q = 0). The complete computational details are available in the supplementary material.

“On-the-fly” excited state TSHD simulations were carried out using Newton-X (v2.4).^56,57 All potential energy surfaces, derivatives, and respective couplings were computed using Turbomole (v7.4)^36,37 at two separate levels of theory [LR-TDDFT(PBE0) and ADC(2)] in two separate sets of simulations. Trajectories were propagated using the velocity Verlet algorithm^58,59 for a maximum time of 5 ps (t_max = 5000 fs) with a time step of 0.5 fs (dt = 0.5 fs). The state-to-state transitions (or “hops”) were simulated using the Hammes-Schiffer-and-Tully fewest-switches algorithm⁶⁰ with the state-to-state couplings estimated via the time-dependent Baeck–An⁶¹ scheme. The excited state TSHD simulations were initiated by vertical projection of a set of initial conditions, generated according to a Wigner distribution with a temperature of 100 K, from the electronic ground state (S₀) into the S₂ Rydberg (3s ← n) state.

To avoid instabilities at/around the S₁/S₀ crossing seam, S₀ ← S₁ surface hops were forced when the S₁/S₀ energy gap fulfilled the criterion $ΔES1−S0<0.1$ eV.

The focus of this work is upon the ultrafast excited state dynamics of cyclobutanone and, accordingly, the ultrafast electron diffraction scattering signal is presented as a transient (i.e., as an excited state − ground state “difference”) signal where appropriate under the nuclear ensemble model. The ground state signal used to generate the transient is predicted from the trajectory surface-hopping dynamics initial conditions, i.e., the nuclear ensemble representing the state of the system at t = 0. No scaling for, e.g., excitation percentage/photolysis yield has been carried out.

Table I shows the energies of the low-lying singlet (S_n; n = [1,…,5]) and triplet (T_n; n = [1,…,5]) electronically excited states of cyclobutanone. Figure S1 shows the molecular orbitals involved in the electronic transitions. The symmetry point group of cyclobutanone (C_s at the S₀ minimum-energy geometry) is used to characterize the symmetries of the electronically excited states (an important facet of the generation and interpretation of the model Hamiltonian used in Sec. III B).

At the S₀ minimum-energy geometry, the S₁ and S₂ excited states are located (1) 4.21 and 5.99 eV, respectively, above the electronic ground state at the ADC(2)/aug-cc-pVTZ level of theory and (2) 4.33 and 6.10 eV, respectively, above the electronic ground state at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory. Both levels of theory give good agreement with the experimental absorption spectrum recorded in the work of Diau et al.¹⁰ in which absorption bands are observed at ∼4.2 and ∼6.2 eV. At both levels of theory, the S₁ and S₂ states have π* ← n (LUMO ← HOMO) and 3s ← n (LUMO + 1 ← HOMO) character, respectively. The Rydberg character of the S₂ state requires a larger basis set to describe accurately its energy and electronic structure, while—in contrast—the valence S₁ state exhibits little to no dependence on basis set (Table S1). Both the S₂ and S₁ states are of A″ symmetry; direct first-order vibronic coupling between these two states is consequently forbidden, and this can be expected to slow down the rate of internal conversion (IC) between the two states.

The ADC(2) calculations give good agreement with the experimentally observed ordering and energies of the electronic states at the Franck–Condon geometry. However, previous work⁶² has highlighted a limitation of the ADC(2) approach for studying non-radiative relaxation pathways of carbonyl-containing molecules as it has a tendency to predict artificial S₁/S₀ crossings along the carbonyl stretching coordinate. These artificial S₁/S₀ crossings arise from (i) a π* ← n potential that is too shallow combined with (ii) a ground state potential that destabilizes too rapidly. To demonstrate this, we present calculations of the two potentials along the carbonyl stretching mode in Fig. S2 at the LR-TDDFT(PBE0), ADC(2), and NEVPT2(12,12) levels of theory. Figure S2 indicates good agreement on the form of the potentials between the LR-TDDFT(PBE0) and NEVPT2(12,12) levels of theory: Although the ground state and π* ← n surfaces come closer in energy, they do not cross. In contrast, and as expected from Ref. 62, there is a clear (and artificial) S₁/S₀ crossing between the two potentials at the ADC(2) level of theory occurring at a C=O bond length of ≈1.6 Å. This low-lying (i.e., accessible) S₁/S₀ crossing point is likely to distort the excited state dynamics crossing from the S₁ to ground state surface and indeed such dynamics are presented in the supplementary material.

The S₁-state minimum-energy geometry is reached from the Franck–Condon point via out-of-plane puckering of the cyclobutane ring and a slight elongation of the C=O bond (1.21–1.26 Å); the carbon–carbon bonds in the cyclobutane ring (C_α and C_β) remain almost unchanged. These structural changes destabilize the electronic ground state (Table S3); it is increased in energy by 1.13 eV relative to the S₀-state minimum-energy geometry at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory. In contrast, the S₁ (π* ← n) state is stabilized, decreasing the energy gap with the S₀ state while increasing the gap with the higher-lying singlet states (S_n; n > 1). The predicted emission energy from the S₁ minimum-energy geometry is ∼3.9 eV, a value that is in excellent agreement with the emission spectrum recorded in the work of Lee et al.⁶³ The broad band observed in the absorption spectrum recorded in the work of Diau et al.¹⁰ is consistent with a short-lived electronically excited state, however, which is indicative of competitive photochemical channels, e.g., non-radiative decay on the femto/picosecond timescale through accessible CIs as discussed by Liu and Fang.¹⁸ 

The S₂-state minimum-energy geometry (Table S4), in contrast, is reached from the Franck–Condon point via contraction of the C=O bond (1.21–1.16 Å) and slight elongation of the carbon–carbon bonds in the cyclobutane ring (C_α and C_β). The effect of these structural changes is a destabilization of the ground and S₁ states; the S₀ state increases in energy by 0.51 eV with a similar (0.58 eV) increase observed for the S₁ state at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory. The absorption spectrum for excitation into this state exhibits a distinct vibronic structure suggestive of both a longer-lived electronically excited state and the activation of vibrational modes on excitation.¹⁰ 

Previous studies, e.g., that of Kao et al.,³ have hypothesized as to the potential role of triplet states in the photochemistry of cyclobutanone; this is our motivation for presenting the details of these states (T_n; n = [1,…,5]) in Table I and the relevant state-to-state spin–orbit couplings (SOCs) in the supplementary material (Table S7). At the Franck–Condon geometry, the lowest-energy electronically excited triplet state (T₁) is located at 3.73 eV above the electronic ground state at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory and is of the same π* ← n character as the S₁ state; consequently, S₁/T₁ SOC will be formally forbidden.^64,65 The T₂ and T₃ states are near-degenerate and located a little under the S₂ (6.10 eV) at 5.98 and 6.02 eV, respectively, above the electronic ground state at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory. The SOC between these triplet states and the low-lying singlet electronically excited states is generally small (Table S7), suggesting that the triplet states are unlikely to play a significant role in the early-time (e.g., $<$2 ps) channels.

Figure 2(b) shows the structure of an S₂/S₁ MECI, while Figs. 2(c)–2(e) show the structures of three S₁/S₀ MECI as located via penalty-function-constrained optimization at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory. Cartesian coordinates are given in the supplementary material. Potential energy surfaces between the Franck–Condon point and each of the S₁/S₀ MECI were calculated via linear interpolation in internal coordinates (LIICs) at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory, and these are also given in the supplementary material.

The S₂/S₁ MECI [Fig. 2(b)] is located 5.85 eV above the S₀ minimum-energy geometry, i.e., ∼0.5 eV below the 3s ← n excitation energy, at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory. Its structure is similar to the Franck–Condon geometry although it features a symmetry-breaking distortion of the cyclobutane ring [seen clearly on inspection of the position of the hydrogen atoms in Fig. 2(b)]. The potential energy surface between the Franck–Condon geometry and this S₂/S₁ MECI is barrierless at LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory (Fig. S3), which would suggest that S₂/S₁ internal conversion should be (ultra)fast in the absence of any additional considerations. However, the symmetries of the two states at the Franck–Condon geometry are such that the interstate coupling is forbidden; even at the distorted S₂/S₁ MECI geometry, it is weak and results predictably in a longer lifetime for the S₂ state.

The three S₁/S₀ MECIs [Figs. 2(c)–2(e)] are located in qualitative agreement (geometrically and energetically) with those reported by Liu and Fang¹⁸ [located by the authors at the complete active space self-consistent field (CASSCF) level of theory]. The first [Fig. 2(c)] one is located along the α-cleavage channel, the second [Fig. 2(d)] along the β-cleavage channel, and the third [Fig. 2(e)] along a concerted α/β-cleavage channel. The energetic ordering of the three S₁/S₀ MECIs (3.40, 2.56, and 5.0 eV, respectively, above the S₀ minimum-energy geometry at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory) follows qualitatively the trend observed for the three S₁/S₀ MECIs in Ref. 18, although it is important to note that the single-reference nature of LR-TDDFT(PBE0) favors charged rather than biradical dissociation along the α- and β-cleavage channels and, furthermore, renders it unable to describe properly the topology/dimensionality of the S₁/S₀ crossing seam. Although all of the S₁/S₀ MECIs are energetically accessible, i.e., they and their barriers are submerged relative to the ∼6.2 eV excitation energy, the description of the potential energy surface around the S₁/S₀ MECIs is likely to be problematic for LR-TDDFT(PBE0) and it is quite possible that this might affect product production by influencing the dynamics through and in the vicinity of the MECI and, subsequently, the motion of the trajectory/wavepacket on the electronic ground state potential energy surface.

Overcoming these aforementioned limitations could be achieved using a multireference (active space) method, e.g., CASSCF/CASPT2 or NEVPT2, for the excited state TSHD simulations. However, the performance of this family of methods is greatly dependent on the choice of active space; an appropriate active space should be large enough to incorporate all of the orbitals required over all of the nuclear configurations explored in the TSHD simulations while not too large so as to render the TSHD simulations computationally costly to the point of intractability. We found an active space of eight electrons in eight orbitals [e.g., NEVPT2(8,8)] unstable with respect to orbital rotation(s) at some geometries (particularly those where the cyclobutane ring was distorted), while a larger active space of 12 electrons in 12 orbitals [e.g., NEVPT2(12,12)] was too computationally expensive to carry out practicably TSHD simulations. Consequently, we elected to carry out our excited state TSHD simulations at the LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory, keeping in mind the aforementioned (although well-understood) limitations and their potential impact on the dynamics (which we discuss in detail). However, we note within the context of the present challenge that other contributors have performed excited state dynamics simulations with multireference (active space) methods, e.g., CASSCF based on an eight-electron-in-eleven-orbital active space⁶⁶ and extended multistate CASPT2 (XMS-CASPT2) based on an eight-electron-in-eight-orbital active space.⁶⁷ 

To (i) identify possible photochemical (fragmentation) channels, (ii) clarify the potential involvement of triplet states, and (iii) establish a mechanism for internal conversion from the initially excited S₂ (3s ← n) state to the lowest-energy singlet electronically excited state (S₁) at early times, we developed a model Hamiltonian and carried out MCTDH simulations. The model Hamiltonian comprised the electronic ground state (S₀) and nine electronically excited states [four singlets (S_n; n = [1,…,4]) and five triplets (T_n; n = [1,…,5])]. The inclusion of singlet states higher in energy than the S₂ (e.g., S_n; n > 2) is motivated by the absence of vibronic coupling between the S₁ and S₂ states, both of which are of A″ symmetry (Table I); here, coupling to higher-lying singlet states of, e.g., alternative symmetries offers the potential of second-order population transfer channels comparable to those identified in other systems where direct population transfer channels are weak or otherwise absent.^68–70 

The model Hamiltonian incorporated eight degrees of vibrational freedom in nuclear configurational space: ν₁, ν₇, ν₁₀, ν₁₁, ν₁₂, ν₁₃, ν₁₅, and ν₂₁, which were selected on the basis of the magnitude of their first-order couplings and the symmetries of the vibrational modes. The degrees of vibrational freedom included in the model Hamiltonian comprised cyclobutane ring puckering (ν₁ and ν₁₂), symmetric and antisymmetric cyclobutane ring breathing (ν₇ and ν₁₀), cyclobutane ring deformation (ν₁₁, ν₁₃ and ν₁₅), and the carbonyl stretching mode (ν₂₁), for consistency with previous work, e.g., that of Kuhlman et al.¹¹ To assess the approximate accuracy of our model Hamiltonian, we present Fig. 3, which compares the experimental¹⁰ and theoretical absorption spectrum for the S₂ state. There is excellent agreement between the two absorption spectra, supporting the assertion that our model Hamiltonian describes accurately the (local) potential energy surface(s). The absorption spectra in Fig. 3 feature a vibrational progression with structured peaks separated by ∼0.1 eV (consistent with the carbonyl out-of-plane wagging and cyclobutane ring puckering vibrational modes previously identified⁷¹ and arising, in the present model, via the inclusion of ν₁₂ into the model Hamiltonian).

Figure 4 shows the transfer of population from the S₂ state over 600 fs post-photoexcitation into the S₂ (3s ← n) state as computed via quantum dynamics. The S₁ ← S₂ population transfer occurs at a rate of ∼1.67 × 10⁻¹² s⁻¹ with the S₁ state population exceeding S₂ within ∼350 fs. This is qualitatively consistent with the decay of the peak associated with the S₂ state in the photoelectron spectrum recorded in the work of Kuhlman et al.¹¹ (the authors report a biexponential fit yielding time constants of ∼350 and 750 fs). The S₁ ← S₂ population transfer in the model Hamiltonian obtains intensity via mixing with the higher-lying singlet electronically excited states, especially the S₄ state (the lowest-lying singlet electronically excited state of A′ character). This occurs most prominently along ν₁₁. While this mode is responsible for the state-to-state coupling, it is ν₁, ν₁₂, ν₁₃, and ν₂₁ that exhibit the largest-amplitude nuclear oscillations and electronically excited state structural changes that are likely to be observed in the ultrafast electron diffraction experiment. The population oscillations observed between the S₂ and S₁ states are associated with the overcoherence of the reduced model Hamiltonian and are, in any case, faster than the temporal resolution of the ultrafast electron diffraction experiments.

S₀ ← S₁ population transfer is not included in the present model Hamiltonian. The high energy of the S₂ state results in an exceptionally hot wavefunction on the electronic ground state that is difficult to converge under the framework of the quantum dynamics simulations. In addition, previous works—and our own quantum chemical calculations—have identified S₁/S₀ MECI at highly distorted geometries beyond the limits of the normal model representation on which the quantum dynamics simulations are predicated (the representation is only valid to small distortions from the equilibrium, e.g., S₀ minimum-energy, geometry). Consequently, to describe more completely the excited state relaxation dynamics, we explore excited state molecular dynamics simulations operating in unconstrained nuclear configurational space through the trajectory surface-hopping approach.

Figure 5 shows the populations of the S₂, S₁, and S₀ states over the first 2 ps post-photoexcitation into the S₂ (3s ← n) state as obtained from 289 on-the-fly surface-hopping trajectories propagated at LR-TDDFT(PBE0)/aug-cc-pVDZ level of theory. A similar set of trajectory surface-hopping dynamics propagated at the ADC(2)/aug-cc-pVDZ level of theory are presented and analyzed in the supplementary material and are of comparative interest. Figure 5 shows S₁ ← S₂ population transfer with a decay constant of ∼356 fs. This is slightly faster than observed in the quantum dynamics simulations, i.e., 50% of population decays from S₂ in 225 fs, which is due to the inclusion of the ground state to which the wavepacket can rapidly decay. This decay constant is in close agreement with the fastest time-constant extracted from a photoelectron spectroscopic study,¹¹ but we do not see any dynamics associated with the ∼700 fs component reported. Interestingly, this slower component is in close agreement with the population kinetics for the dynamics performed using potentials calculated at ADC(2) level of theory, shown in the supplementary material.

While the timescales between the quantum dynamics and TSH simulations suggest similar dynamics, further analysis is required to assess this in more detail. To achieve this, we transform the first 500 fs of excited state molecular dynamics from Cartesian coordinates into a normal mode representation similar to Ref. 72. Figure 6 shows the normal modes active during the TSH simulations as a bar chart of the number of times each normal mode features in the top 6 of the largest displacements during a trajectory’s excited state dynamics. The red bars represent those normal modes included in the model Hamiltonian. This shows close agreement between the modes that appear most frequently in the TSH simulations and those included in the quantum dynamics. The most notable exception is ν₅, which exhibits larger amplitude motion due to the flat nature of the potential but does not act as either tuning or coupling mode and therefore does not strongly influence the excited state dynamics.

Figure 7 shows average C_α and C_β bond lengths for the structures where each trajectory hops from the S₂–S₁ (a) and S₁–S₀ (b) states. For the former, there is a clear cluster around 1.5–1.6 Å consistent with the ground state structure and therefore close to the Franck–Condon geometry, as expected from the optimized S₂/S₁ discussed above. In contrast, for the S₁–S₀ hopping geometries, there is a significant change, with the majority of hops occurring for C_α bond lengths $>$2.2 Å. Using the geometries provided in the supplementary material, the C_α CI occurs when the C_α bond length is 2.35 Å, with very little corresponding change along the C_β bond. This suggests, in agreement with previous work, that crossing from the S₁–S₀ occurs primarily at the CI exhibiting a C_α bond break.

Figure 8 shows the fractional population of the photoproducts formed from the TSH trajectories. This indicates ∼20% of the excited states form the C₂H₄ + CH₂CO, i.e., decay via the C2 channel, while 2.5% forms the C3 products. The ring-open species are formed, but are very short-lived and either contribute forming either the C2 or C3 products or undergo bond reformation to form vibrationally excited ground state cyclobutanone. The formation of C2 is comparable to but lower than other excited state dynamics simulations performed at a higher level of theory reported in Ref. 67 (34%) and by the authors of the work of Trentelman et al.¹⁵ who reported 43% of yield experimentally. The major discrepancy in our simulations occurs for the C3 channels, which is $>$60% in these previous works. The near-absence of the C3 channel is associated with the multireference character of the potential in this region and the bias of single reference methods for charged rather than biradical bond breaking. While this does not significantly increase the energies of the C_α and C_β CIs (see Fig. S4), it does increase the energy of the double bond breaking CI, making the formation of the C3 channel challenging. To assess this, we also perform dynamics using trajectories in the T₁ state, performed using unrestricted Kohn–Sham facilitating the description of biradical character. These were initiated at random from trajectories populating the S₁ state. Importantly, these show a much higher formation of C3 photoproducts (C3: 53%, C2: 5% and ring-open: 17%) consistent with previous experiments.¹⁵ 

Figure 9 shows the time-resolved electron diffraction simulations arising after photoexcitation of cyclobutanone into the S₂ state. Figure 9(a) shows the electron diffraction scattering signal as calculated, while Fig. 9(b) is convolved along the temporal axis with a Gaussian kernel (FWHM = 150 fs) to reproduce the effect(s) of the finite temporal resolution of the proposed electron diffraction experiment. Figures 9(c) and 9(d) show time-resolved pair-distribution function (PDF) maps with and without temporal broadening, respectively, and were produced via sine transformation of the modified scattering intensity maps in Figs. 9(a) and 9(b), respectively.

The modified scattering intensity maps [Figs. 9(a) and 9(b)] show two strong negative (∼1 and 9 Å⁻¹) and two positive (∼2.7 and 7.5 Å⁻¹) features but do not reveal the richness of the dynamics that reflect the complex photochemistry of cyclobutanone, in part due to the incoherent/stochastic nature of the photochemical processes taking place.

A deeper understanding of the structural changes can be established from Figs. 9(c) and 9(d), which show the time-resolved PDF. For clarity, Fig. 10 shows the PDF(t = 0 fs), PDF(t = 2000 fs), and the ΔPDF(t = 2000 fs).The PDF acquired for the 289 initial conditions exhibits three peaks at ∼1.5, ∼2.5, and ∼3.0 Å. The first corresponds to the neighboring C–C and C=O distances, the second corresponds to C–C distance on the opposite side of the cycle, and the final peak corresponds to the C–O distances, which are not directly bonded. The transient indicates primarily a loss of the first two peaks associated with dissociation and the formation of the C2 products. We note here that despite the aforementioned differences in the photoproduct formation, the transient scattering and PDFs are very similar to those in Ref. 67. This highlights the challenge in disentangling the exact photoproduct formation of these systems due to the overlapping bands.

In this work, we have carried out quantum and excited state trajectory surface-hopping molecular dynamics simulations to study the electronically excited state relaxation mechanisms and electronic ground state dynamics of cyclobutanone post-photoexcitation into the S₂ Rydberg (3s ← n) state. Our focus has been upon translating these simulations to predict the experimental observables associated with the ultrafast electron diffraction experiments and ultimately to answer the question: Can excited state dynamics simulations be predictive? However, even for small molecules such as cyclobutanone, certain approximations in the underlying computational chemical methods are required, which will influence the outcome of such simulations: We highlight this in comparison between the present work and other works related to the same challenge.^{66,67,73–83} Consequently, in this section we discuss the relaxation mechanism observed in our simulations as well as potential sources of error and how we expect that these will influence interpretation of the experimental observables.

Within this challenge, the objective has been to translate excited state dynamics into experimental observables. The importance of this cannot be understated. In many cases, collaboration between experimental and theoretical studies focuses upon the comparison of quantities that are easy to calculate, such as electronic state population kinetics. These kinetics are then compared to experimentally extracted timescales and agreement is taken as accuracy of the simulations. However, as shown in this work the LR-TDDFT (356 fs) and XMS-CASPT2⁶⁷ (335 fs) dynamics provide very similar decay kinetics but different predictions of photoproducts influencing the experimental signal. While slightly slower, the excited state dynamics performed using ADC(2) potentials also occurs on a comparable timescale (∼700 fs), but owing to the artificial crossing along the C=O bond stretch,⁶² the excited state decay occurs via a completely different mechanism.

This highlights one of the key choices that has to be made for any excited state dynamics simulation, i.e., the electronic structure method to simulate the excited state potentials. As discussed above, the single reference methods used in both this LR-TDDFT and ADC(2) exhibit challenges in describing the dynamics between electronically excited states and the ground state. While multireference methods such as CASPT2 will overcome these limitations, this gives rise to a significant computational burden, which, while possible for small molecules such as cyclobutanone, would be prohibitive for larger systems. In addition, their performance strongly depends on the choice of the active space, which is likely to be dynamic, i.e., it changes during the excited state dynamics, necessitating the use of a larger active space. Active space methods without perturbation theory, such as CASSCF, reduce the computational expense, but they do not treat dynamic correlation effects. While the influence of this will be system specific, for cyclobutanone this appears to speed up the excited state decay.^66,77

Our simulations indicate that after excitation of the 3s ← n Rydberg state, the system relaxes within 1–2 ps to form a broad range of photoproducts. Decay of this initially excited S₂ state occurs with a time constant of ∼350 fs. This is in good agreement with the fastest kinetics reported from previous time-resolved photoelectron experiments presented in the work of Kuhlman et al.¹¹ However, we note that Ref. 11 also reports a strong contribution from a slower time component, ∼750 fs. This is not observed within the LR-TDDFT population kinetics, but it is in very close agreement with the ADC(2) kinetics presented in the supplementary material. The exact origin for this difference between LR-TDDFT and ADC(2) is unclear; however, analysis of the hopping geometries indicates a flatter potential in the case of the latter, which permits a slightly wider spread of the trajectories in nuclear configuration space. Importantly, in both cases despite the small nuclear displacement required to reach the crossing point, the internal conversion from S₂ to S₁ is comparatively slow due to the symmetry-forbidden nature of the transition.

Once populated, the S₁ undergoes a large structural distortion, primarily along the C_α bond, consistent with previous work.^10,18 This drives the population to be rapidly transferred into the electronic ground state to form very vibrationally hot species. The fast nature of the population transfer from the S₁ to the ground state means that population of the S₁ does not exceed ∼30%. The excited state molecular dynamics consider only the dynamics within the singlet manifold. To assess the potential influence of the intersystem crossing and the triplet states, our quantum dynamics include the low-lying excited triplet states. These simulations indicate negligible amount of intersystem crossing into the triplet manifold, leading us to conclude that this channel will be unable to compete with internal conversion rates found here.

Although the excited state dynamics and high energy of excitation lead to some highly distorted geometries, our simulations point to the formation of the photoproducts being determined as the trajectory passes through the CI between the first electronically excited state and the ground state via a ring-opened intermediate. This is consistent with the conclusions from Ref. 10, whose authors demonstrated that motion away from the CI branching space leads to all of the observed photoproducts. Herein lies the most significant approximation within our work, as neither of the single reference method used will capture the biradical nature of the photoproducts. Indeed, although all of the CIs identified¹⁸ in previous work have been found within the LR-TDDFT(PBE0) framework and exist at accessible energies, i.e., below the excitation energy, our simulations show a much lower fraction of photoproduct formation than previous theoretical^66,67,73,74 and experimental¹⁵ works. This appears to most strongly affect the C3 photoproducts, which only form in 2.5% of the trajectories. In contrast, trajectories in the T₁ state, performed using unrestricted Kohn–Sham facilitating the description of biradical character, initiated at random from trajectories populating the S₁ state, show a much higher formation of C3 photoproducts (C3: 53%, C2: 5%, and ring-open: 17%) consistent with previous experiments.¹⁵ The motion through the CI and therefore the potential shape in this region is likely to be critical in determining the branching ratio of the photoproducts. Here, it may not only be a limitation of the single reference methods used but also a condition of the excited state dynamics. As stated in Sec. II C, our present dynamics attempts to avoid instabilities in the multi-configurational region, near the degeneracy by enforcing a hop to the ground state when the S₁–S₀ energy gap became $ΔES1−S0<0.1eV$⁠. While this avoids the explicit motion through the CI, the enforced earlier transition may also promote populated transfer closer the cyclobutanone structure, encouraging reformation of vibrationally hot cyclobutanone, rather than the photoproducts.

From these excited state molecular dynamics simulations, the ultrafast electron diffraction observable shows distinct changes and by studying the time-resolved PDF, these are largely associated with a loss in intensity for interactions at 1.5 and 2.5 Å, arising from dissociation. Despite the rich dynamics and the distinct changes observed, the time-resolved scattering curves show very little distinct dynamics largely associated with the incoherent nature of the dynamics and the comparatively low temporal resolution (150 fs).

Finally, a logical question would be ask if the limitations discussed above can be overcome within the present framework, i.e., without adopting a multireference wavefunction method, which would prove challenging for larger systems. Here, it is important to stress that while the relative yields of photoproducts formed appears somewhat at odds with previous works, all major reported products are generated, i.e., the full nuclear configuration space has been sampled. Excited state simulations have previously been used to simulate the experimental observables associated with structurally sensitive techniques of electron⁸⁴ and x-ray diffraction.³⁵ Importantly, in both of these works the outcomes of the trajectory-based dynamics were used as a basis to fit to experimental data and deliver an interpretation. Recently, forward-optimization for mapping trajectory basis functions onto time-resolved data has been shown to be highly effective.⁸⁵ While the use of a fit means that such an approach may not be classed as fully predictive, in both cases an excellent agreement between experiment and theory was achieved, providing deep insight into the dynamics observed.

The supplementary material contains additional electronic structure calculations, expansion coefficients for the Hamiltonian, and computational details of for the quantum dynamics. Excited state dynamics performed using ADC(2) potentials and key optimized geometries are also included.

This research made use of the Rocket High Performance Computing service at Newcastle University and the Viking High Performance Computing service at the University of York. We also acknowledge the COSMOS Program Grant (Grant No. EP/X026973/1). T.J.P would like to thank the EPSRC for an Open Fellowship (Grant No. EP/W008009/1).

The authors have no conflicts to disclose.

J. Eng: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Methodology (equal); Visualization (equal); Writing – review & editing (equal). C. D. Rankine: Data curation (equal); Formal analysis (equal); Visualization (equal); Writing – review & editing (equal). T. J. Penfold: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (lead); Methodology (equal); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal).

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