Ionization processes can lead to the formation of radical cations with population in several ionic states. In this study, we examine the dynamics of three radical cations starting from an excited ionic state using trajectory surface hopping dynamics in combination with multiconfigurational electronic structure methods. The efficiency of relaxation to the ground state is examined in an effort to understand better whether fragmentation of cations is likely to occur directly on excited states or after relaxation to the ground state. The results on cyclohexadiene, hexatriene, and uracil indicate that relaxation to the ground ionic state is very fast in these systems, while fragmentation before relaxation is rare. Ultrafast relaxation is facilitated by the close proximity of electronic states and the presence of two- and three-state conical intersections. Examining the properties of the systems in the Franck-Condon region can give some insight into the subsequent dynamics.

## I. INTRODUCTION

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^{9} for neutral excited states where fast relaxation from higher excited states to S_{1} is assumed followed by emission from the S_{1} 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_{0} to examine fragmentation of the uracil radical cation in the past, assuming that barriers on the ground state D_{0} are responsible for the dissociation pathways.^{10} 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_{0}? Is ultrafast radiationless decay to D_{0} 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.^{19} The excitation of CHD to its S_{1} 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.^{3} 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.^{3}

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.^{24} In order to interpret the results, we needed a good understanding of the ionization and fragmentation steps, and we investigated those in detail.^{10} An assumption that we made was that fragmentation of the cation occurs after relaxation to D_{0}. In addition, we have examined the potential energy surfaces of several states in the cation including CIs,^{25} and we have simulated dynamics and spectra using the Multiconfigurational Time-Dependent Hartree (MCTDH) method.^{26} 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.

## II. METHODOLOGY

### A. Trajectory surface hopping

The reference geometry for all TSH simulations is the corresponding optimized neutral ground state (S_{0}) 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.^{30} All trajectories were run using the local diabatization method.^{31} 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.^{30} 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.^{34} 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_{2} adiabatic state in all systems, except HT where we have also run a set of trajectories starting from D_{1}. For CHD and uracil, ionization simulations and experiments indicate that D_{2} 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_{1} and D_{2} 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^{30} 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_{0} minimum with a least-squares fit method, also part of the *Newton-X* module.^{38}

### B. Quantum chemistry calculations

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 *Gaussian*03 program.^{39} Both isomers possess C_{2} symmetry in their S_{0} and D_{0} minima.

The location of CIs and the TSH dynamics were carried out using the complete active space self-consistent field method^{40} (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^{43} (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^{44} 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^{45} (EOM-IP-CCSD) was used combined with the 6-311+G(d) basis set.

In the case of uracil, both the S_{0} and D_{0} 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.^{25} 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.^{25} 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*.^{39} The EOM-IP-CCSD calculations were performed using *Q-Chem*,^{47} while the *Gamess* software package^{48} was used for the MCQDPT2 calculations.

The geometry data of optimized minima and minimum energy CIs can be found in the supplementary material.^{30}

## III. RESULTS AND DISCUSSION

### A. Electronic structure of the HT-CHD system

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_{0} 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^{49} 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.

State (Sym) . | Character . | CASSCF . | MCQDPT2 . | EOM-IP-CCSD . | Expt.^{49}
. |
---|---|---|---|---|---|

D_{0} (B) | π_{3} | 0 | 0 | 0 | 0 |

D_{1} (A) | π_{2} | 2.07 | 1.83 | 1.95 | 1.97 |

D_{2} (B) | π_{1} | 2.98 | 2.85 | 3.00 | 3.31 |

State (Sym) . | Character . | CASSCF . | MCQDPT2 . | EOM-IP-CCSD . | Expt.^{49}
. |
---|---|---|---|---|---|

D_{0} (B) | π_{3} | 0 | 0 | 0 | 0 |

D_{1} (A) | π_{2} | 2.07 | 1.83 | 1.95 | 1.97 |

D_{2} (B) | π_{1} | 2.98 | 2.85 | 3.00 | 3.31 |

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_{0} minimum are holes in the *π* orbitals while D_{2} and D_{3} 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,^{50} 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 *σ*_{2−3} and $\sigma 2\u22123*$ (orbitals located at the C^{2}—C^{3} 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_{2} state is 1.49 eV above the experimental value.

State (Sym.) . | Character . | CASSCF(5,6) . | RASSCF . | MCQDPT2 . | EOM-IP-CCSD . | exp^{50}
. |
---|---|---|---|---|---|---|

D_{0} (A) | π_{2} | 0 | 0 | 0 | 0 | 0 |

D_{1} (B) | π_{1} | 2.56 | 2.65 | 2.55 | 2.71 | 2.45 |

D_{2} (A) | σ_{2} (σ_{2-3})^{a} | 4.49 | 3.15 | 3.51 | 2.92 | 3.0 |

D_{3} (B) | σ_{1} | 3.69 | 3.41 | 3.55 |

State (Sym.) . | Character . | CASSCF(5,6) . | RASSCF . | MCQDPT2 . | EOM-IP-CCSD . | exp^{50}
. |
---|---|---|---|---|---|---|

D_{0} (A) | π_{2} | 0 | 0 | 0 | 0 | 0 |

D_{1} (B) | π_{1} | 2.56 | 2.65 | 2.55 | 2.71 | 2.45 |

D_{2} (A) | σ_{2} (σ_{2-3})^{a} | 4.49 | 3.15 | 3.51 | 2.92 | 3.0 |

D_{3} (B) | σ_{1} | 3.69 | 3.41 | 3.55 |

^{a}

In the CASSCF(5,6) calculation D_{2} is a hole in the *σ*_{2−3} orbital that is depicted in Figure 3(h).

This error in D_{2} 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^{5}—C^{6} 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 *σ*_{2−3} and *σ*_{5−6} with the corresponding antibonding orbitals to an (7,8) active space is still not adequate and yields too high energies for D_{2}. Even adding perturbation theory after CASSCF is not adequate as can be seen from the energy of D_{2} at the MCQDPT2 level which still has an error of 0.5 eV (see Table II). Describing properly, the D_{2} 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.

### B. Important features on the PESs of the HT-CHD system

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.^{30} The state energies are given relative to the minimum of the CHD D_{0} state.

Several two-state CIs between D_{1}/D_{0} as well as D_{2}/D_{1}, and one three-state CI involving D_{0}, D_{1}, and D_{2} have been optimized. The two-state CIs between states *I*, *J* are labeled according to the involved electronic states D*IJ*CI (D01CI and D12CI). The C^{5}—C^{6} distances (denoted as R) and C^{1}—C^{2}—C^{3}—C^{4} 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 *σ*_{5−6}, and the $\sigma 5\u22126*$ orbital. Those CIs that account for the opening of the ring beyond HT in the direction of isomerization around C^{2}—C^{3} have the subscript “HT”. Furthermore, a CI with intermediate structure is found and labeled D12CI_{M}.

Structure . | R . | ϕ
. |
---|---|---|

D01CI_{CHD} | 2.30 | 15.8 |

D12CI_{CHD} | 1.78 | 10.4 |

D01CI_{HT} | 4.53 | −2.0 |

D12CI_{HT} | 4.74 | −0.5 |

D12CI_{M} | 3.54 | 23.4 |

3CI^{a} | 4.87 | −1.2 |

Structure . | R . | ϕ
. |
---|---|---|

D01CI_{CHD} | 2.30 | 15.8 |

D12CI_{CHD} | 1.78 | 10.4 |

D01CI_{HT} | 4.53 | −2.0 |

D12CI_{HT} | 4.74 | −0.5 |

D12CI_{M} | 3.54 | 23.4 |

3CI^{a} | 4.87 | −1.2 |

^{a}

This geometry is optimized at the MRCIS level.

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_{2} at the D_{0} 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.

### C. Electronic structure and important features on the PESs of uracil

In previous studies, we have investigated the ionic states of uracil and CIs between them.^{25} The electronic structure and energies of the ionic states were calculated accurately using MRCI, MCQDPT2,^{25} and EOM-IP-CCSD^{26} 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.^{26} 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_{0} 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_{0} minimum as well as the experimentally measured values^{51} can be found in Table IV. The previously obtained MCQDPT2 energies are also shown.^{25} 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.

. | Character . | CASSCF . | MCQDPT2^{25}
. | EOM-IP-CCSD . | Expt.^{51}
. |
---|---|---|---|---|---|

D_{0} (A″) | π_{C=C} | 0 | 0 | 0 | 0 |

D_{1} (A′) | n_{O,2} | 0.78 | 0.33 | 0.69 | 0.6 |

D_{2} (A″) | π_{NO} | 0.83 | 0.97 | 1.05 | 1.2 |

D_{3} (A′) | n_{O,1} | 1.46 | 1.45 | 1.65 | 1.65 |

The structures that could be reached after ionization are depicted with their corresponding state energies in Figure 6. These are the S_{0} minimum, the D_{0} minimum and the two-state (2CI, between D_{0}/D_{1}) and a three-state CI (3CI, between D_{0}/D_{1}/D_{2}), 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_{1} and D_{2} at the S_{0} minimum. Another difference is the energetic difference between the energy of D_{2} 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.

### D. Dynamics: Reactivity

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_{1} while Set2 includes trajectories from D_{2}. 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^{5}—C^{6} distance. Employing the distance as indicator for the isomerization means that isomerization also includes the rotation around all the bonds, not only around C^{5}—C^{6}. 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_{1} 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^{2}—C^{3} with values above 180°. Angles up to 400° indicate that more than a full rotation occurs.

When starting from D_{2}, the majority of the trajectories undergo an isomerization to ctc-HT, in contrast to the trajectories from D_{1} 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^{5}—C^{6} 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^{5}—C^{6} distance indicative of the opening reaction occurs close to D_{1} → D_{0} jumps and the bond breaking would then eventually occur on D_{0}. D_{1} → D_{0} 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.

### E. Dynamics: Radiationless decay

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

#### 1. HT

Figure 8(a) shows the adiabatic state populations for HT when dynamics are initiated on D_{1}. The corresponding accumulated jumping probabilities can be found in Figure S5(a).^{30} It can be seen that initially there are a few hops from D_{1} to D_{2} until around 25 fs. After this time, the dynamics are determined by hops D_{1} → D_{0}, i.e., the population in D_{0} decreases about the same amount the population in D_{1} increases. At 300 fs, the population in D_{0} is almost 1 and D_{1} is depleted. The D_{2} 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_{2} are shown in Figure 8(b) and Figure S5(b),^{30} respectively. D_{2} is again depopulated fast after the start of the trajectory by jumps D_{2} → D_{1}. After around 30 fs D_{1} → D_{0} hops occur with rapidly increasing frequency. Therefore, D_{1} only shows an intermediate population with a peak of almost 60% at around 60 fs. At 300 fs D_{2} is almost completely depopulated, about 20% of population is in D_{1} and 80% in D_{0}, 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^{30}). 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_{2} state which determines the force applied to the molecule. This vector at the S_{0} minimum geometry in mass weighted coordinates is shown in Figure S13^{30} and is denoted as **g**_{2} here. Besides **g**_{2}, 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**_{12}) and coupling vector **h**_{12}. **h**_{12} is defined as 〈Ψ_{1}|∇*H ^{e}*|Ψ

_{2}〉, where Ψ

_{I}is the electronic wavefunction for state

*I*and

*H*is the electronic Hamiltonian, and is related to the nonadiabatic (derivative) coupling by

^{e}**f**

_{12}=

**h**

_{12}/Δ

*E*

_{21}for exact wavefunctions.

^{52}These two vectors evaluated at the S

_{0}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.

#### 2. CHD

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

NMA (shown in Figures S9 and S10 of the supplementary material^{30}) shows that normal modes involving the C—C and C = C vibrations are very active in the initial dynamics, and especially C^{2}—C^{3} and C^{5}—C^{6}. The dynamics starts on the D_{2} state which is a state with a hole primarily on the *σ*_{2−3} orbital, so the participation of the corresponding bond is expected.

More insight can be obtained by examining the **g**_{2}, **g**_{12}, and **h**_{12} vectors, shown in Figure S14.^{30} **g**_{2} involves primarily C^{2}—C^{3} stretching in agreement with the character of the state, and again this explains the importance of this mode in the dynamics. **g**_{12} 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_{2}/D_{1} (D12CI_{CHD}) is very close to the FC region, and the **g**_{12}, and **h**_{12} vectors shown in Figure S14 are also representative of the ones at the CI.^{30} 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.

#### 3. Uracil

Finally, uracil is described in Figures 8(d) and 9. In uracil, the decay of the D_{2} 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_{1}. This can also be seen by fast increase of D_{2} → D_{1} hops in the first 20 fs. However, there are also many back hops that significantly reduce the number of net D_{2} → D_{1} hops. After 10 fs, there is a fast increase of the number of D_{1} → D_{0} and D_{2} → D_{0} hops, so that there is a very fast increase of population in D_{0} and the D_{1} population has a maximum of only 30%. There is a certain number of trajectories remaining in D_{1}, and its population decreases only very slowly after 20 fs. After 50 fs, the accumulated jumping probabilities of D_{2} → D_{1} and D_{2} → D_{0} stay constant, with only D_{1} → D_{0} jumps still increasing. At 120 fs, there is 80% of the population in D_{0} and still 20% in D_{1}. A very interesting observation here is that unlike CHD and HT in uracil there are many D_{2} → D_{0} jumps. In fact the accumulated jumping probability from D_{2} → D_{0} reaches more than 30%. Furthermore, there are many back hops between D_{0} and D_{1} after 15 fs, see Figure 9(b).

NMA (shown in Figures S11 and S12 of the supplementary material^{30}) shows that initially normal modes involving CN vibrations play the most important role, along with CO stretches. Again, the **g**_{2}, **g**_{12}, and **h**_{12} vectors are shown in Figure S15. **g**_{2} and **g**_{12} are in plane vectors and they involve CN and CO vibrations. The coupling **h**_{12} 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**_{12} and **h**_{12} vectors in uracil have also been published previously at conical intersections, and they are qualitatively similar to the vectors in the FC region.^{25} Since a substantial amount of D_{2} → D_{0} hops occur in uracil, the **g**_{02}, and **h**_{02} vectors corresponding to these transitions are also shown in Figure S14.^{30} The main difference here is that the coupling is an in plane motion since it couples two states of the same symmetry.

### F. Comparison of TSH and MCTDH dynamics

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.^{26} 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_{0}/D_{2} 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_{2} 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.^{26}

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_{2}-D_{1} gap is only 0.05 eV at the CASSCF level but 0.36 eV at the EOM-IP-CCSD level. Nevertheless decay from D_{2} occurs similarly in the two studies. Another difference is the energy of the 3CI, which is 0.01 eV below the D_{2} energy at S_{0} minimum and 0.08 eV above the D_{2} energy at the EOM-IP-CCSD level. Nevertheless D_{2}-D_{0} 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.

### G. Correlations between properties and decay times

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_{2} population, and we use the time *τ*_{2} 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_{2}/D_{1} and the D_{2}/D_{0} transitions since it is evident in Figure 9 that there is a significant amount of D_{2}/D_{0} hops. The decay of D_{2} 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_{2}-D_{1} (Δ*E*_{21}) at the initial S_{0} 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_{2} to D_{1} is faster for CHD than for uracil. So clearly the gap cannot explain the observed dynamics.

System . | τ_{2}
. | ΔE_{21}
. | g_{2}
. | f_{12}
. | |f_{12} ⋅ g_{2}|/10^{4}
. | |g_{12} ⋅ g_{2}|/10^{4}
. |
---|---|---|---|---|---|---|

HT | 46 | 0.91 | 0.04 | 2.38 | 0.4 | 4.6 |

Uracil | 12 | 0.05 | 0.03 | 6.07 (4.38^{a}) | 5.5 (125.0^{a}) | 1.4 (46.0^{a}) |

CHD | 4 | 0.50 | 0.06 | 2.50 | 10.2 | 124.5 |

System . | τ_{2}
. | ΔE_{21}
. | g_{2}
. | f_{12}
. | |f_{12} ⋅ g_{2}|/10^{4}
. | |g_{12} ⋅ g_{2}|/10^{4}
. |
---|---|---|---|---|---|---|

HT | 46 | 0.91 | 0.04 | 2.38 | 0.4 | 4.6 |

Uracil | 12 | 0.05 | 0.03 | 6.07 (4.38^{a}) | 5.5 (125.0^{a}) | 1.4 (46.0^{a}) |

CHD | 4 | 0.50 | 0.06 | 2.50 | 10.2 | 124.5 |

^{a}

Values in parentheses correspond to the D_{2}/D_{0} pair in uracil.

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_{2} state in mass weighted coordinates at the S_{0} minimum geometry (*g*_{2}, 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_{0} minimum geometry. Table V shows the magnitude of the derivative coupling **f**_{12}, where we have approximated it with **f**_{12} = **h**_{12}/Δ*E*_{21}. Although this relation is exact for exact eigenfunctions, it only gives the dominant CI component of **f**_{12} 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**_{12} and the coupling **h**_{12} (which is related to **f**_{12}). 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**_{2}). Then the dot product of **g**_{2} with the derivative coupling (**f**_{12} ⋅ **g**_{2}) is calculated to give the overall effect of how the direction of **f**_{12} is coupled to the initial momentum. Similarly, we can also take the dot product between **g**_{12} and **g**_{2}. These two dot products are shown in Table V.

The dot product that correlates best with the observed rates is the **g**_{12} ⋅ **g**_{2}. This has a much larger value for CHD correlating with the much faster decay. In uracil, the decay depends on both D_{2} → D_{1} and D_{2} → D_{0} transitions. The dot products for these two transitions are quite different with D_{2} → D_{0} 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**_{12} ⋅ **g**_{2} 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_{2} → D_{1} and D_{2} → D_{0} transitions. The coupling D_{1}/D_{0} 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**_{2} for the D_{2}/D_{0} 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_{2} → D_{0} transition the gap is very small for the D_{2} → D_{1} transition, and the dynamics will depend on a balance between those two. As Figure 9 shows, more hops occur between D_{2} → D_{1} rather than D_{2} → D_{0}, indicating that the very small gap is important.

The final conclusion is that examining properties initially (at the S_{0} 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.

### H. Three-state conical intersections and D_{2} → D_{0} transitions

Another interesting aspect of the calculations to analyze is how three-state deactivation occurs when dynamics start on D_{2}. In all three systems, we started the dynamics on D_{2} and deactivation to the ground state requires coupling to D_{1} and to D_{0}. 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_{2} → D_{0} 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_{2} → D_{0} hops. Naturally, the question that arises is what determines the involvement of direct D_{2} → D_{0} hops.

The magnitudes of the couplings *h _{IJ}* at the S

_{0}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

_{0}/D

_{2}coupling is an order of magnitude larger than the D

_{0}/D

_{1}and D

_{1}/D

_{2}ones. This leads to a very favorable coupling directly from D

_{2}to D

_{0}and that is a major reason why hops occur directly between these states.

System . | h_{01}
. | h_{02}
. | h_{12}
. |
---|---|---|---|

HT | 0.09 | 0.08 | 0.08 |

CHD | 0.10 | 0.07 | 0.05 |

Uracil | 0.02 | 0.13 | 0.01 |

System . | h_{01}
. | h_{02}
. | h_{12}
. |
---|---|---|---|

HT | 0.09 | 0.08 | 0.08 |

CHD | 0.10 | 0.07 | 0.05 |

Uracil | 0.02 | 0.13 | 0.01 |

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_{0} at the S_{0} minimum, but also energetically in between D_{1} and D_{2}. 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_{0} at the S_{0} minimum geometry, 0.59 eV below D_{1}, 1.5 eV below D_{2}, and therefore, not as near to D_{2} 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_{0} 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_{2}–D_{1} and D_{1}–D_{0} hops occur close to the seam as the energy gap suggests, while a small percent occurs at larger gaps. The D_{2}–D_{0} hops occur further away from the seam. This is reasonable since the space where the seam occurs (3*N* − 6 − 5 for a system with *N* atoms) is much smaller than the conformational space available to the molecule (3*N* − 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*_{02} is about an order of magnitude larger than *h*_{12} and *h*_{01}, it will lead to a larger derivative coupling *f*_{02} = *h*_{02}/(*E*_{2} − *E*_{0}) even for relative larger gaps compared to the situation in D_{1} → D_{0} and D_{2} → D_{1}. In general, the combination of very strong coupling between D_{0}/D_{2} and the accessibility of the region around the 3CI lead to the importance of D_{2} → D_{0} hops in uracil. The coupling is strong enough that the molecule does not need to access the seam of three-state conical intersections.

## IV. CONCLUSIONS

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_{2} to D_{0} is very important and contributes significantly to the decay of D_{2} 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_{2} and D_{0} is much larger than the couplings between neighboring states (D_{1}/D_{0} and D_{2}/D_{1}).

## Acknowledgments

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.

## REFERENCES

*ab initio*electronic structure program, release 7.0, Columbus, 2012.