Benzophenone Ultrafast Triplet Population: Revisiting the Kinetic Model by Surface-Hopping Dynamics

The photochemistry of benzophenone, a paradigmatic organic molecule for photosensitization, was investigated by means of surface-hopping ab initio molecular dynamics. Different mechanisms were found to be relevant within the first 600 fs after excitation; the long-debated direct (S1 → T1) and indirect (S1 → T2 → T1) mechanisms for population of the low-lying triplet state are both possible, with the latter being prevalent. Moreover, we established the existence of a kinetic equilibrium between the two triplet states, never observed before. This fact implies that a significant fraction of the overall population resides in T2, eventually allowing one to revisit the usual spectroscopic assignment proposed by transient absorption spectroscopy. This finding is of particular interest for photocatalysis as well as for DNA damages studies because both T1 and T2 channels are, in principle, available for benzophenone-mediated photoinduced energy transfer toward DNA.


Computational Details
In this study, the ab-initio multiconfigurational quantum chemistry method CASSCF (Complete Active Space Self Consistent Field) S1 as implemented in MOLCAS 8 S2 was employed, averaging over the lowest-lying 3 singlet or triplet states (SA-3-CASSCF). The calculation of several energy gradients is necessary for each step of each trajectory. The analytical gradients of each state in an interval of 0.25 eV from the active state were calculated in parallel and the Cholesky decomposition was used to approximate the two-electron integrals, in order to significantly reduce the computational effort (by ca. a factor 3 on our local linux cluster based at the Université de Lorraine). Spin-orbit couplings were calculated based on the SA-CASSCF wavefunctions using the RASSI module of MOLCAS, which calculates spin-orbit matrix elements according to the AMFI (atomic mean field integrals) formalism. S3 The values of the spin-orbit coupling elements were found to be of the order of 20 cm -1 , in agreement with previous computational studies (see ref. 21 in the main text). 65 initial conditions (Cartesian coordinates and atomic velocities) were generated by sampling from a Wigner distribution. S4,S5 As input, a set of vibrational frequencies and the corresponding normal mode vectors were provided by a frequency calculation performed on the Franck−Condon geometry at the CASSCF(12,11)/ANO-S-VDZP level of theory. The calculated kinetic and potential contributions to the total energy are of the same order: the kinetic energy is 0.087 ± 0.019 a.u. and the potential energy is 0.099 ± 0.019 a.u. A 0.5 fs time step was selected for each trajectory, while the electronic wavefunctions were propagated with a 0.02 fs step. The dynamic was run in the full diagonal representation including non-adiabatic and spin-orbit coupling using the SHARC code. S6 To calculate the non-adiabatic interaction between the different states we considered the wavefunction overlaps ⟨ ( )| ( )⟩, within the local diabatization formalism. S7 The spin-orbit couplings were explicitly calculated at each time step, considering all multiplet components. During the dynamics, the velocity is rescaled to adjust the kinetic energy after a surface hop and energy-based electronic decoherence correction with a parameter of 0.1 Hartree S8 was applied. On the one hand, we note that more sophisticated electronic decoherence corrections are available. S9-S12 Nevertheless, they are adding significant extra cost to the simulations or require considerable implementation effort. On the other hand, the Granucci−Persico algorithm S8 is a conceptually and computationally simple means of estimating electronic decoherence, and we therefore used it. All details concerning the SHARC code can be found at https://sharc-md.org.

Kinetic Models Parameters
Two kinetic models were applied, called parallel and serial kinetic models. The parallel model takes into account the simultaneous formation of T 1 and T 2 from S 1 , as shown by the surface hopping dynamics study, and it is formulated as follows: where t is time, k is the rate constant value and τ is the lifetime value. The serial kinetic model considers the formation of either T 1 or T 2 (or T 3 ) from S 1 , as usually postulated in the fitting procedures of the transient absorption experimental data, based on signal strengths. It can be formulated as follows: For the serial kinetic model, a full decomposition of the triplet contribution is given in Figure S1. Moreover, all lifetime values and relative fitting errors are given in Table S1, including the contribution from the singlet low-lying excited state.
As can be seen, the sum of T 1 and T 2 contributions accounts for the whole intersystem crossing process, while T 3 can be neglected. Figure S1. Global fitting to a single exponential rate law for T 1 , T 2 , T 3 , and their sum. For the parallel kinetic model, the rate constant values and relative fitting errors are given in Table S2. Table S2. Parallel kinetic model: rate constant values and relative fitting errors for the S 1 →T 1 (k 1 ) and S 1 →T 2 (k 2 ) pathways.
Pathway k (fs -1 ) Δk (fs -1 ) Δk (%) S 1 → T 1 0.00111531 1.129•10 -5 1.01 S 1 → T 2 0.00029563 1.002•10 -5 3.39 It should be noted that the relative fitting errors shown in Tables S1 and S2 take into account the collected data from the sampled trajectories, but do not include the systematic errors due to the chosen level of theory and to the applied surface hopping algorithm. A comparison of the two proposed kinetic models up to 5 ps is shown in Figure  S2a,b. As it can be seen (also in Figure 1 of the main text) the two curves almost coincide up to 600 fs, while the population plateau for the triplet states is reached at larger times for the serial kinetic model. Nevertheless, while the overall population of the theoretically oriented parallel model adds up to one, this is not the case for the experimentally oriented serial model. Indeed, the serial model is intended to reproduce the experimental fitting of signal strength that (unlike the overall population) does not have to add up to one. In order to get a deeper insight into the singlet-triplet and triplet-triplet transitions, Figure S2c shows net transfers between couples of electronic excited states, based on the number of hopping events recorded along the trajectories. We found, confirming our models, net transfers from S 1 to both T 1 and T 2 , with a slight preference for the S 1 →T 2 channel, i.e. the indirect mechanism. Once the triplet manifold is populated, an equilibrium is established between T 1 and T 2 states, with a net transfer (14 hops) from T 2 and T 1 , indicating an emerging equilibrium (the number of T 1 →T 2 and T 2 →T 1 hops is comparable but not identical, since the equilibrium first has to be established).

Figure S2
. Serial (a) and parallel (b) kinetic models including the dynamics data for the first 600 fs and the analytical curves up to 5 ps. Scheme of the net transfers between S 1 , T 1 and T 2 electronic states, based on the number of hopping events: 86 S 1 →T 1 hops, 70 T 1 →S 1 hops, 349 S 1 → T 2 hops, 327 T 2 → S 1 hops, 50 T 2 →T 1 hops, 36 T 1 → T 2 hops (c). A net transfer is the difference between forward and backward hopping events.

Single Point Calculations of Geometries Along the Potential Energy Surface
The CASSCF(12,11) benchmark on Franck−Condon and S 1 minimum geometries was performed for different basis sets (Table S3)  The values reported in Table S4 refer to benzophenone in water solvent, taken into account by a polarizable continuum model (PCM) through state averaged calculations at the CASSCF level. 110.64 Figure S3 shows the agreement between CASSCF(12,11)/ANO-S-VDZP and SS-CASPT2(12,11)/ANO-L-VDZP along the excited state MEP. As can be seen S 1 , T 2 and T 1 energies are close (interval of 0.35 eV) at the T 1 minimum geometry. We note that the CASPT2 level of theory is more accurate than the CASSCF level of theory. Nevertheless, presently CASPT2 dynamics are not feasible, since the present implementations of the CASPT2 analytical energy gradient are impractical for molecules of the size of BP. The same can be said about the use of the CASPT2 numerical energy gradient.

Wavefunction Amplitude Graphs of the 39 Trajectories Populating the Triplet manifold
The following graphs show the wavefunction amplitude as a function of the time (fs), for all 39 trajectories where a successful population of the triplet manifold was observed. All electronic states considered for the simulation (the singlets S 0 , S 1 , S 2 and the triplets T 1 , T 2 , T 3 ) are shown. The wavefunction amplitude is calculated as the sum of the absolute squares of the MCH (molecular Coulomb Hamiltonian) coefficients for each state. For triplets, multiplet components are summed up. Most of the trajectories were stopped when the sum of the wavefunction amplitudes of a certain triplet state -including multiple components -equal to one; while 15 trajectories were calculated for additional time, in order to study the kinetic equilibrium between triplet states, and check the eventuality of re-crossing to populate back the singlet state (observed only once).

Geometrical analysis of the 39 Trajectories Populating the Triplet Manifold
A geometrical analysis was performed, showing the evolution of d CO , Φ plan and θ pyr (see Scheme 1 in the main text) with time. Figure S4. Distance d CO as a function of time. Figure S5. Improper dihedral Φ plan as a function of time. Figure S6. Improper dihedral θ pyr as a function of time.

Franck−Condon Structure at the CASSCF(12,11)/ANO-S-VDZP Level of Theory
Cartesian coordinates of the ground state minimum (in Ångström):

Description of the Transient Absorption Spectrum Simulation
The transient absorption spectrum (see Figure 2c in the main text) was simulated directly from non-adiabatic dynamics. For each trajectory at each time step we considered vertical energy differences from the current diagonal state to all the other states. Oscillatory strengths between current state and all the other electronic states, i.e. absorption or emission intensity, were directly calculated by the MOLCAS software. Oscillatory strengths were set negative for transition to lower energy states as compared to the current state, and positive elsewhere. This in turn simulates stimulated emission and excited state absorption, respectively. The set of vertical transitions, in eV, was convoluted in the energy domain placing a Gaussian function of full-width at half-length (FWHL) of 0.1 eV, and subsequently converted to nm. Regrouping all the results for all the time steps in a single matrix gave the final map, plotted with the gnuplot program. We note that this procedure corresponds to the calculation of only the inhomogeneous part of the transient absorption spectrum, since it is based on only the energy gaps between electronic states. In order to include the homogeneous component of the transient absorption spectrum, the calculation of the dynamical Franck−Condon factors would be required. S13,S14 Concerning the interpretation of the transient absorption spectroscopy experiments, we note that for both − FC and S 1 minimum − geometries T 3 is higher in energy at the CASSCF level than at the CASPT2 level. Nevertheless, as also mentioned by Sergentu et al. (ref. 21 of the main text), at the CASPT2 level in the FC region the SOC value between S 1 and T 3 was found to be not significant, while SOC values of 24 and 22 cm -1 were found for S 1 -T 2 and S 1 -T 1 states, respectively (ca. 20 cm -1 in our study). Therefore, the fact that T 3 is high in energy at FC should not affect the reactivity observed in our study, where T 3 never plays a significant role. The same holds for the S 1 minimum energy region, since in this case our trajectories are populating the singlet manifold, and therefore they are not used for the interpretation of the transient absorption spectroscopy experiments, that do imply population of the triplet manifold.
As can be seen in Figures 5 and 7 of ref. 21 main text, the CASPT2 MEP predicts that the T 2 -T 3 almost degeneracy is left when evolving from FC to the S 1 -T 1 -T 2 crossing region, corresponding to our T 1 -T 2 dynamic equilibrium. In this case, T 3 rises in energy, at both CASSCF and CASPT2 levels of theory. Our considerations about the interpretation of transient absorption spectroscopy experiments do not refer to the FC region, but after triplet population, when finally T 1 and T 2 are almost degenerate and T 3 is indeed upper in energy. Therefore, even considering CASSCF deficiencies, we expect an almost identical T 1 -T 3 and T 2 -T 3 energy difference (once the triplet manifold is populated and T 1 -T 2 equilibrium established). Hence, since T 1 →T 3 corresponds to a (n,π*) transition and T 2 →T 3 corresponds to a (π,π*) transition, the T 2 →T 3 transition results in an optically much brighter (higher oscillator strength) transition than the T 1 →T 3 transition.