Surface hopping investigation of the relaxation dynamics in radical cations

Radical cations are created through ionization processes when high energy photons are present. Ionization is often used as a probe in various pump-probe techniques applied to examine dynamics in molecules. Time-resolved photoelectron spectroscopy^1,2 and strong-field dissociative ionization are two such techniques which can be used to probe excited state molecular dynamics. Strong-field dissociative ionization leads to the formation of fragments that can be detected with the help of time-of-flight mass spectrometry (TOFMS). The fragments that are formed depend on several factors, including the initial neutral electronic state as well as the final ionic state and the geometrical distortion of the molecule at the time of fragmentation. Thus, those fragments give a hint to the underlying dynamics.^3,4

In order to understand and interpret strong-field ionization and fragmentation experiments, we need to understand the various steps involved in the process, that is, the ionization and the fragmentation steps. In previous studies, we focused on the ionization process induced by strong fields and found that the intuitive thinking of ionization based on Koopmans’ theorem does not always apply in this case.^5,6 More recently, coincidence measurements of photoelectrons and the corresponding cations using velocity map imaging were used to explore whether ionization occurs directly or indirectly.^7,8 It was found that the energy gap between dissociative and non-dissociative states is an important factor in the ionization to excited states.

While we have made progress in our understanding of ionization in strong fields, the fragmentation step and the dynamics of the cations leading to it are less understood. Fragmentation can take place directly by ionization to a dissociative cationic state or indirectly after deactivation from the initial cationic state to the ground state using the available energy for dissociation. Indirect fragmentation assumes fast relaxation to the ground state before dissociation and seems reasonable based on the high density of ionic states in the radical cations. It is similar to Kasha’s rule⁹ for neutral excited states where fast relaxation from higher excited states to S₁ is assumed followed by emission from the S₁ state, although the time scale for spontaneous emission in neutral molecules is much longer than the time scale for dissociation in cations. We have used the hypothesis of fast radiationless decay to D₀ to examine fragmentation of the uracil radical cation in the past, assuming that barriers on the ground state D₀ are responsible for the dissociation pathways.¹⁰ Nevertheless, a more detailed investigation on how fast relaxation occurs in radical cations is necessary in order to check the validity of this hypothesis. This is the subject of the current study.

The main questions we want to address in this work are: How does reactivity/dissociation on an excited state compete with relaxation to D₀? Is ultrafast radiationless decay to D₀ the norm for the ionic states? Are there any common features in the dynamics of cations? Can we predict the dynamics based on the properties of states upon ionization? Although the dynamics of radical cations have been studied in many theoretical studies in the past,^11–18 this work tries to focus specifically on the above questions. In order to address these questions, we investigate the dynamics of three cations, that have already been of interest in our work previously. These are the 1,3,5-hexatriene (HT) and 1,3-cyclohexadiene (CHD) cations which are connected by ring-closure/opening processes, and the cation of the RNA base uracil. The structures are shown in Figure 1.

The conversion of CHD to HT has been extensively studied as a prototype photoisomerization reaction.¹⁹ The excitation of CHD to its S₁ state eventually leads to a conical intersection (CI) that facilitates the isomerization and leads to a branching ratio for the two isomers. In previous work, we attempted to control the ratio of the two isomers using ultrashort laser pulses.³ If cis-HT is reached other isomers of HT can be formed as well by torsional motions around the single and double bonds. The fragmentation patterns seen in TOFMS for HT and CHD are very different, with many fragments produced in HT while mostly the parent ion is seen in CHD. This pattern was used as a detection scheme in our previous work, and we explained the differences in fragmentation by the spacing between ionic states in the two cations.³ 

The second system we are investigating here is the uracil radical cation. Uracil is an important biological molecule and its excited states have been studied extensively the last decade.^20–23 We recently used strong-field dissociative ionization to distinguish between two different mechanisms that had been proposed for the observed excited state lifetimes.²⁴ In order to interpret the results, we needed a good understanding of the ionization and fragmentation steps, and we investigated those in detail.¹⁰ An assumption that we made was that fragmentation of the cation occurs after relaxation to D₀. In addition, we have examined the potential energy surfaces of several states in the cation including CIs,²⁵ and we have simulated dynamics and spectra using the Multiconfigurational Time-Dependent Hartree (MCTDH) method.²⁶ This wealth of information on this system can serve as a guide and benchmark in our current study. Of particular interest will be the comparison between the MCTDH and Trajectory Surface Hopping (TSH) dynamics, as described below.

In this manuscript, we describe the dynamics on the doublet ionic states of the radical cations, CHD, HT, and uracil using TSH. Sec. II describes the electronic structure applied as well as the TSH details. A description of the electronic states and important features of their potential energy surfaces for the three cations follows. Finally, the results of the dynamics are analyzed and compared between the systems to see if global lessons can be learned.

The reference geometry for all TSH simulations is the corresponding optimized neutral ground state (S₀) minimum. Initial geometries and velocities were generated with uncorrelated and correlated Wigner distributions of 150–300 geometries using the Newton-X program.^27–29 This assumes a harmonic oscillator for each nuclear coordinate which is taken to be a normal mode. Correlated and uncorrelated Wigner distributions were tested and found to give similar results. We will present results from uncorrelated distributions for HT and uracil and from correlated distributions for CHD. The results for correlated and uncorrelated Wigner distributions in CHD are almost identical, and the uncorrelated results are shown in Figure S3 of the supplementary material.³⁰ All trajectories were run using the local diabatization method.³¹ This approach does not use derivative couplings and may be more stable than the dynamics when the derivative coupling is used. The approach was tested initially by running TSH using either local diabatization or derivative coupling and it was found that very similar results were obtained. The results of test runs with both methods are shown in Sec. 4 of the supplementary material.³⁰ Non-adiabatic effects were taken into account by the fewest-switch criterion by Hammes-Schiffer and Tully.^32,33 The momentum after a surface hopping was adjusted according to the momentum direction. Decoherence effects were accounted for by using the method of Persico and Granucci with a parameter of α = 0.1 hartree.³⁴ A time step of 0.5 fs was used and a total time of 300 fs was investigated. In the uracil cation, the dynamics were run only for 120 fs due to the high computational cost of these simulations.

The initial population is on the D₂ adiabatic state in all systems, except HT where we have also run a set of trajectories starting from D₁. For CHD and uracil, ionization simulations and experiments indicate that D₂ is populated in strong-field ionization experiments, and this justifies our initial conditions. In HT we have trajectories starting from two different states so that we examine the effect of the different initial population on the dynamics. The number of trajectories are as follows: 150 trajectories were launched for HT in both D₁ and D₂ simulations, 186 for CHD and 150 for uracil. Trajectories that failed before the end of the simulation time (300 fs for CHD, HT, and 120 fs for uracil) were not included in the population plots but they were analyzed to determine whether they contribute to reactivity. 32, 42, and 12 trajectories failed for HT, CHD, and uracil, respectively. In some cases, a trajectory that fails may be an indication of reactivity, since a broken bond may cause the quantum chemistry employed to fail. Plots of populations including failed trajectories are shown in Figure S4³⁰ in order to demonstrate that neglecting these trajectories does not affect the decay rates we will discuss.

The adiabatic state populations and the accumulated jumping probabilities were calculated. Accumulated jumping probabilities are obtained as the number of jumps D_i → D_j minus D_j → D_i divided by the number of trajectories. So, they give the net probability of a trajectory having performed a jump after a given amount of time.

Additionally, we have performed Normal Mode Analysis (NMA)^35–37 which is implemented in the Newton-X program. For each trajectory and time t, the Cartesian difference vector between the current and a reference geometry is taken which is then projected onto the normal modes. This can be done in selected time intervals of the trajectories. The new vector in the normal mode basis is averaged over all trajectories to obtain an averaged trajectory, and the standard deviation of the averaged trajectory is calculated. The individual normal modes contribute according to their importance to that standard deviation. Large individual standard deviations mean that the corresponding normal mode dominates the nuclear motions in the considered time interval. To remove translation and rotation before carrying out the NMA, the structures of each trajectory were aligned with the S₀ minimum with a least-squares fit method, also part of the Newton-X module.³⁸ 

Optimizations of minima and calculations of vibrational frequencies for CHD and HT were carried out with Density Functional Theory (DFT) using B3LYP/6-31G(d) and the Gaussian03 program.³⁹ Both isomers possess C₂ symmetry in their S₀ and D₀ minima.

The location of CIs and the TSH dynamics were carried out using the complete active space self-consistent field method⁴⁰ (CASSCF) and the cc-pVDZ basis set. Specifically, the state-averaged version (SA-CASSCF) as implemented in the Columbus package was used.^41,42 For HT, the active space includes six π orbitals and the corresponding five electrons (denoted (5,6)), while for CHD, four π orbitals and two σ orbitals are included in the active space. Because of the limitations of CASSCF in CHD (which will be explained in Sec. III A) restricted active space SCF⁴³ (RASSCF) calculations were performed as well. The electronic part of the dynamics simulations was done using CASSCF and RASSCF for HT and CHD, respectively, without symmetry restrictions.

CIs for CHD and HT were optimized at the CASSCF(5,6) level. Multi-Reference Configuration Interaction with single excitations out of the active space (MRCIS) was used for the location of the minimum of the three-state CI seam. The location of the three-state CI was conducted with the algorithm for optimizing three-state CIs⁴⁴ implemented in a modified version of Columbus.^41,42

The accuracy of the CASSCF/RASSCF energies at the critical points on the Potential Energy Surfaces (PESs) was tested using high level electronic structure methods which include dynamical correlation. The Multi-Configurational Quasi-Degenerate Perturbation Theory through second order (MCQDPT2) method was used to correct perturbatively the SA3-CASSCF(5,6) energies. This approach was used with the same basis set as CASSCF, i.e., cc-pVDZ. In addition, the Equation-of-Motion Coupled-Cluster method with single and double excitations for the calculation of ionization potentials⁴⁵ (EOM-IP-CCSD) was used combined with the 6-311+G(d) basis set.

In the case of uracil, both the S₀ and D₀ minima (C_s symmetry) were optimized and vibrational frequencies were obtained with B3LYP/aug-cc-pVTZ, which, according to Ref. 46, gives results very close to experimental data. The IPs were calculated using EOM-IP-CCSD/6-311+G(d) in this case as well. MCQDPT2 values are taken from prior work.²⁵ The cationic states in TSH were calculated using CASSCF(13,10) with the active space including all eight π orbitals and two lone pairs on the oxygen atoms, and averaging over four states. The cc-pVDZ basis set was used with these calculations. Previous calculations have shown that this active space is sufficient to describe the low energy ionic states in the uracil cation.²⁵ From the latter reference, we also obtained the geometries of the CIs to be used here.

For all CASSCF/MRCI calculations, the basis set cc-pVDZ was used and Columbus^41,42 was employed. All B3LYP calculations were performed with Gaussian.³⁹ The EOM-IP-CCSD calculations were performed using Q-Chem,⁴⁷ while the Gamess software package⁴⁸ was used for the MCQDPT2 calculations.

The geometry data of optimized minima and minimum energy CIs can be found in the supplementary material.³⁰ 

The initial HT conformation is ccc-HT. For simplicity, we will omit the stereo description and in the remaining text HT refers to ccc-HT. The numbering of the atoms of HT and CHD can be found in Figures 1(a)-1(c). The first three cationic states of HT at the S₀ minimum correspond to holes in the three binding π orbitals. The CASSCF(5,6), MCQDPT2, and EOM-IP-CCSD configurations and state energies in comparison to the experimentally observed energies are shown in Table I, and the corresponding orbitals are depicted in Figure 2, respectively. The first IP obtained with CASSCF(5,6) is 7.55 eV. The first IPs in the EOM-IP-CCSD and MCQDPT2 calculations are 8.37 eV and 8.39 eV, respectively, in very good agreement with the experimental value of 8.45 eV given in Ref. 49. The experimental IPs were measured for ctc-HT, but the IPs of all HT isomers are similar⁴⁹ and, thus, the ctc-HT values can be used as a reference. Even though CASSCF does not reproduce the IP very well, it does a much better job at describing the relative energies between ionic states. As can be seen in Table I, the CASSCF excitation energies are close to the EOM-IP-CCSD, MCQDPT2, and experimental values.

For the CHD isomer, the most important orbitals participating in ionic states are the four π orbitals originating from the double bonds as well as σ orbitals. The first two cationic states at the S₀ minimum are holes in the π orbitals while D₂ and D₃ are produced by electrons removed from σ orbitals. The orbitals are shown in Figure 3 while the energies calculated at various levels of theory along with experimental values are shown in Table II. The experimental first IP is 8.25 eV,⁵⁰ while EOM-IP-CCSD predicts a value of 8.21 eV, MCQDPT2 8.35 eV, and CASSCF 7.77 eV. As in HT, the methods involving dynamical correlation give very accurate IPs while CASSCF has an error of 0.5 eV. The (5,6) active space of CHD includes all the π orbitals as well as σ₂₋₃ and $σ2−3*$ (orbitals located at the C²—C³ bond, see Figure 3) in order to represent the first three states. This active space gives the correct qualitative ordering of states but the energy of the D₂ state is 1.49 eV above the experimental value.

This error in D₂ cannot be improved even when employing MRCIS calculations and going beyond that level is not possible in the TSH studies. Furthermore, the σ and σ^* orbitals located at the C⁵—C⁶ bond are the orbitals necessary to represent the ring-opening from CHD to HT or the reverse process. So, the σ orbitals responsible for the ionic states are different from the ones responsible for the reaction, making the (5,6) active space inadequate for our purposes. Including both σ₂₋₃ and σ₅₋₆ with the corresponding antibonding orbitals to an (7,8) active space is still not adequate and yields too high energies for D₂. Even adding perturbation theory after CASSCF is not adequate as can be seen from the energy of D₂ at the MCQDPT2 level which still has an error of 0.5 eV (see Table II). Describing properly, the D₂ state requires more dynamical correlation and larger active spaces than is possible to use within TSH. For this reason, we found an alternative approach to this problem. We found that a RASSCF calculation where we included all the σ orbitals (except the core) in the restricted active space could bypass the problem arising from the mixing of σ orbitals. The RASSCF included all the 14 σ orbitals as doubly occupied with only single excitations allowed out of them, and the four π orbitals included in the complete active space. This RASSCF calculation preserves the order of the cationic states and also yields state energies comparable to the experimental data and for these reasons is being used in the TSH studies of CHD. Table II compares the energies at the CASSCF, RASSCF levels with the EOM-IP-CCSD and demonstrates that RASSCF gives results close to the high level methods and experiment.

Figure 4 summarizes all the important points found on the global PESs of the CHD-HT cations. All energies are given at the MCQDPT2(5,6)/cc-pVDZ level, while the energies obtained at the level of theory used in the TSH dynamics are shown in parentheses. The same diagram at the CASSCF(5,6)/cc-pVDZ and EOM-IP-CCSD/6-311+G(d) levels is shown in Figures S1 and S2.³⁰ The state energies are given relative to the minimum of the CHD D₀ state.

Several two-state CIs between D₁/D₀ as well as D₂/D₁, and one three-state CI involving D₀, D₁, and D₂ have been optimized. The two-state CIs between states I, J are labeled according to the involved electronic states DIJCI (D01CI and D12CI). The C⁵—C⁶ distances (denoted as R) and C¹—C²—C³—C⁴ dihedral angles (denoted as ϕ) of all the CIs are listed in Table III. These two geometrical parameters were chosen since they represent the ring opening/ring closing process (R) and the isomerization of HT towards its trans isomer (ϕ). Examining these parameters, we can categorize the CIs. Two CIs have structures between CHD and HT; those have the subscript “CHD”. The active space for those CIs contains two π, two π^*, the σ₅₋₆, and the $σ5−6*$ orbital. Those CIs that account for the opening of the ring beyond HT in the direction of isomerization around C²—C³ have the subscript “HT”. Furthermore, a CI with intermediate structure is found and labeled D12CI_M.

Possible pathways after ionization of HT are (a) ring closure and (b) isomerization to one of the other conformations of HT. The latter includes rotation around the various bonds of HT to form cct-HT, ctc-HT, or ttt-HT. The cct-HT and ctc-HT local minima have energies ca. 0.3 eV below ccc-HT while ttt-HT is the global minimum and its energy is 0.67 eV below ccc-HT.

The CASSCF energies (shown in parentheses in Figure 4 for the geometries of HT, to the right of the arrow) used in the dynamics differ from the MCQDPT2 on average by 0.2 eV with the maximum difference being 0.3 eV. The RASSCF values describing the CHD dynamics (shown in parentheses in Figure 4 for the geometries of CHD to the left of the arrow) give somewhat larger errors, with an average of 0.3 eV and a maximum of 0.5 eV. The largest differences occur for D12CI_CHD, and the state D₂ at the D₀ minimum. These differences may have some effects in the dynamics of CHD, but RASSCF still is the best method to balance accuracy and efficiency for this system.

In previous studies, we have investigated the ionic states of uracil and CIs between them.²⁵ The electronic structure and energies of the ionic states were calculated accurately using MRCI, MCQDPT2,²⁵ and EOM-IP-CCSD²⁶ methods. These energies can be compared with the current CASSCF approach used for the dynamics. We have also studied the dynamics and photoelectron spectrum using MCTDH.²⁶ The detailed previous studies provide an excellent way to benchmark the current work. Specifically, the photoelectron spectrum produced with MCTDH was in very good agreement with the experimental spectrum. As we will see, the present TSH dynamics produce similar results as MCTDH, indicating that the current approach represents the dynamics well.

Initially, the previously found geometries are used here and the energies are recalculated using the electronic structure methods applied in this work. The reference geometry, the S₀ minimum, and its vibrational normal modes are taken from Ref. 26. The first four cationic states of uracil correspond to electron holes in π_C=C (A″), n_O,2 (A′), π_NO (A″), and n_O,1 (A′) (see Figure 5). The atom numbering used in the following is shown in Figure 1(d). The corresponding state energies obtained with SA4-CASSCF(13,10) and EOM-IP-CCSD at the S₀ minimum as well as the experimentally measured values⁵¹ can be found in Table IV. The previously obtained MCQDPT2 energies are also shown.²⁵ EOM-IP-CCSD/6-311+G(d) gives the best agreement with the experimental results, with the first IP being only 0.1 eV higher than experiment and the other IPs having even smaller errors.

The structures that could be reached after ionization are depicted with their corresponding state energies in Figure 6. These are the S₀ minimum, the D₀ minimum and the two-state (2CI, between D₀/D₁) and a three-state CI (3CI, between D₀/D₁/D₂), whose geometries are taken from Ref. 25. Similarly to the CHD-HT system, there are two- and three-state CIs present in this system which will affect the dynamics as will be discussed below. Figure 6 shows the energies at the MCQDPT2 level while in parentheses are the CASSCF energies which are used in the TSH dynamics. The average difference between the two methods is 0.2 eV. The main difference appears to be that CASSCF underestimates the gap between D₁ and D₂ at the S₀ minimum. Another difference is the energetic difference between the energy of D₂ upon ionization and the 3CI. This difference is only 0.01 eV at the CASSCF level, while it is 0.51 eV at the MCQDPT2 level. Quite interestingly, however, these differences do not seem to have a big impact on the dynamics, and the population decays are similar to previous studies using higher level of theory, as will be discussed below.

We now discuss the results obtained from the time dependent treatment of these systems. We will first examine whether any reaction/fragmentation takes place during the dynamics. Although this is a basic question, we have to be careful since it can only be addressed in a limited way in this work. In the CHD-HT system, we can explore whether the conversion from one form to the other takes place during the dynamics and if it does whether it happens on the excited state. In uracil, fragmentation involves breaking of the ring which cannot be described properly in our studies since the active space does not include σ antibonding orbitals. Furthermore, since the low excited states are valence and not dissociative, it is unlikely that fragmentation will occur during the 120 fs of the dynamics we explored. So we only address this question for CHD-HT here while the main conclusions on relaxation are described in Sec. III E.

We first examine dynamics of HT. The reaction that can occur here is closing of the ring to form CHD. Alternatively, the cis form of HT can convert into the trans forms which are energetically more stable. In order to examine how different excited states may react we examined two sets of trajectories. Set1 includes trajectories launched from D₁ while Set2 includes trajectories from D₂. In both cases, CASSCF(5,6) calculations were averaged over the first three states. To determine whether isomerization or ring closure takes place, the angle ϕ and distance R (see Figure 1) were followed along the trajectories. Ring closure was considered when R was smaller than 2 Å and isomerization for values above 4.5 Å since the isomerization includes a change of the dihedral angle but also a huge change in the C⁵—C⁶ distance. Employing the distance as indicator for the isomerization means that isomerization also includes the rotation around all the bonds, not only around C⁵—C⁶. The two parameters for the optimized CIs are listed in Table III.

In Figure 7, the evolutions of R and ϕ are depicted. Trajectories that isomerize to trans HT are colored orange, those in which the ring closes are green, and the remaining ones in which no notable conformational change takes place (i.e., only vibrations) are blue. When starting on D₁ only six trajectories show ring closure, with one of them only temporary closing the ring and another failing before 300 fs. Out of the isomerization trajectories, those that show a permanently increasing dihedral angle perform isomerization around C²—C³ with values above 180°. Angles up to 400° indicate that more than a full rotation occurs.

When starting from D₂, the majority of the trajectories undergo an isomerization to ctc-HT, in contrast to the trajectories from D₁ where the number of trajectories with isomerization is much lower than those showing just vibrational motion. However, the number of trajectories with ring closure is very small as well, with only two trajectories forming CHD in this case.

In CHD, most of the trajectories only perform different vibrational motions. Only six trajectories elongate the C⁵—C⁶ bond toward the formation of HT, three of which fail before 300 fs (as can be seen in Figure 7(f)). The increase of the C⁵—C⁶ distance indicative of the opening reaction occurs close to D₁ → D₀ jumps and the bond breaking would then eventually occur on D₀. D₁ → D₀ jumps occur near the seam related to D01CI_CHD.

In summary, our studies show that the conversion between CHD-HT is not very common in these dynamics. Up to 6/150 (4%) trajectories in HT produced CHD and 6/186 (3%) trajectories in CHD opened the ring.

We now focus on the relaxation dynamics and examine how the population decays from the initial state to the others.

Figure 8(a) shows the adiabatic state populations for HT when dynamics are initiated on D₁. The corresponding accumulated jumping probabilities can be found in Figure S5(a).³⁰ It can be seen that initially there are a few hops from D₁ to D₂ until around 25 fs. After this time, the dynamics are determined by hops D₁ → D₀, i.e., the population in D₀ decreases about the same amount the population in D₁ increases. At 300 fs, the population in D₀ is almost 1 and D₁ is depleted. The D₂ state is only marginally populated. However, the most important part of the dynamics takes place in the first 100 fs.

The adiabatic state populations and the accumulated jumping probabilities when the initial population of HT is on D₂ are shown in Figure 8(b) and Figure S5(b),³⁰ respectively. D₂ is again depopulated fast after the start of the trajectory by jumps D₂ → D₁. After around 30 fs D₁ → D₀ hops occur with rapidly increasing frequency. Therefore, D₁ only shows an intermediate population with a peak of almost 60% at around 60 fs. At 300 fs D₂ is almost completely depopulated, about 20% of population is in D₁ and 80% in D₀, although there is still population dynamics observed at the end of the simulation time.

In order to determine the modes that are actively participating in the dynamics, we performed NMA (shown in Figures S7 and S8 of the supplementary material³⁰). This analysis indicates that the normal modes that are more active initially are the ones involving C = C and C—C stretching as well as CCH bending motions. At longer times, lower frequency modes become more important.

We can gather more information about the motion at the time immediately after 0 by examining important gradients in the Franck-Condon (FC) region. The first important gradient to consider is the gradient of the D₂ state which determines the force applied to the molecule. This vector at the S₀ minimum geometry in mass weighted coordinates is shown in Figure S13³⁰ and is denoted as g₂ here. Besides g₂, there are two vectors that are important in describing conical intersections and nonadiabatic transitions between states; these are the difference between the gradients of the two adiabatic states involved in the transition (g₁₂) and coupling vector h₁₂. h₁₂ is defined as 〈Ψ₁|∇H^e|Ψ₂〉, where Ψ_I is the electronic wavefunction for state I and H^e is the electronic Hamiltonian, and is related to the nonadiabatic (derivative) coupling by f₁₂ = h₁₂/ΔE₂₁ for exact wavefunctions.⁵² These two vectors evaluated at the S₀ minimum for HT are also shown in Figure S13. Overall, all three vectors involve motion of the C = C and C—C bonds and CCH angles, justifying the modes observed in NMA and indicating that these modes are important for initial motion away from the FC region and for nonadiabatic transitions.

The state populations for CHD are plotted in Figures 8(c). Accumulated jumping probabilities are shown in Figure S5(c).³⁰ The decay from D₂ to D₁ is extremely fast, much faster than for HT. The population on D₂ is removed in less than 10 fs. Decay from D₁ however is much slower. At 300 fs, about 10% of the population is still on D₁ while the remaining is on the ground state.

NMA (shown in Figures S9 and S10 of the supplementary material³⁰) shows that normal modes involving the C—C and C = C vibrations are very active in the initial dynamics, and especially C²—C³ and C⁵—C⁶. The dynamics starts on the D₂ state which is a state with a hole primarily on the σ₂₋₃ orbital, so the participation of the corresponding bond is expected.

More insight can be obtained by examining the g₂, g₁₂, and h₁₂ vectors, shown in Figure S14.³⁰ g₂ involves primarily C²—C³ stretching in agreement with the character of the state, and again this explains the importance of this mode in the dynamics. g₁₂ involves this stretching but also the other CC bonds are involved. The coupling involves more motion along the C = C bonds. In CHD, the CI between D₂/D₁ (D12CI_CHD) is very close to the FC region, and the g₁₂, and h₁₂ vectors shown in Figure S14 are also representative of the ones at the CI.³⁰ The close proximity of the CI and the rapid dynamics in this system lead to the FC modes being mainly responsible for the nonadiabatic transitions and decay.

Finally, uracil is described in Figures 8(d) and 9. In uracil, the decay of the D₂ population is intermediate between HT and CHD. The population decreases very fast and is totally depleted at around 60 fs. Most of the initial jumps go to D₁. This can also be seen by fast increase of D₂ → D₁ hops in the first 20 fs. However, there are also many back hops that significantly reduce the number of net D₂ → D₁ hops. After 10 fs, there is a fast increase of the number of D₁ → D₀ and D₂ → D₀ hops, so that there is a very fast increase of population in D₀ and the D₁ population has a maximum of only 30%. There is a certain number of trajectories remaining in D₁, and its population decreases only very slowly after 20 fs. After 50 fs, the accumulated jumping probabilities of D₂ → D₁ and D₂ → D₀ stay constant, with only D₁ → D₀ jumps still increasing. At 120 fs, there is 80% of the population in D₀ and still 20% in D₁. A very interesting observation here is that unlike CHD and HT in uracil there are many D₂ → D₀ jumps. In fact the accumulated jumping probability from D₂ → D₀ reaches more than 30%. Furthermore, there are many back hops between D₀ and D₁ after 15 fs, see Figure 9(b).

NMA (shown in Figures S11 and S12 of the supplementary material³⁰) shows that initially normal modes involving CN vibrations play the most important role, along with CO stretches. Again, the g₂, g₁₂, and h₁₂ vectors are shown in Figure S15. g₂ and g₁₂ are in plane vectors and they involve CN and CO vibrations. The coupling h₁₂ is an out of plane motion since it couples states of different symmetries. The in plane motions found here involve similar motions to the ones involved in the g vectors. There is one a″ mode that is actively participating in the dynamics as well. The g₁₂ and h₁₂ vectors in uracil have also been published previously at conical intersections, and they are qualitatively similar to the vectors in the FC region.²⁵ Since a substantial amount of D₂ → D₀ hops occur in uracil, the g₀₂, and h₀₂ vectors corresponding to these transitions are also shown in Figure S14.³⁰ The main difference here is that the coupling is an in plane motion since it couples two states of the same symmetry.

The dynamics on the uracil radical cation were studied by us previously using the MCTDH approach and a vibronic coupling Hamiltonian model including 10 degrees of freedom.²⁶ Quite interestingly the population decay is very similar between the two approaches, MCTDH and TSH. Figure 10 shows the populations obtained from the two approaches, and one can see the similarities in the overall trends. Dynamics of the uracil radical cation using the MCTDH approach have also shown the importance of direct D₀/D₂ coupling. This indicates that the main features of the dynamics do not change much when more than 10 degrees of freedom are included. It also shows that the semiclassical TSH approach is sufficient to predict the main effects of the dynamics. The MCTDH work predicted the photoelectron spectrum of uracil very well, providing further support that the simulated dynamics reproduce the actual dynamics in the cation accurately. The comparison between TSH and MCTDH gives us further confidence in analyzing the current dynamics.

A feature that is different between the two plots in Figure 10 is the initial decay close to t = 0. Although the D₂ decay in TSH shows an exponential decay, this is not true for MCTDH initially. The initial behavior has a Gaussian shape while later it becomes exponential. MCTDH uses only 10 degrees of freedom so the motion stays rather coherent initially, while the addition of the remaining bath degrees of freedom increases dissipation leading to the exponential decay. In the extreme case of very few degrees of freedom with bound potentials that are directly involved in the dynamics, the population will oscillate. This was shown in our previous publication using MCTDH where we showed population decays with three up to 10 degrees of freedom.²⁶ 

The modes identified here to be actively involved in the dynamics were included in the MCTDH study as well. The six most important modes used in MCTDH are three modes involving CN stretches, a C = C stretch and the two carbonyl stretches. In that case, the important modes had to be identified prior to the dynamics so that the reduced dimensionality PES would be built. In the current study, there is no a priori assumption about the important modes, and rather analysis of the results identifies the important modes. This is an advantage of TSH compared to reduced dimensionality models.

It should be pointed out that the two studies differ in the underlying electronic structure theory in addition to the dynamical approach. The MCTDH study used EOM-IP-CCSD/6-311+G(d) while the current TSH study uses CASSCF/cc-pVDZ. There are some differences in the energetics predicted by these two methods, but based on the results these differences do not significantly affect the dynamics. The differences can be seen in the FC region where the D₂-D₁ gap is only 0.05 eV at the CASSCF level but 0.36 eV at the EOM-IP-CCSD level. Nevertheless decay from D₂ occurs similarly in the two studies. Another difference is the energy of the 3CI, which is 0.01 eV below the D₂ energy at S₀ minimum and 0.08 eV above the D₂ energy at the EOM-IP-CCSD level. Nevertheless D₂-D₀ transitions occur in both cases. As will be discussed below, the reason may be that the coupling plays a bigger role than the exact energy of 3CI.

Although the dynamics are fast for all the cations studied here, there are several differences in the details of how they occur. Can we explain the differences in the dynamics by some of the characteristics of the states? It would be very informative if we could use the energy gap between the states or the coupling or another property to guide us in predetermining the rates for decay. We are interested to see whether there are any correlations of the properties of the molecule initially and the rate of decay, keeping in mind that we only have three molecules, thus a very small number of cases to make general statements. Nevertheless, some interesting observations can be made. We focus on the decay of the D₂ population, and we use the time τ₂ at which the population is 1/e as the quantitative measure of the dynamics. We attempt to correlate that time to various properties. In uracil, we use properties for both the D₂/D₁ and the D₂/D₀ transitions since it is evident in Figure 9 that there is a significant amount of D₂/D₀ hops. The decay of D₂ will depend on both of these transitions so we take into account both of them in the following discussion.

The first intuitive idea is that the gap between states is important in the nonadiabatic rate. The gap between D₂-D₁ (ΔE₂₁) at the initial S₀ geometry is shown in Table V. The gap is much smaller in uracil compared to CHD and HT (at least an order of magnitude), yet the decay from D₂ to D₁ is faster for CHD than for uracil. So clearly the gap cannot explain the observed dynamics.

Another factor contributing to fast decay, especially in fast decay initially to move away from the FC region, is the initial slope of the state. The magnitude of the gradient of the D₂ state in mass weighted coordinates at the S₀ minimum geometry (g₂, shown in Table V) is higher for CHD correlating with the fastest decay in this system, but it does not correlate well with the decay for the other two molecules since it predicts that HT would be somewhat faster than uracil.

Another crucial property for nonadiabatic transitions is the coupling between the two states, which is examined here by the magnitude of the derivative coupling at the S₀ minimum geometry. Table V shows the magnitude of the derivative coupling f₁₂, where we have approximated it with f₁₂ = h₁₂/ΔE₂₁. Although this relation is exact for exact eigenfunctions, it only gives the dominant CI component of f₁₂ for CASSCF wavefunctions used here.^52,53 The derivative coupling is maximum for uracil, again not correlating with the observed rates.

As discussed above, two critical vectors in nonadiabatic transitions are the gradient difference g₁₂ and the coupling h₁₂ (which is related to f₁₂). In the surface hopping approach, the nonadiabatic transitions are proportional to the dot product of the derivative coupling vector and the velocity vector. For simplicity here, the velocity can be approximated by the initial force (given by g₂). Then the dot product of g₂ with the derivative coupling (f₁₂ ⋅ g₂) is calculated to give the overall effect of how the direction of f₁₂ is coupled to the initial momentum. Similarly, we can also take the dot product between g₁₂ and g₂. These two dot products are shown in Table V.

The dot product that correlates best with the observed rates is the g₁₂ ⋅ g₂. This has a much larger value for CHD correlating with the much faster decay. In uracil, the decay depends on both D₂ → D₁ and D₂ → D₀ transitions. The dot products for these two transitions are quite different with D₂ → D₀ being more effective, so a combination of these two values (1.4 and 46) will be responsible for the observed decay constants. HT has the smallest value of the dot product correlating with the slowest decay.

The dot product f₁₂ ⋅ g₂ is also very small for HT, again in agreement with a slow decay. In uracil, this dot product shows very different values for the D₂ → D₁ and D₂ → D₀ transitions. The coupling D₁/D₀ is perpendicular to the plane of the molecule, giving a very small dot product with the initial momentum which is in plane. Physically that means that vibrational motion has to be redirected from in plane to out of plane for the transition to occur making the transitions inefficient. On the other hand, the derivative coupling has a very high overlap with g₂ for the D₂/D₀ transition in uracil since both vectors are in plane and strongly coupled. This highlights once again the importance of symmetry in uracil. One could attempt to generalize that this redirection of vibrational motion will be a problem in transitions involving states of different symmetries in general. The dynamics in uracil are driven by an interplay of many factors. Even though the coupling is much larger for the D₂ → D₀ transition the gap is very small for the D₂ → D₁ transition, and the dynamics will depend on a balance between those two. As Figure 9 shows, more hops occur between D₂ → D₁ rather than D₂ → D₀, indicating that the very small gap is important.

The final conclusion is that examining properties initially (at the S₀ minimum) can give some information about the relative magnitudes of decay rates in molecules similar to the ones we study here. These correlations seem to work here because the conical intersections are not far away from the FC region and the dynamics are very fast. In cases where the CI is very far away from the FC region examining properties in the FC region will probably not work as well. Intuitive ideas, such as the expectation that a small energy gap between electronic states will lead to the fastest decay, may lead to wrong predictions. Examining both the magnitude and the direction of the coupling and the gradient difference vectors and how they are related to the slope of the PES provides a much better insight.

Another interesting aspect of the calculations to analyze is how three-state deactivation occurs when dynamics start on D₂. In all three systems, we started the dynamics on D₂ and deactivation to the ground state requires coupling to D₁ and to D₀. How this deactivation occurs however is quite different between the systems. For HT and CHD, only jumps between neighboring states can be observed, whereas uracil clearly allows for direct jumps D₂ → D₀ especially in the first 50 fs. This allows for a faster deactivation than in HT, although CHD relaxes even faster than uracil even without D₂ → D₀ hops. Naturally, the question that arises is what determines the involvement of direct D₂ → D₀ hops.

The magnitudes of the couplings h_IJ at the S₀ minimum between all pairs of states for the three systems are shown in Table VI. The couplings are similar in magnitude between all pairs of states for CHD and HT. In uracil, however, the D₀/D₂ coupling is an order of magnitude larger than the D₀/D₁ and D₁/D₂ ones. This leads to a very favorable coupling directly from D₂ to D₀ and that is a major reason why hops occur directly between these states.

The presence and energetics of the three-state CIs and the region around them are also important. In the case of uracil, the three-state minimum CI is 0.81 eV above D₀ at the S₀ minimum, but also energetically in between D₁ and D₂. The geometry of the three-state CI does not deviate from the equilibrium geometry of uracil dramatically. For these reasons, the region around 3CI is easily accessible. Although there is a three-state CI that could be accessed in HT, this is 1.48 eV above D₀ at the S₀ minimum geometry, 0.59 eV below D₁, 1.5 eV below D₂, and therefore, not as near to D₂ as in uracil. The energy of the HT 3CI is close to D01CI_HT and its geometry is closer to D12CI_HT which poses a larger geometrical change with respect to the S₀ minimum of HT.

Figure 11 shows the energy gap D_I–D_J when the D_I → D_J hops occur in uracil. The D₂–D₁ and D₁–D₀ hops occur close to the seam as the energy gap suggests, while a small percent occurs at larger gaps. The D₂–D₀ hops occur further away from the seam. This is reasonable since the space where the seam occurs (3N − 6 − 5 for a system with N atoms) is much smaller than the conformational space available to the molecule (3N − 6), so statistically it is more likely to not occur at the seam. The difference in the magnitude of the couplings again has an effect here. Since h₀₂ is about an order of magnitude larger than h₁₂ and h₀₁, it will lead to a larger derivative coupling f₀₂ = h₀₂/(E₂ − E₀) even for relative larger gaps compared to the situation in D₁ → D₀ and D₂ → D₁. In general, the combination of very strong coupling between D₀/D₂ and the accessibility of the region around the 3CI lead to the importance of D₂ → D₀ hops in uracil. The coupling is strong enough that the molecule does not need to access the seam of three-state conical intersections.

Nonadiabatic dynamics of three radical cations (CHD, HT, and uracil) starting from an excited ionic state have been examined using trajectory surface hopping. The aim of this study was primarily to examine how radicals that we have previously investigated using strong fields behave after generated on an excited state and to test assumptions often made in pump probe techniques where ions are observed. We specifically wanted to examine the applicability of fast radiationless decay before dissociation in ionic states, as well as whether there are general rules describing ionic dynamics. Even though the three radical cations studied here cannot represent every type of radical cation, there is one common characteristic that we want to highlight: Here, and in most radicals, the energy gaps between all states including the ground state are quite small, while in most neutral molecules, the ground state is well separated from the excited states. As a result radiationless decay to the ground state is expected to be much faster in radicals.

The results indicate that relaxation to the ground ionic state is indeed very fast in these systems, while fragmentation before relaxation is rare. Based on these findings, it is reasonable to assume that fast radiationless decay to the ground state can be expected in radical cations when there are no directly dissociative states involved. Even though the systems studied here did not involve directly dissociative states, we can assume that in that case dissociation may precede relaxation.

While the overall decay is always fast, the details of the dynamics and the lifetimes of the states depend on the specific systems. Examining the properties of the systems in the FC region can give some insight in the subsequent dynamics. Several factors play a role in the nonadiabatic transitions, and focusing only on one property, such as the energy gap between electronic states, may lead to wrong predictions. The gradient of the PES as well as the coupling between states are also important. The best insight however is provided when examining both the magnitude and the direction of the g and h vectors and how they are related to the slope of the PES, since this shows whether the initially excited modes can be projected onto the coupling facilitating transitions.

While in HT and CHD only transitions between neighboring states dominate the dynamics, in uracil direct decay from D₂ to D₀ is very important and contributes significantly to the decay of D₂ population. One reason for this direct decay is the presence of an energetically accessible three-state CI seam, but even more important is the fact that the coupling between D₂ and D₀ is much larger than the couplings between neighboring states (D₁/D₀ and D₂/D₁).

We gratefully acknowledge the support from the Department of Energy under Award Nos. DE-FG02-08ER15983 and DE-FG02-08ER15984. We thank Felix Plasser for help in setting up the trajectories with Newton-X. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1053575.