The excited state dynamics of o-nitrophenol (o-NP) has been a topic of great interest because it produces atmospherically significant photofragments, such as HONO, NO and OH radicals after excitation by UV light. [1][2][3] Experimental techniques like femtosecond time resolved photoelectron and ionization spectroscopy (TRPES and TRIS) 4 and ultrafast electron diffraction (UED) 5 can be, and have been, used to study these processes. [6][7][8][9][10] Informative as these methods are, spectroscopic measurements such as TRPES cannot resolve structural changes (reaction mechanisms) and structural probes such as UED cannot resolve changes of electronic state character or changes in energy. Therefore, they must be compared with theoretical calculations to obtain a more complete picture of when and how the important photophysical and photochemical processes, such as internal conversion (IC), excited state intramolecular proton transfer (ESIPT), intersystem crossing (ISC), and fragmentation, occur. 6,11 The dynamics of o-NP after excitation have been studied numerous times before using femtosecond time resolved techniques. [6][7][8][9][10] The gas-phase absorption spectrum of o-NP has two broad bands, centered at about 331 and 259 nm. [12][13][14][15] Most of the previous experimental pump-probe studies excited o-NP within the lower absorption peak using pulses between 300-400 nm. [7][8][9] Ernst et al. performed TRPES measurements with a pump pulse at 350 nm found two very short time constants, 14 fs, 130 fs and a much longer one. 8 In another TRPES study, Nitta et al. 9 used 400 nm pump and 29.5 eV probe, and also observed a very fast decay of about 60 fs. There are two recent TRPES studies where shorter wavelength pump pulses were used, 250-270 nm, leading to excitation of the second absorption peak. 6,10 Ciavardini et al. performed TRPES using a 268 nm pump and 53 nm probe. They measured four time constants from their signals: 74 fs, 127 fs, 2.3 ps, and ≥ 10 ps. 10 In our previous work using 255 nm pump and 153 nm probe to measure a combination of TRPES and TRIS signals, we found decay time constants of about 120 fs. 6 Very recently a UED probe has also been applied to o-NP after being excited by 330 nm pulses. 7 The experimental signal, which originates from time-dependent changes in interatomic distances, shows a delay of approximately 400 fs between the onset of features assigned to nitro group rotation and time-zero. By comparing with theory (ab initio multiple spawning, AIMS), this time has been assigned to decay to the ground state through conical intersections. In summary, all these experiments measure an ultrafast decay of the measurement observable, but there are differences which require further interpretation.
In order to interpret the experimental measurements, theoretical studies are needed. Such studies have found that ESIPT, leading to the aci-nitrophenol isomer, is important in the excited state dynamics, as well as passage through conical intersections leading to the ground state. 7,[9][10][11]16,17 Most previous theoretical calculations on the mechanisms for photochemistry were initiated on the first excited state, corresponding to the first absorption peak. 7,9,11,16,17 These studies find that the dynamics is governed by ESIPT followed by internal conversion to the ground state or by intersystem crossing. Rotation of the HONO group out of the plane is an important coordinate leading to the conical intersection with the ground state. While most studies used calculations of stationary points and potential energy surface pathways, there were two studies employing on the fly excited state dynamics. Xu et al. used a global switching surface hopping algorithm where a modified Landau-Zener formula is used to calculate both internal conversion and ISC 17 . They found that most of the population goes to triplet states through ISC followed by ESIPT. Nunes et al. used ab initio multiple spawning 7 and found contradictory results. They found that ESIPT occurs on S 1 , while they determined that ISC does not take place on such fast timescales because of the small spin-orbit coupling (SOC). Other theoretical studies on the ground singlet and triplet surfaces suggest that dissociation to form NO, OH, HONO, most likely occurs on the triplet surface. 2,18 While there are many experiments pumping on either absorption band, theoretical dynamics studies have focused on the lower energy peak. 11,16,17,19 In this work, we study the dynamics starting from the high energy peak for a better comparison with the corresponding experiments. Furthermore, o-NP gives the perfect opportunity to examine the effect of excitation wavelength on the excited state dynamics. So, we also carry out dynamics starting from the lower energy peak and compare the two sets of results. We make use of trajectory surface hopping (TSH) 20 combined with floating occupation molecular orbital-complete active space configuration interaction (FOMO-CASCI) 21 . As implemented in TeraChem (v1.9) 22 , this method offers extremely fast single point calculations by leveraging the power of graphical processing units (GPUs) that allows the simulation of a large number of trajectories for long timescales in a short amount of time. Furthermore, because FOMO-CASCI is a multi-reference method, it can describe radiationless decay to the ground state as well as bond dissociation qualitatively accurately.
The FOMO-CASCI method 21 is used in this work in order to be able to run dynamics for longer times. FOMO-CASCI has been tested before for nonadiabatic dynamics using TSH and has been found to perform similar to Complete Active Space Self-Consistent Field Theory (CASSCF) while being much faster. 23 In FOMO-CASCI, orbitals are obtained from a finite-temperature self-consistent field calculation (or fractional occupation number SCF, FON-SCF) and then used in CASCI. In the FON-SCF method, the occupation number of a given orbital i changes continuously according to the energy of the orbital ε i . TeraChem implements FOMO-CASCI using a Gaussian broadening scheme:
where the occupation number n i of orbital i is obtained by the integration of the distribution function up to the Fermi level, ε f ,
The broadening coefficient β, which is analogous to an electronic temperature, is a parameter that affects the performance of the methodology. It is crucial that FOMO-CASCI accurately captures the energies and oscillator strengths of the first four singlet excited states in nitrophenol to run accurate relaxation dynamics on them. To tune the accuracy of FOMO-CASCI, the optimal broadening parameter, basis set, and active space must be chosen.
We benchmarked FOMO-CASCI against CASSCF 24 and multistate complete active space second order perturbation theory (MS-CASPT2). 25,26 CASSCF with an active space of 14 electrons in 10 orbitals with the cc-pVDZ basis set, averaging over 6 singlet states, SA6-CASSCF(14,10)/cc-pVDZ, was chosen as the standard to benchmark FOMO-CASCI. The orbitals for both CASSCF and FOMO-CASCI are shown in Supplementary Material (SM) (SM: Figure 2). Additional MS-CASPT2 calculations using the same active space and 6
states were also performed to compare with FOMO-CASCI. An IPEA shift of 0.25 a.u. was used for the MS-CASPT2 calculations. Initially, the best choice for the optimal broadening parameter (β) was determined by examining how the excitation energies change with β (SM:
Figure 1). The optimum value was selected to be 0.2 for the four singlet excited states on which the dynamics occur. More details on selecting these parameters can be found in SM (Section I). The energies and oscillator strengths for the excited states in the equilibrium geometry of the ground state (S 0 -NP), optimized with SA6-CASSCF(14,10)/cc-pVDZ, were calculated using FOMO-CASCI (β = 0.2) with various basis sets and active spaces. These results are shown in SM (Section I). The cc-pVDZ basis set along with an active space of 12 electrons in 10 orbitals showed the best results and was chosen for the remaining calculations.
In order to further validate the efficacy of FOMO-CASCI in nonadiabatic dynamics, S 0 -S 1 hopping geometries from the singlet trajectories were used as starting geometries to find the unique minimum energy conical intersection geometries (MECIs) using SA6-CASSCF(14,10)/cc-pVDZ. Reaction paths were then produced from the S 0 -NP geometry to each MECI using linear interpolation of internal coordinates (LIIC) using both CASSCF and FOMO-CASCI. The results indicated that the pathways are qualitatively similar, including the character of the first excited singlet state (S 1 ) at the MECI and the energy barriers to reach them (see Section III E and Figure 7 in SM).
All geometry optimizations and single point energy calculations were performed with OpenMolcas (v24.1) 27 and all FOMO-CASCI calculations were performed with TeraChem (v1.9). 22
The ground state o-NP (S 0 -NP) was optimized using SA6-CASSCF(14,10)/cc-pVDZ. The vibrational normal modes of S 0 -NP were sampled using a Wigner distribution 28,29 to capture the geometrical phase space of the ground state vibrational wavefunction assuming a harmonic oscillator along each normal mode. 1450 initial conditions were obtained to generate the absorption spectrum. Excited state energies and oscillator strengths for these geometries were calculated using FOMO-CASCI with a broadening parameter of 0.2 and an active space of 12 electrons in 10 orbitals with the cc-pVDZ basis set (FOMO-CASCI(0.2)(12,10)/cc-pVDZ). The spectrum was plotted by convoluting the oscillator strengths with Gaussian functions with a full width half maximum of 0.25. A large FWHM was chosen to smooth out the spectrum to extract the two most prominent peaks. For comparison, the absorption spectrum was also calculated at the CASSCF and MS-CASPT2 level using 200 geometries (SM: Section II). The experimental spectrum is also shown in SM (Figure 3).
In order to simulate pump-probe experiments, the initial conditions for the dynamics are obtained using geometries corresponding to the excitation window used in experimental work. The energy of the pump pulse used in our previous work 6 was 4.9 eV which is close to the absorption maximum of the second peak in the experimental absorption spectrum of o-NP (5.0 eV). Therefore, initial conditions for the dynamics were chosen from the second peak of the theoretical spectrum, within the energy window of 6.2 ± 0.1 eV (6.2 eV being the maximum at the FOMO-CASCI level), and a cut off for the oscillator strength of 0.06. Of the initial 1450 geometries, 380 were found to have excitation energies within this selection window and an oscillator strength larger than 0.06. Due to limitations in computational resources, only 165 out of the 380 initial conditions were randomly selected and run: 50 exciting to the third excited singlet state (S 3 ) and 115 to the fourth excited singlet state (S 4 ). Although S 4 is the bright state vertically, distortions along the normal modes can increase the oscillator strength for S 3 as well, either by intensity borrowing due to normal modes that break symmetry, or when the ordering of nπ * and ππ * states changes. So the absorption peak has contributions from both adiabatic states. An alternative approach would be to run all trajectories within the energetic window and then weight by the oscillator strengths, but this approach would require many more trajectories and would most likely not change the results.
To compare and contrast the dynamics of o-NP initiated by excitation to the second absorption peak with previous literature, which had focused on exciting to the first absorption band, we also run dynamics starting on the first peak. For this, initial conditions were selected from the first peak of the theoretical spectrum within the 5.2 ± 0.1 eV window with an oscillator strength cutoff of 0.13. This selection led to 51 geometries, with 18 of them populating the S 1 state and 33 the S 2 state.
The excited state dynamics were studied using Trajectory Surface Hopping (TSH), 20 a semi-classical method where nuclear dynamics are treated classically with the velocity-Verlet algorithm, while electronic states and nonadiabatic couplings are handled quantum mechanically. The Fewest Switches Surface Hopping (FSSH) method, 30 which incorporates nonadiabatic couplings (NACs), was used to determine the probabilistic hop between states.
Following a hop, the momentum was rescaled along the derivative coupling vector to conserve total energy, while frustrated hops left the momentum unchanged. Additionally, the decoherence correction by Persico and Granucci 31 with a decoherence factor of 0.1 Hartree was applied to account for the loss of coherence between quantum states. Two sets of trajectories were run for 500 fs using the described TSH method and the initial conditions described above to simulate excitation to the two peaks. When using complete active space methods, sudden changes in the ordering of orbitals in the active space tend to cause large momentary deviations in the total energy of the system. Therefore, to monitor energy conservation, NewtonX calculates the deviation of the total energy of the system compared to the initial total energy and that of the previous timestep. If either deviation is more than 0.5 eV, the trajectory is terminated. The trajectories that terminated before the 500 fs mark were restarted using a time step of 0.25 fs. Despite this adjustment, some trajectories continued to fail with approximately 20% of the trajectories initiated on S 1 /S 2 and 40% of trajectories initiated on S 3 /S 4 remaining active at the 500 fs mark. In the subsequent analysis of state populations over time, these trajectories were included up to the point of failure. In order to examine the effect of including these trajectories in our analysis vs not including them at all we show comparisons in SM, Section IIIC.
The TSH simulations were performed using Newton X (v2.4) 32 and TeraChem (v1.9) 22 using an in-house interface between the two software packages.
Although our simulations do not include the ISC directly in the dynamics, we calculated the spin-orbit coupling (SOC) and used it to estimate the ISC rate and probabilities. This estimate was calculated for trajectories that were on the S 1 surface after relaxation from S 3 /S 4 states, since the initial relaxation is very fast. We only estimated ISC for trajectories that were initiated on the second absorption peak, although it is likely that similar SOC will be found in the trajectories on S 1 initiated on the first absorption peak. Trajectories were then initiated on the triplet surfaces based on the value of these rates.
SOC calculations between S 1 and T 1 -T 3 were performed at 50 fs intervals for the dura-tion of each trajectory that stayed on S 1 using an average of states CASSCF with 6 singlet and 3 triplet states, SA(6S,3T)-CASSCF(14,10)/cc-pVDZ. The SOC was calculated using CASSCF wavefunctions for the singlet and triplet states as implemented in Molcas. 33,34 In this implementation, the Spin-Orbit Hamiltonian has been approximated by a one-electron effective Hamiltonian, which also avoids the calculation of multicenter integrals. 35 We performed the SOC calculation using CASSCF wavefunctions, since we could not do so using FOMO-CASCI. Nevertheless, we checked that the results are not affected by this. We monitored the character of states, and energy differences between S 1 and triplet states for all the geometries that we calculated ISC rates. The results were consistent between the two methods on most cases. In few cases that the character of S 1 was different at the CASSCF level compared to FOMO-CASCI we did not use these values any further. Table V in SM reports the information relevant to ISC.
To evaluate the probability of ISC (p ikt ) from singlet state i to a triplet state k, at each sampling point t, we employ the Landau-Zener formula, [36][37][38][39][40][41]
with
∆E it-kt is the energy gap between the S 1 (i = S 1 ) and the respective triplet state, and
is the spin-orbit coupling to that state, both of which were calculated with CASSCF here. dt is the time step that we used to calculate the SOC, which was 50 fs to avoid running too many CASSCF calculations. We tested the effect of this timestep by comparing the results to dt = 5 fs for one trajectory using energies at every 5 fs. The comparison is shown in SM and demonstrates that it is not a too severe approximation. Using this formula, we calculated the cumulative probabilities for ISC between S 1 and triplet states T 1 , T 2 , T 3 , along the duration of all trajectories that remained on S 1 . The total probability of intersystem crossing (P T ) occurring at least once to any of the three tested triplet states (k = 1, 2, 3) from S 1 (i = 1) was computed across all sampling points in the trajectory (t = 1 to n) as follows:
In addition, we used this expression to derive an effective rate at each time step (of 50 fs) that we used to decide where to launch the triplet trajectories. The effective rate we used to decide whether to launch trajectories on triplet surfaces is,
Other more formal expressions of rates for intersystem crossing also depend on the square of SOC. [42][43][44] Fedorov et al. 45 have used an effective coupling SOC/∆E which is dimensionless, thus related to probability, while the property we chose has units of rate (s -1 ).
For each singlet trajectory staying on S 1 for more than 50 fs, the point of maximum effective rate γ ikt was selected as the starting geometry for each triplet trajectory, while disregarding any trajectories whose point of maximum coupling fell below 5000 cm -2 /eV.
Converting eV to cm -1 and using the speed of light, the cutoff is equal to 1.85 × 10 10 s -1 .
This cutoff was based on some reasonable values for the SOC and ∆E that we would expect ISC to start to occur (for example, SOC = 40 cm -1 and ∆E = 0.3 eV gives a value close to the cutoff), and it results in a p ikt = 0.6%. The velocities of the atoms at the selected points were used as the starting velocities on the triplet surfaces.
In the S 3 /S 4 dynamics, out of the 77 trajectories that remained on S 1 , 47 showed ISC rates high enough to cross to the triplet manifold. Accordingly, 13 trajectories were started on the second excited triplet state (T 2 ) and 34 on the first excited triplet state (T 1 ). Further information on the starting points of the triplet trajectories can be found under Section V of SM. All trajectories on triplet surfaces were carried out using FOMO-CASCI(0.2)(12,10)/cc-pVDZ with a timestep of 0.5 fs. Trajectories that failed before 200 fs were rerun with a timestep of 0.25 fs. The trajectories were allowed to run for 1-2 ps, in order to see if dissociation occurs, although many failed before that time. The same parameters and algorithms used for the singlet states were applied to the triplet states using NewtonX (v2.4) and TeraChem (v1.9).
The UV absorption spectrum of o-NP consists of excitations to the first four excited states of the molecule. Two of these states are bright ππ * states, while the other two are nπ * states with very small absorbance. The orbitals involved in these excitations are shown in Figure 1. Table I shows the energies of these states calculated at different levels of theory in this work compared with previous experimental studies. [12][13][14] MS-CASPT2 predicts S 2 and S 4 to be the bright states contributing to the two absorption bands. The vertical excitation energies however are blue shifted by about 0.7 eV. This error is very sensitive to the choice of active space and states included in the multi-state description. The gap between the two bright states is about 1 eV experimentally, and it is reproduced well by MS-CASPT2.
CASSCF has a much larger error for the two bright states, 1.5 eV for the first bright state and 1.8 eV for the second. The gap between them is also larger than experiment, 1.4 eV.
Furthermore, the first bright state is predicted to be S
3 . FOMO-CASCI also has a large error compared to experiment; the error, however, is smaller than CASSCF, 1.3-1.4 eV blueshifted compared to experiment. It also reproduces the gap between the two bright states, and predicts the first bright state to be S 2 , in agreement with MS-CASPT2. So, even though FOMO-CASCI is an approximation of CASSCF it seems to give better results than CASSCF for the ππ * states. There are no experimental values for the nπ * states, but comparison with MS-CASPT2 shows that the FOMO-CASCI errors for these states are similar to the errors for ππ * states. On the other hand, the CASSCF energies for the nπ * states are much closer to the MS-CASPT2, since dynamical correlation is not as important for these states as it is for the ππ * states. The overall result of these comparisons is that FOMO-CASCI should give a more balanced description between the nπ * and ππ * states and possibly be closer to the MS-CASPT2 dynamics compared to CASSCF.
Figure 2 shows the absorption spectrum calculated from the Wigner sampled geometries using FOMO-CASCI(0.2)(12,10)/cc-pVDZ. The absorption spectrum was also calculated at the CASSCF and MS-CASPT2 level, and these spectra are shown in SM (Figure 4). Two peaks were observed, just as in the experimental spectrum, both with a blue shift of about 1.2 eV. The separation between the two peaks is about 1 eV in both theory and experiment.
The first peak corresponds to the first ππ * state, which is associated with S 1 or S 2 in the calculated spectrum, while the second peak is associated with the second ππ * state, which for these geometries can be S 3 or S 4 . Even though at vertical excitation the bright states are S 2 and S 4 , when we include the zero point vibrational motion through the Wigner distribution the character of the states can switch. This is because either the states switch character along displacements, as illustrated in Figure 1, or vibrations along some modes break the planar symmetry leading to mixing and intensity borrowing. That is what we see in the spectrum where all states have some contribution to the peaks.
The theoretical maxima of the absorption spectra are included in Table I. These values are a more direct way to compare to the experimental absorption maxima, rather than the vertical excitations that we used in the previous section. It is very interesting to note that the theoretical absorption maxima are not exactly the same as the vertical excitations. The absorption maxima are red-shifted compared to vertical excitation energies by as much as 0.2 eV.
Initial conditions were selected from either the first or the second peak of the theoretical spectrum. These two separate excitations correspond to different experimental pump wavelengths. A main focus here is to examine the dependence of the dynamics on the excitation wavelength.
TABLE I. Energies (E) and oscillator strengths (f) of the first four singlet excited state with respect to the ground state (S 0 ) for the S 0 -NP geometry using various methods compared with gas phase experimental results (EXP). 12,13 The peak maxima (M), for the theoretical spectra calculated using 200 Wigner geometries with each method are also reported. FOMO-CASCI: FOMO-CASCI(0.2)(12,10)/cc-pVDZ, CASSCF: SA6-CASSCF(14,10)/cc-pVDZ, CASPT2: SA6-CASSCF(14,10)/MS-CASPT2/cc-pVDZ . EXP 12,13 CASPT2 CASSCF FOMO-CASCI State E (eV) E (eV) f M (eV) E (eV) f M (eV) E (eV) f M (eV) S 1 -4.32 0.00 -4.41 0.00 -5.09 0.00 -S 2 3.9 4.62 0.14 4.44 5.06 0.00 -5.21 0.10 5.2 S 3 -4.88 0.00 -5.43 0.16 5.36 5.98 0.00 -S 4 5.0 5.68 0.35 5.58 6.81 0.23 6.76 6.38 0.25 6.2
The semi-classical simulation of o-nitrophenol's response to excitation by a pump laser pulse was conducted using TSH. In these simulations, trajectories initiated from S 1 and S 2 correspond to a pump pulse at 3.9 eV, which aligns with the first UV absorption peak (Table I) and reflects the conditions used in recent UED experimental work with a 330 nm (3.8 eV) pump pulse. 7 In contrast, trajectories starting from S 3 and S 4 simulate excitation by a 5.0 eV pump pulse, resembling the conditions of a recent TRPES experiment conducted with a 255 nm (4.9 eV) pump pulse. 6 Nonadiabatic dynamics under these two excitation conditions were analyzed to understand how the initial electronic state influences relaxation kinetics and product distributions.
For the high-energy scenario, 165 trajectories were initiated from S 3 and S 4 with an initial population distribution of approximately 70% on S 4 and 30% on S 3 . For the low-energy case, 51 trajectories were initiated on S 1 and S 2 , with 65% on S 2 and 35% on S 1 . The populations for all states as a function of time are shown in Figure 3. In both cases, state populations were interpreted using a kinetic model for a sequential process (Equations 7 -10). 46 The populations of the rapidly decaying excited states were grouped as the "reactant", N 2 , that feeds into the intermediate S 1 state, with population N 1 , which subsequently decays to the ground state S 0 , population N 0 , which is treated as the "product." In Equations 8 -10, A 1 , A 2 , A 3 are pre-exponential constants, t 1 and t 2 are time constants, C is a constant and t is time. In the high-energy dynamics, we combined the populations of S 2 , S 3 , and S 4 as N 2 , while for the lower energy dynamics the population of S 2 corresponds to N 2 .
The time constants derived from fitting the variation of the populations of the electronic states with time to the above equations are shown in Table II, while all constants are shown 0 100 200 300 400 500 Time (fs) 0.0 0.2 0.4 0.6 0.8 1.0 Population S 3 /S 4 S 0 S 1 S 2 S 3 S 4 S 2 +S 3 +S 4 0 100 200 300 400 500 Time (fs) 0.0 0.2 0.4 0.6 0.8 1.0 Population S 1 /S 2 S 0 S 1 S 2 FIG. 3. Populations of states with fits for the trajectories started on S 3 /S 4 (left) and S 1 /S 2 (right).
in SM. In the high-energy dynamics, the combined decay of S 2 , S 3 , and S 4 into S 1 occurred with a time constant of 51 fs. In contrast, the low-energy dynamics showed an even faster decay of S 2 into S 1 , characterized by a time constant of 20 fs. The subsequent decay from S 1 to S 0 is slower in both cases, with a time constant 330-430 fs. Given the relatively long lifetime of S 1 in both cases, there is the possibility of competing processes, such as intersystem crossing, provided there is sufficient spin-orbit coupling and a small energy gap between S 1 and nearby triplet states. It should also be noted that there is greater uncertainty for the timescale for S 1 decay and S 0 rise compared to the very fast decay of the initial states. This is because over time many trajectories fail and are excluded from the population thus decreasing the statistics.
Although the decay kinetics from S 1 to S 0 were consistent across both scenarios, the resulting photoproduct distributions varied significantly, as discussed next.
TABLE II. Time constants for population on singlet surfaces for trajectories started from S 3 /S 4 and S 1 /S 2 . Time Constant S 3 /S 4 S 1 /S 2 t 1 (fs) 50.9 20.0 t 2 (fs) 422.9 330.0 D. Excited state intramolecular proton transfer (ESIPT) 0 100 200 300 400 500 Time (fs) 1 2 3 4 5 NO-H Distance (Å) S 3 /S 4 H transferred H not transferred 0 100 200 300 400 500 Time (fs) 1 2 3 4 5 NO-H Distance (Å) S 1 /S 2 H transferred H not transferred FIG. 4. NO-H distance as a function of time for all singlet trajectories for dynamics on S 3 /S 4 (left) and S 1 /S 2 states (right). Based on this distance, 27% of trajectories showed ESIPT on S 3 /S 4 dynamics, and 70% of trajectories on S 1 /S 2 . Dots represent the points where hops to the ground state occur.
ESIPT is defined by the transfer of proton from the OH group to the NO 2 group leading to aci-nitrophenol, while the electronic wave packet is in an excited state. Figure 4 illustrates the NO-H bond length, which monitors ESIPT, for singlet trajectories initiated in each of the two absorption peaks. ESIPT occurs at geometries that remain close to planar, as can be seen in SM, Figure 5.
For trajectories starting on S 3 /S 4 , the length of the NO-H bond is initially approximately 2 Å, which corresponds to the enol structure before proton transfer. Figure 4 displays the bifurcation of these trajectories, revealing that about 27% undergo ESIPT. In contrast, for the trajectories initiated from S 1 /S 2 , ESIPT is observed in 70% of the trajectories. These differences suggest that the initial electronic state preparation has a pronounced effect on the branching pathways on the S 1 surface.
ESIPT occurs on S 1 prior to reaching the S 0 -S 1 conical intersection for all trajectories that exhibit ESIPT. Figure 5 shows the time it takes for ESIPT to occur, measured as the average time it took the H atom to travel 1.5 Å away from the O atom of the OH group. ESIPT occurs mostly in less than 100 fs for the high energy trajectories, and even faster, less than 50 fs, for the low energy trajectories. The figure also makes clear that ESIPT takes place on the excited state, preceding nonadiabatic transitions to the ground state.
E. The S 1 -S 0 conical intersections Nonadiabatic transitions to the ground state are facilitated by S 1 -S 0 conical intersections.
Using hopping geometries from trajectories initiated on S 3 and S 4 as initial guesses, we identified five unique minimum energy conical intersection (MECI) structures, shown in The lowest energy MECI is CI1, which involves ESIPT and the character of S 1 is ππ * . CI1 and CI2 are similar, since they both involve ESIPT and the S 1 state is ππ * . Their main difference is the rotation of the NO 2 H group. CI2 should be reached first, since the H is still between the O and NO 2 groups, while in CI1 it has rotated by almost 180°. So, even though CI1 has lower energy, CI2 is probably easier accessible. The other three MECIs do not involve ESIPT and the S 1 state has nπ * character. CI3 and CI4 differ by the NO 2 rotation again, while CI4 involves a different rotation of the OH group. In that sense, they are just different conformers of each other.
FOMO-CASCI predicts the energies of MECI that involve the ππ * character (CI1 and CI2) with similar accuracy to that of CASSCF, but those that involve the nπ * character are destabilized considerably compared to CASSCF. We have also seen a difference in nπ * energies for the vertical excitations, as seen in Table I. The decreased accuracy of FOMO-CASCI in predicting nπ * states compared to CASSCF has been discussed before. 21 Quite interestingly, FOMO-CASCI agrees better with CASPT2, having smaller differences for both ππ * and nπ * states. This is consistent with what we found for the vertical excitations.
TABLE III. Character and energy of the S 1 state for significant geometries and all unique MECIs. S 1 -NP: The minimum energy o-nitrophenol conformation on S 1 . S 1 -aci-NP: The minimum energy ESIPT o-nitrophenol conformation on S 1 . The methods listed are CASSCF: SA6-CASSCF(14,10)/cc-pVDZ, CASPT2: SA6-CASSCF(14,10)/cc-pVDZ/MS-CASPT2, FOMO-CASCI: FOMO-CASCI(0.2)(12,10)/cc-pVDZ. CASSCF CASPT2 FOMO-CASCI Geometry Character E (eV) E (eV) E (eV) ESIPT S 0 -NP n → π * 4.37 4.32 5.09 No S 1 -NP n → π * 3.13 3.41 4.00 No S 1 -aci-NP π → π * 3.53 3.39 3.84 Yes CI1 π → π * 3.01 3.58 3.87 Yes CI2 π → π * 3.57 3.55 3.87 Yes CI3 n → π * 3.62 4.04 4.76 No CI4 n → π * 3.36 4.30 4.94 No CI5 n → π * 3.55 4.31 4.74 No
MECI similar to CI1 and CI2 that involve ESIPT have been previously reported in the literature [7][8][9]16 . The pathways connecting the FC region to these MECI have also been calculated and have been shown to be barrierless.
7,8 However, the other MECIs, without ESIPT, have not been reported before. The interpolated potential energy surfaces connecting the Franck-Condon region to each MECI are shown in SM: Section IV. The interpolated reaction pathways calculated using FOMO-CASCI and CASSCF are qualitatively similar, suggesting similar dynamics predicted by the two methods. The MECIs without ESIPT (CI3, CI4, CI5) are almost barrierless. The MECI leading to CI2 was shown to be barrierless if an intermediate planar S 1 minimum was used. 7
In order to understand better which MECI are mostly responsible for the hops observed in TSH, we chose two important coordinates to depict the geometries for MECI and hops.
These coordinates are the NO-H distance (ESIPT coordinate) and the CC-NOH dihedral angle. Figure 7 displays the hopping geometries and MECI geometries as a function of these coordinates for the two sets of trajectories initiated on S 3 /S 4 , as well as S 1 /S 2 . This allows a better insight into which MECI is mostly responsible for the nonadiabatic transitions.
For the former set, which begins at higher energy levels, the majority of transitions to the ground state are concentrated around the CI3 geometry. At this point, no ESIPT occurs, but the NO 2 group is slightly rotated (∼20 degrees) out of the plane. For the low energy excitation, CI2, where ESIPT has occured, is preferred. This agrees with the dynamics that showed that more ESIPT is occurring when starting from the low energy absorption peak.
Franck-Condon excitation to S 3 /S 4 leads to a fast decay to S 1 followed by slower decay to the ground state, as shown in Figure 3. We calculated SOC for the trajectories that remain on S 1 in order to evaluate if they were possible to transition to the triplet manifold via ISC.
In this study, we adopted an approximate method to incorporate the effects of ISC. We calculated the SOC along the S 1 trajectories and the energy gap between S 1 and nearby triplet states. Figure 8 presents a scatter plot of these values for all trajectories using CASSCF values, while Table V in SM has a detailed list. Several trajectories demonstrated notably high SOC values, particularly between S 1 -T 1 . Concurrently, the small gap between these states resulted in high ISC rates. Similar high coupling values were also observed for certain S 1 -T 2 cases. We calculated the effective rate along the S 1 trajectories using Eq. 6, and identified the point with the highest value to launch a triplet trajectory. These points are highlighted in Figure 8 with black circles. An example of the coupling along a trajectory is given in SM (see SM: Figure 8).
As shown in Figure 8, ISC is predicted mostly from S 1 to T 1 for SOC values of 60-70 cm -1 . The S 1 state has nπ * character, while T 1 is ππ * , and this occurs in both ESIPT and non-ESIPT geometries. In some cases, T 2 is almost degenerate to S 1 and in this case ISC is predicted to T 2 , even though the SOC is somewhat smaller, 25-60 cm -1 . The lower values of SOC are when S 1 is a ππ * state coupled to T 2 with nπ * character. In general, ISC is predicted from S 1 regardless of whether ESIPT has occurred, although the SOC is smaller for structures with ESIPT. As a reference, we have calculated SOC at the two S 1 minima.
At S 1 -NP the SOC between S 1 and T 1 is 64 cm -1 , while at S 1 -aci-NP the highest SOC is between S 1 and T 2 at 24 cm -1 (see Table VI in SM).
Cumulative probabilities for ISC from S 1 to any triplet state have also been calculated using the Landau Zener formula in Eqs. 3,5. These results are tabulated in SM, Table V, and they indicate that the probabilities for ISC are on average 3-4% per 100 fs, or 30-40% within 1 ps.
A total of 47 trajectories were initiated on the triplet manifold, with 13 starting on T 2 and 34 on T 1 . 16 of the 47 trajectories already had ESIPT at the starting geometry. Some population from T 2 transferred to T 1 , although most did not. The populations are shown in the SM (see SM: Figure 10). Interestingly, within the first 100 fs there is some ESIPT back and forth from the OH group to the NO 2 group, but the majority of trajectories retain their structure. As a result, these trajectories suggest that the majority of ESIPT occurs on S 1 before ISC, while minimum ESIPT occurs after ISC.
Two trajectories that were run on the triplet manifold stood out from the rest due to the dissociation of previously predicted species, NO and OH. Figure 9, shows representative steps for the mechanisms followed by those trajectories.
For the trajectory that led to OH dissociation, ESIPT on S 1 occurred before ISC. The for all the geometries that were selected as the starting points of the triplet trajectories. All calculated at the SA(6S,3T)-CASSCF(14,10)/cc-pVDZ level. Circled in black are the states on which each triplet trajectory was started. largest calculated coupling was observed between S 1 and T 1 (57 cm -1 ) at 193.5 fs along the singlet dynamics with an energy gap of 0.3 eV and an effective rate of 2.1 × 10 10 s -1 . Barrierless dissociation of OH was observed after several hundred fs on the triplet surface.
It is important to note that OH dissociation requires ESIPT prior to ISC. At the point of simulated ISC, the character of S 1 was nπ * and the character of T 1 was ππ * .
For the trajectory that led to NO dissociation, following internal conversion to S 1 there is no ESIPT. The largest coupling was observed to T 2 (27.3 cm -1 ) at 291.5 fs with ∆E = 0.04 eV and an effective rate of 5.5 × 10 10 s -1 . Once on T 2 , it rapidly relaxes to T 1 and the NO 2 group cyclizes by attaching one of its oxygen atoms to the same carbon atom on the benzene ring that the NO 2 group is attached to, forming an intermediate, which we previously called the "spiro intermediate". 6,18 Then, NO dissociates from the spiro bicyclic intermediate in about 650 fs on the triplet surface. The energies and motion along the two trajectories are shown in Figure 9.
These trajectories reveal that dissociation occurs primary on the triplet surfaces. This is in agreement with previous hypotheses. 2,6 On the other hand, we only saw two events in our simulations, which indicates that dissociation is probably slower than our simulation times.
In our previous work, the NO dissociation was observed very fast, much faster than what we see here, while we did not observe any OH dissociation. 6 The cause of the latter is that OH dissociation was not visible with our probe because of the high ionization potential of the fragments. For the NO dissociation, the experiment and theory agree qualitatively, but not quantitatively.
This study provides an excellent model system to examine the effect of excitation wavelength on the subsequent dynamics. Initial excitation to S 3 /S 4 leads to very fast relaxation to S 1 (with a time constant of about 50 fs). Relaxation from S 1 to S 0 is much slower, with a lifetime of about 400 fs in both sets of dynamics, i.e. regardless of the initial excitation. This would suggest that the dynamics are independent of excitation, but a closer look at the hopping geometries as shown in Figure 7 shows that the radiationless decay occurs through different pathways. When starting on S 3 /S 4 most of the nonadiabatic transitions occur before ESIPT through conical intersections resembling CI3. On the other hand, initial excitation to S 1 /S 2 leads to nonadiabatic transitions to the ground state after ESIPT has occurred, through a conical intersections resembling CI2. Similar dynamics on S 1 have been found in previous studies when starting from the low energy peak. 7,11 Furthermore, the rate of ESIPT itself was found to be highly dependent on the excitation wavelength.
The difference observed on the S 1 relaxation depending on the initial excitation is quite intriguing. It can, however, be explained by considering the motion of the wavepacket focusing particularly on the coordinate that facilitates ESIPT. Previous dynamical studies showed that ESIPT is best facilitated by the two O approaching each other. 7 In order to examine this motion we have plotted the O-O distance for the first 50 fs of the dynamics for both sets of trajectories. Figure 10 shows the average O-O distance and standard deviation for all trajectories. Dynamics starting on S 1 show a coherent motion of the wavepacket leading to a contraction of that distance followed by ESIPT. On the other hand, when starting on the higher states, there are several nonadiabatic crossings occurring before reaching S 1 , and the energy may be more dissipated to other degrees of freedom. So, the O-O distance does not contract as much and there is less ESIPT. Interestingly the standard deviation is not very different between the two cases, indicating that a different direction may be taken for the S 3 /S 4 dynamics. The percentage of ESIPT is strongly correlated with the O-O distance.
As a result different conical intersections are reached for decay to the ground state.
Studies centered around Franck-Condon excitation close to 3.9 eV (lower energy experimental peak in Table I) [7][8][9]11,16,17 , suggest that the most prominent relaxation pathway is ESIPT followed by the rotation of the newly formed HONO group out of the plane of the ring that facilitates either the relaxation to the ground state (S 0 ) or ISC to T 1 . The combined UED and AIMS study by Nunes et al. 7 correlate the time dependent changes in geometrical parameters seen in UED to the evolution on nonadiabatic surfaces from AIMS to propose decay to S 0 through a CI that resembles CI2 in Figure 6. This is also corroborated in the work by Ernst et al. 8 and Nitta et al. 9 through TRPES experiments coupled with potential energy scans along important reaction coordinates. In the current study, ESIPT was observed at sub 50 fs, which agrees with all of the studies cited above. The decay from S 1 to S 0 in our work, 400 fs, agrees with the UED experiments, but it is slower than the TRPES measurements, that have a decay time of 130 fs, and the AIMS dynamics that predict a decay time of 250 fs.
The decay times in our work come from the dynamics on singlet surfaces only, so the decay time is overestimated, since it only includes the radiationless transition to the ground state, but it does not include ISC. Even though we include ISC estimates, we do not predict the time it takes to reach the triplet states, since ISC is not included directly in the surface hopping algorithm. Nevertheless, the SOC couplings and ∆E we calculate clearly suggest that ISC should be occurring, and we have estimated that there is a 30-40% probability that ISC occurs within 1 ps. The AIMS study did not observe significant SOC between S 1 and any nearby states in the triplet manifold, 7 whereas Xu et al. predicted a very fast ISC (occurring in less than 50 fs). 17 Our results seem to be somewhere in between those two extremes.
The occurrence of ISC is further supported by the dissociation products. As discussed in our previous work, 6 but also shown in the current trajectories, NO dissociation can only occur on the triplet surfaces. In our previous study we calculated the barriers to the dissociation of NO on the neutral first excited singlet and triplet states as well as the ground ionic state.
The barrier to dissociation of NO through the "spiro intermediate" was the smallest on the first excited triplet state, indicating that dissociation of NO was most likely to occur in the first excited triplet state and unlikely in the ionic surfaces. It is reinforced in this study as well that detection of the NO channel is a clear indication that ISC is occurring.
ISC from S 1 to triplet states is predicted here for both isomers of o-NP on the S 1 surface, regardless of whether ESIPT has occured or not. Given that less ESIPT occurs when starting on S 3 /S 4 , we should expect that the NO fragmentation is more prominent when exciting to the higher energies, while OH dissociation should be more important for the lower energy excitations.
In our previous work where o-NP was pumped at 4.9 eV, we detected a time-dependent photo-ion signal with a mass of 109 amu, indicative of the NO elimination reaction. While the triplet trajectory demonstrated NO elimination, it does not explain the statistical correlation between the parent o-nitrophenol ion and the hydroxyphenoxy ion, the elimination product. Instead, it offers evidence for a qualitatively feasible mechanistic pathway that is both energetically and kinetically sound. The formation of NO must occur within 250 fs in the experiment, as the TRPES signal decays entirely by this time. Conversely, the NO elimination identified in this research occurs roughly at times closer to 1 ps. However, we only have one trajectory in our simulation leading to this product, so one would need a much bigger statistical sample to explain the experiment. Furthermore, our inclusion of ISC is rather approximate, so this will also introduce an error.
The work by Cheng et al. 2 and Wei et al. 3 presents experimental evidence for the elimination of OH radicals. They measured the density of OH radicals following the excitation of o-nitrophenol at its major absorption peaks. Their observations showed nascent OH formation. Through density functional theory calculations for the Franck-Condon geometry and key intermediates involved in the OH elimination mechanisms, they predicted that OH elimination is most likely to occur on the triplet surface after ESIPT and an out-of-plane rotation of the OH group on aci-nitrophenol, which aligns with our findings. However, the time resolution in this study was insufficient to resolve ultrafast processes occurring on femtosecond timescales. Thus, definitive experimental information regarding the timescale of OH elimination from o-nitrophenol remains unavailable. Other studies have suggested that OH is produced on the hot ground state, 9 but our study supports its formation on the triplet surfaces. In this work we examined the excited state dynamics of o-NP when excited either on the first or second absorption peak. The results are summarized in Figure 11. The results indicate that excitation to the higher energy peak (which corresponds to electronic states S 3 and S 4 ) leads to fast decay to S 1 (in about 50 fs) followed by a much slower decay from S 1 to S 0 (about 400 fs). The timescales for the decay to S 0 is similar regardless of whether the first or second peak was excited, indicating that pump probe experiments can give similar decay rates.
Nevertheless, the pathways for decay are different depending on the initial excitation.
Both dynamics showed a bifurcation on S 1 with some trajectories undergoing ESIPT while others did not. The bifurcation was different however depending on the initial excitation.
A. Stolow, A. E. Bragg, and D. M. Neumark, "Femtosecond time-resolved photoelectron spectroscopy," Chemical Reviews 104, 1719-1758 (2004).