Coherent Mixing of Singlet and Triplet States in Acrolein and Ketene: A Computational Strategy for Simulating the Electron–Nuclear Dynamics of Intersystem Crossing

We present a theoretical study of intersystem crossing (ISC) in acrolein and ketene with the Ehrenfest method that can describe a superposition of singlet and triplet states. Our simulations illustrate a new mechanistic effect of ISC, namely, that a superposition of singlets and triplets yields nonadiabatic dynamics characteristic of that superposition rather than the constituent state potential energy surfaces. This effect is particularly significant in ketene, where mixing of singlet and triplet states along the approach to a singlet/singlet conical intersection occurs, with the spin–orbit coupling (SOC) remaining small throughout. In both cases, the effects require many recrossings of the singlet/triplet state crossing seam, consistent with the textbook treatment of ISC.

O ur focus in this article is the description of coupled electron−nuclear dynamics in molecules with spin−orbit coupling induced superpositions of electronic states of different spin multiplicities. This is analogous to attochemistry, 1−3 where the superpositions are created in laser experiments.
In photochemical dynamics, radiationless decay�internal conversion (IC) between states of same spin multiplicity and intersystem crossing 4−6 (ISC) between states of different spin multiplicity�is conventionally described in terms of descent through electronic states linked by potential energy surface crossings/conical intersections. 7−11 . 12,13 Many nonadiabatic dynamics methods, e.g., surface hop, describe this effect by moving on one electronic state surface at a time. 6,14 Here we explore an alternative formalism 15 where radiationless transitions are continuous and take place on a single time-dependent potential energy surface as opposed to being discontinuous (e.g. via a surface hop) at a crossing. Thus, a singlet + triplet state superposition will generate a gradient for nuclei that is different from either the singlet or triplet alone (and that will evolve with time), and the subsequent nuclear motion will be determined by the mixed state. At the start of the dynamics, in addition to the intrastate gradient components, there is an interstate gradient component that arises from the mixing. At later times, there is an instantaneous gradient that arises from the electron dynamics (as demonstrated before in the glycine 16 and benzene 17

cations).
Intersystem crossing in photochemistry is known to be relatively slow. 18 In this work, we will see this effect as a slow buildup of a mixed singlet + triplet state after repeated passage across the crossing region. The reaction path for ISC will therefore be a time-dependent mixture of singlet and triplet states determined by the solution of the time-dependent Schrodinger equation (TDSE).
We now discuss the formalism employed; namely, the 1 electron average spin−orbit coupling Hamiltonian and how this is included in an Ehrenfest dynamics framework.
Spin−orbit coupling (SOC) mixes electronic states with different spin multiplicities. In the context of photochemical dynamics, SOC leads to intersystem crossing (ISC), and the dynamics of the ISC is expected to be governed by the size of the SOC between the states involved. SOC is a relativistic effect 19 whose contribution can be approximated by the socalled Breit−Pauli term, 20 and contribution to the molecular Hamiltonian in atomic units can be expressed as where the nuclei are indexed with capital letters and the electrons with lowercase letters, Z A is the charge on nucleus A, r iA is the distance between nucleus A and electron i, r ij is the interelectron distance between electrons i and j, l̂is the orbital angular momentum operator, Ŝis the spin operator, and c is the speed of light.
To avoid the computationally expensive two-electron term we use the approximation of Koseki et al. 21,22 wherein the twoelectron term is incorporated into the one-electron term using an empirical effective nuclear charge (Z A,ef f ) defined in ref 21: SOC also acts to break the degeneracy of multiplet microstates (doublets, triplets, etc.) (for example, the well documented splitting of the sodium 2 P state into 2 P 1/2 and 2 P 3/2 by approximately 2 meV. 6 This effect is relatively small for light atoms, and, as is common practice,we treat the multiplet microstates (e.g., m s = −1, 0, +1 in a triplet) as one "average" state for the purposes of dynamics. 23−26 In this work we employ this averaged triplet approximation and test it by evaluating the triplet splitting at points of high singlet/triplet mixing along the trajectories using a relativistic calculation 27 in the SI. In our dynamics simulation, we use the Ehrenfest method 15,28−31 which describes the electronic wave function in a Complete Active Space Configuration Interaction (CAS-CI) formalism. The Ehrenfest electronic state is expressed as a linear combination on the basis of all Slater determinants formed in the active space with zero z spin component (so that both singlet and triplet states can be represented). A related approach is exact factorization, 32,33 which can also be adapted to include ISC as has been described by Talotta et al. 34 In this article we will examine the intersystem crossing dynamics for acrolein and ketene (shown in Scheme 1) which exhibit different mechanistic effects of ISC. The potential energy surface topology of these simple organic carbonyl− allylic systems has been well characterized with a wide range of other electronic structure and dynamics methods. 26,35−45 In acrolein, the spatial/orbital configurations of singlet and triplet are different, whereas in ketene they are the same (Scheme 1). As a result, we have a large difference in spin−orbit couplings: ∼65 cm −1 in acrolein versus ∼0.1 cm −1 in ketene at the S 1 minima. According to El-Sayed's rule 46 we should therefore expect the yield of triplet in acrolein to be greater.
Our focus in this article is the mechanism of ISC on a mixed state. We will take an exploratory approach that is similar to the reaction path idea that is central to chemistry. Our computations are intended to explore the singlet/triplet crossing region of the potential energy surface in a way similar to that in which one may explore the transition region of a thermal reaction. We now elaborate on this idea briefly.
We start from the fact that passing through a singlet/triplet crossing is a relatively rare event (like passing through a transition state 47−51 ). In this way we can draw an analogy with the reaction path idea 52 that is central to mechanistic chemistry. 53 The reaction path approach is most commonly used with transition states. 53,54 It can be realized in dynamics (a dynamical reaction path) by starting a trajectory with the initial momentum in the direction of negative curvature. 48,55 Alternatively, one could start a trajectory at the "reactants", with an initial momentum in the direction of the transition state. The results presented in this work are an excited state analogue of this technique; we start from an excited singlet minimum (our "reactants") (labeled S 1 in both systems) with activation of the vibrational normal mode that corresponds to motion toward the S/T crossing seam (our "transition state")�i.e. the normal mode approximately parallel to the gradient difference vector at the minimum energy point on the crossing seam ( Figure 1). The essential difference, as we shall show subsequently, is that when the trajectory reaches the crossing seam, it does not proceed to "products" but rather results in a mixture of singlet and triplet states, and after reaching a "turning point" proceeds backward across the seam again. This occurs many times with the population of triplet increasing on each passage of the crossing seam. Thus, we believe that our single "reaction path" trajectory yields some mechanistic insight of ISC. Scheme 1. Electronic States and Their Dominant Configurations from a SA-CASSCF(6,5)/6-31G* Calculation at the S 1 Minimum Geometry for Acrolein and Ketene a a We illustrate the M s = 1 configuration of the triplets for visual clarity (in the dynamics we use the M s = 0). Our S 1 and T 1 states for ketene are the same as those Xiao e. al., 37 with their work characterizing them as π ⊥ → π || * states, where the π ⊥ orbital is the second and π || * is the fifth. In this work we are not attempting to calculate the rate of triplet formation; in our formalism ISC presents itself as a slow buildup of the triplet component of the mixed state.
We now summarize the main chemical concepts associated with ISC as seen in acrolein and ketene. Scheme 2a shows that the acrolein dynamics proceeds with many S 1 /T 2 recrossings resulting in a continuous increase in the triplet component of the mixed state that develops. For ketene, we have an important effect of the mixing of singlet and triplet on the nuclear motion of the singlet state. We sketch this effect in Scheme 2b,c. Ketene has both a singlet/singlet conical intersection and a triplet/singlet crossing near the starting minima. As we shall show, the mixed S 1 /T 1 state created by spin−orbit coupling interacts in a different way with S 0 and thus changes the dynamics of the singlet−singlet nonradiative decay (internal conversion) at the S 0 /S 1 conical intersection (Scheme 2c). The numerical details will be provided subsequently.
We now discuss how we include SOC within the Ehrenfest approximation of the TDSE for electron motion as devised by Vacher et al. 15 In eq 3 t n is the current timestamp, t n−1 is the previous timestamp, and H t ( ) n is the matrix representation of the CI Hamiltonian in a basis of Slater determinants within a CAS-CI 56 expansion. We use Slater determinants because we will mix different S z = 0 states when the SOC is included. The Ehrenfest vector A t ( ) n contains the weights of the Slater determinants in the time-dependent wave function.
The SOC magnitude matrix elements involve only the oneelectron operator (eq 2) and are collected in the matrix t ( ) n in the same basis as H t ( ) n (in fact, it turns out to be convenient to evaluate these matrix elements in the CI eigenvector basis and transform back�the full implementation details are provided in SI), and we recompute this SOC matrix at every time step.
Once this electronic propagation step is done, the updated Ehrenfest vector A can then be used to compute the molecular properties at each time step�namely the energy, gradient, and Hessian of the electronic wave function. 15 These quantities are then used to propagate the nuclear wave function. We wish to emphasize that the time-dependent Ehrenfest potential is not a Because of the spin−orbit coupling term, A t ( ) n will evolve to become a mixture of singlet and triplet states. Since Slater determinants are not, in general, eigenstates of the spin operator Ŝ2, to interpret the results, it is necessary to consider a projection of A t ( ) n onto possible states that are either singlet or triplets. The simplest choice is projection of the Ehrenfest vector onto a set of spin pure Hartree−Waller 57 functions or Clifford Algebra spinors. 58 The Hartree−Waller (HW) functions 59 are defined as a singlet Or a triplet combination of Slater determinants where |φ 1 A ⟩ and |φ 2 A ⟩ are two determinants that have the same spatial orbitals A but with all the spins inverted. Any singlet adiabatic state must be a linear combination of singlet HW functions, and any triplet must be a combination of triplet HW functions.
As previously discussed, upon inclusion of SOC effects, the vector A t ( ) n may have contributions from both triplet and singlet HW functions. If we sum the populations of all pairs of eq 5 and eq 6, we obtain the total triplet and singlet populations of the mixed state.
Finally, if one inspects the individual coefficients of |φ 1 A ⟩ and |φ 2 A ⟩ in A t ( ) n one can distinguish two cases. The coefficients will be equal if the coupling for spatial arrangement A is either a singlet or triplet (this is the situation in acrolein, as we shall show because we mix a n−π* singlet with a π−π* triplet). However, if individual coefficients of |φ 1 A ⟩ and |φ 2 A ⟩ in A t ( ) n are not equal then both singlet and triplet states are present with the same spatial arrangement. In this case, the singlet triplet mixed state has the same orbital components A (this is the case in ketene where singlet and triplet states are both π−π*).
In this article both systems (acrolein and ketene) have been studied with the second order Ehrenfest 15 dynamics method.
The electronic wave function is propagated in a CAS-CI formalism starting from a converged CASSCF calculation (described above) as implemented in a development version of Gaussian DV J05. 60 This code also computes the gradient and Hessian of the Ehrenfest state. The nuclei are propagated on the Ehrenfest surface using a fifth order predictor corrector scheme of Schlegel 61 et al. (also, as implemented in Gaussian). For both systems we have used an integrator step size of 10 −2 bohr amu · (leading to an average time step of ∼0.06 fs). We started the nuclear dynamics (a single representative trajectory) with an initial energy of 20 kcal/mol in acrolein and 24 kcal/mol in ketene in their respective singlet−triplet gradient difference vector (GDV) normal modes. The initial energy in those modes was chosen after some experimentation so that the crossing region was adequately explored because the "reaction path" is somewhat curved.
We shall present the data for 2 simulations for each system with the same initial conditions but for the inclusion of spin− orbit coupling effects. The differences in the resultant nuclear and electronic trajectories illustrate the mechanistic effects of the SOC.
We will now provide the numerical details that support the qualitative picture outlined in Scheme 2.
In both acrolein and ketene, the dynamics is started on the S 1 state from the S 1 minimum with an initial momentum of 20 kcal/mol in S 1 /T 2 GDV and 24 kcal/mol in S 1 /T 1 GDV, respectively (Scheme 2). From there the trajectory rapidly approaches the singlet/triplet crossing. What happens next The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter depends on the shape of the crossing; in acrolein the S 1 /T 2 crossing is slightly peaked, and we see that one requires many recrossings of the seam (discussed in Figure 2aii) to reach a significant population of the triplet. In contrast, in ketene, the S 1 /T 1 crossing is sloped, and the surfaces are almost parallel since the orbital components of the two states are the same (Scheme 1) and the triplet population remains very small (to be discussed in Figure 2bi). This is consistent with the textbook 7,8 description of ISC as a slow process requiring many recrossings of the singlet−triplet seam (Figure 2aii+bii), from which one would expect a slow transfer of population ( Figure 2ai+bi) onto the triplet at every passage of the crossing seam. For acrolein, as the molecule passes through the S 1 /T 2 seam (Figure 2aii,aiii), singlet population is transferred onto the triplet (Figure 2ai). As a result, propagation takes place on an effective surface governed by the evolving superposition of the singlet and triplet states (Scheme 2). In the simulation without SOC effects, there is no triplet population as the molecule simply rolls around the S 1 minima well. The associated effect of S 1 /T 2 mixing on geometry changes is discussed in the SI ( Figure S2). There we show a C−C−C bending normal mode (SI Figure S2) becomes enhanced at around 200 fs, consistent  Figure 1) versus NM4 or NM1 (components of S 1 /S 0 GDV). Both trajectories start at the S 1 minima at (0,0). The blue curve is with the spin−orbit interaction disabled, and the orange curve is with the SOC included. The trajectories are truncated at 150 fs for clarity. c: The population of 2 dominant Slater determinants with π−π* orbital configuration (|φ 1 A ⟩ and |φ 2 A ⟩ in eqs 5 and 6). Note the splitting after ∼110 fs demonstrating that we are on a mixed singlet/triplet state (due to the phase difference in HW functions). d: Population of the dominant HW functions, S 0 HW, etc. (corresponding to Scheme 1). e, f: Same as c and d, but for trajectory with SOC disabled.
with the sharp rise in triplet population in Figure 2ai and the difference in primary orbital occupancy (S 1 vs T 2 ). Internal conversion from T 2 to T 1 is also possible and could be captured by our dynamics but is not observed in our results on the time scale we have simulated.
In ketene, the nuclear dynamics leads to a singlet−singlet conical intersection with a nearby singlet/triplet seam (a 2 singlet +1 triplet system has previously been documented by Granucci et al. 62 ). To reach the S 0 /S 1 seam, one must excite the associated S 0 /S 1 GDV normal mode, a mixture of NM1 and NM4 (Figure 1). The amplitude of these modes is enhanced with SOC included, similar to acrolein, and thus, the dynamics of the passage across the S 0 /S 1 conical intersection is changed by spin−orbit coupling. To quantify this effect, we show the paths traced out by the trajectory in normal mode space for the first 150 fs (Figure 3a,b). These plots allow us to explore the correlation of 2 important coupling normal modes of the trajectory, as was done previously 63 by Lasorne et al. Since there is a transfer of kinetic energy between the x-axis normal mode (NM7 primary component of T 1 /S 1 GDV) and the y component (NM1 and NM4 components of S 0 /S 1 GDV), we observe a Lissajous figure-like trace. These plots show how we reach the singlet/singlet seam�via a slow transfer of energy to the necessary singlet/singlet GDV normal mode.
As a reference point, we shall address the dynamics without a SOC first. As shown in Figure 3a,b,e,f the trajectory slowly takes the system to the S 1 /S 0 seam after ∼370 fs (Figure 3f).
In Figure 3a−d we show the effect of including SOC for ketene; despite a small triplet population (Figure 2bi) it acts to enhance the transfer of kinetic energy from NM7 to NM1 (Figure 3a,b) and NM4 leading to a faster approach to the S 0 / S 1 seam (within ∼120 fs, Figure 3d). This transfer of kinetic energy between normal modes is reflected in Figure 3a,b (orange): the trajectory curves with inclusion of SOC gain more y-axis amplitude than the ones with SOC effects disabled over the same time frame. The small admixture of T 1 has thus enhanced the S 1 /S 0 internal conversion by placing more momentum in the S 0 /S 1 GDV modes. Notice the rapid oscillation of the diabatic dominant singlet HW populations (Scheme 1) in Figure 3d,f after passing through the conical intersection. This is characteristic of electron dynamics (charge migration). 64 In fact, the oscillations in the electronic state give rise to the oscillations in the nuclear motion, and the motion of the nuclei induces the electronic dynamics.  (Figure 3c) the 2 populations diverge due to the phase difference between singlet and triplet HW functions, which clearly illustrates that our dynamics is running on a mixed state.
We have described the mechanism of intersystem crossing as a superposition of singlets and triplets created under the influence of spin−orbit coupling. In this description, the composition of the mixed state controls the nuclear dynamics, rather than the textbook 7,8 description where one propagates on either a pure singlet or pure triplet surface after the ISC event. Nevertheless, the qualitative ideas about state populations involved with intersystem crossing do not change: the coupling of the singlet and the triplet is controlled by passing through the seam of the singlet/triplet conical section and many recrossings are necessary for even partial population transfer. The distinguishing mechanistic feature of this work is that the nuclear trajectory is driven by the unique gradient of the coherent superposition of singlet and triplet, which is distinct from the gradient of either singlet or triplet.
In this work, we have reported results for two chemical examples: acrolein with a relatively large spin−orbit coupling and ketene where the spin−orbit coupling is much smaller. In acrolein we saw a large population of triplet, and an effect on geometry (see SI). However, in ketene, the population of triplet was small, but the gradient of the mixed state strongly influenced the singlet/singlet conical intersection approach. El Sayed's rule suggests that ISC is much faster in acrolein than in ketene, and indeed this agreement is a validation of our approach.
On a point of method, one might expect the surface hop method to yield reasonable results for acrolein but much less so for ketene. We base this prediction on the fact that SOC (and triplet population) is so small in ketene that a majority of trajectories would simply never make the hop onto the triplet�hence never experience the extra excitation of the S 0 / S 1 GDV normal mode that leads to the faster arrival at the singlet/singlet seam. Of course, one could regard the dynamics presented here as a first step in a hierarchy that might involve many trajectories with sampling, either classical, quantum mechanical, 65,66