Excited state dynamics of cis,cis-1,3-cyclooctadiene: UV pump VUV probe time-resolved photoelectron spectroscopy

Understanding the dynamics of photoexcited molecules is of great importance to many areas of chemistry, physics, and biology, including light harvesting,^1–4 energy conversion,⁵ and the photoprotection of DNA.^6–8 The coupled electron nuclear dynamics initiated by photoabsorption can lead to internal conversion,^9–13 isomerization,^14–18 and dissociation.^19–23 While these dynamics can be quite complicated and difficult to follow, particularly for larger molecules with many degrees of freedom, systematic studies of similar molecules can be very helpful in elucidating essential features in the dynamics and a deeper understanding of the connection between the structure and dynamics.

Here, we extend earlier work on simple organic molecules, in an effort to understand isomerization and internal conversion in systems with a few C–C double bonds. The most basic unit to examine photoisomerization is the smallest molecule with a single double bond: ethylene, which has been studied extensively.^24–29 Beyond ethylene, many biological chromophores, however, such as retinal, include a more extended linear or circular conjugated system with several double bonds.³⁰ 

Several measurements and calculations have extended our understanding of the isomerization dynamics in small conjugated organic molecules. However, these studies have focused on relatively rigid small molecules. Two such systems that have drawn considerable interest in past decades are trans-1,3-butadiene (BD, C₄H₆)^31–36 and 1,3-cyclohexadiene (CHD, C₆H₈).^18,37–42 We aim at extending our understanding to larger, more flexible systems, where we might expect qualitative differences in the internal conversion and isomerization dynamics. A natural choice that extends the previous studies is cis,cis-1,3-cyclooctadiene (cc-COD, C₈H₁₂), which is similar to CHD, but larger and more flexible. In both CHD and cc-COD, there are two alternating double bonds, but the remaining cyclic system imposes different constraints on the dynamics. In CHD, the dynamics are controlled by conical intersections (CIs) with the ground state that can lead to either relaxation or photochemistry to produce hexatriene, which has been studied extensively with both time-resolved spectroscopies and diffraction measurements.^42–44 In cc-COD on the other hand, there can be cis–trans isomerization, as well as photochemistry to other products.^45–47 

Compared with CHD, the isomerization dynamics of cc-COD have not been well studied by either spectroscopic or structural probes. Fuß et al. performed UV pump–IR probe ionization measurements on cc-COD and established basic time scales for the dynamics and possible photochemical products, but did not provide a detailed picture of the dynamics.⁴⁵ 

In this work, we present time-resolved photoelectron spectroscopy (TRPES)^48–53 measurements of cc-COD. We interpret the experimental results together with electronic structure calculations for the TRPES signal and trajectory surface hopping (TSH)^54,55 calculations reported separately.⁵⁶ Our measurements show rapid internal conversion in less than 100 fs. The photoelectron spectra as functions of pump probe delay can be interpreted in terms of ionization to multiple states of the molecular cation that behave differently along the reaction coordinate associated with the neutral dynamics.

Our experiments make use of light from an amplified Ti:sapphire laser system, which produces 30 fs, 1 mJ pulses at a central wavelength around ∼780 nm and a repetition rate of 1 kHz. UV (∼261 nm, ℏω = 4.75 eV) and vacuum ultraviolet (VUV) (∼156 nm, ℏω = 7.92 eV) pulses are generated by the third harmonic generation in crystals and non-collinear-four-wave-mixing (3ω + 3ω − ω = 5ω) in argon, respectively.^57–59 The pump–probe scheme has been described in detail elsewhere.^47,60,61 Figure 1 shows our UV pump pulse spectrum together with the calculated absorption spectrum of cc-COD.

The UV and VUV pulses are combined together with a molecular beam in the interaction region of a coincidence velocity map imaging (VMI) apparatus,⁶² and we measure the momentum resolved photoelectron yield as a function of pump–probe delay. The pump–probe time-zero is determined by TRPES measurements in ethylene, which has a rapid decay time.^12,24,63 By fitting the ethylene signal with a Gaussian and an exponential function, we acquired the time-zero and instrument response function (IRF) of the UV–VUV cross correlation, which is ∼80 fs. The details of the fitting results are included in  Appendix A.

The VMI spectrometer employs a continuous molecular beam, an electrostatic lens, a micro-channel plate (MCP), and a phosphor based position-sensitive detector. The electrostatic lens consists of a standard repeller–extractor–ground electrode lens sitting inside a μ-metal sheet cylinder for magnetic shielding. The molecular nozzle with a diameter of 30 μm is set parallel to the surface of the repeller and extractor plates about 3 mm away from their edges. The pump and probe beams propagate parallel through the plates perpendicular to the molecular nozzle. At the end of the VMI spectrometer, there is a position-sensitive 2D detector, consisting of two 40 mm active diameter MCPs in the Chevron configuration together with a P47 phosphor screen.

The chamber base pressure is ∼10⁻⁷ Torr. The running pressure of the reaction chamber is kept at ∼10⁻⁶ Torr. cc-COD is purchased from Sigma-Aldrich (98%) and used without further purification. A multi-functioning sample manifold is integrated into the VMI chamber, which allows for both liquid and gas samples. The manifold contains a liquid sample cell, a delivery line, and a molecular nozzle inside the reaction chamber. Several valves are used for pumping and sample pressure regulation. All of the components are made of stainless steel and copper in order to prevent leaking and contamination. cc-COD, which is a liquid under room temperature, is loaded into the sample cell, pumped on while frozen (“freeze pumped”) in order to remove air from the sample, and then regulated with a precision leak valve between the sample cell and the nozzle. The pressure behind the nozzle is kept at several Torr, which generates a continuous molecular beam through the nozzle to the VMI chamber.

The photoelectrons are projected to the detector, and the fluorescence signal from the phosphor is recorded by using a fast conventional frame based CMOS camera (BASLER acA2000-340 km).⁶⁴ Both the pump and probe pulses are linearly polarized with the polarization direction perpendicular to the ToF direction. An Abel inversion algorithm is used to recover the full 3D vector momentum distribution of the photoelectrons based on the measured 2D projection with a cylindrical symmetry about the laser polarization axis. The photoelectron spectrum (yield as a function of energy) is then calculated from the 3D momentum distribution.

In order to interpret the experimental measurements, we calculated the TRPES using information from our previously published trajectory surface hopping calculations.⁵⁶ For a more in depth understanding of which states are involved in the TRPES, and how they influence the dynamics, we have also performed calculations along linear interpolation paths. Details on these calculations are provided here.

We calculated energies and Dyson norms along pathways leading to conical intersections between S₁ and S₀. The structures of conical intersections are obtained from our previous theoretical study on cc-COD.⁵⁶ Using three important conical intersection geometries (as described below), we constructed linear interpolations in internal coordinates (LIICs) from the Franck–Condon (FC) geometry to the CIs using the extended multi-state complete active space with second order perturbation (XMS-CASPT2) method^65–67 and the cc-pVDZ basis set.⁶⁸ For neutral states, a (6,6) active space was used, whilst a (5,6) active space was used for cationic states. The XMS-CASPT2 calculations were performed using the corresponding complete active space self-consistent field (CASSCF)⁶⁹ reference wavefunction with an imaginary shift of 0.2 a.u. and the single-state single-reference (SS-SR)⁷⁰ contraction scheme. The energies of three states for the neutral and five states for the cation were calculated. The XMS-CASPT2 calculations were performed using the Bagel package.^71,72

We calculated the Dyson norms from S₁ and S₂ states to all of the five cationic states along the aforementioned LIICs at the CASSCF/cc-pVDZ level³⁶ using MOLPRO 2015.^73,74 The same active spaces were used in the Dyson calculations as those in the energy calculations. Care was taken to ensure the calculation of correct Dyson norms for each point of the LIICs as the neutral excited states switch character near the FC region at the CASSCF(6,6) level.

Using trajectories from non-adiabatic dynamics simulations, the photoelectron signal S for each trajectory and time step t can be calculated assuming that the transition is vertical, i.e., the nuclear wavefunctions of the initial and final states are identical. The photoelectron spectrum at each delay is then given by

where ΔV_IF = V_F − V_I is the difference of the adiabatic electronic energies of the initial neutral I and final ionic F states and ε_k is the kinetic energy of the ejected electron. k is the momentum of the ejected electron, ℏω is the probe photon energy, and c is the speed of light. D_IF is the photoelectron dipole matrix element, which is given by

Here, $μ$·u is the scalar product of the dipole operator and the unit vector along the laser polarization axis. $ψIFd$ is the Dyson orbital, defined as the overlap between the initial neutral electronic state and the final state of the cation after ejection of the electron (assuming the photoelectron ejection is fast, the state of the cation does not interact with the outgoing electron). In the current work, we ignore the wavefunction of the ejected electron and approximate |D_IF|² with the square of the Dyson norm. This approximation has been used very often in the literature, and several studies have shown that it has a small effect on the photoelectron spectrum.^12,75,76

In order to calculate the photoelectron spectrum, we selected trajectories at the lower edge of the absorption spectrum from our previously reported TSH calculations using CASSCF wavefunctions⁵⁶ (corresponding to the 4.75 eV pump pulse). Time-dependent Dyson norms were calculated for each trajectory with a time step of 10 fs. We used the same level of theory [SA3-CASSCF(4,3)/cc-pVDZ] to calculate the neutral wavefunctions as that we used in the original TSH simulations, and we used SA5-CASSCF(3,3)/cc-pVDZ to calculate the cationic wavefunctions. Note that, at this level of theory, D2 state has the same character as the D3 state in the previously calculated LIICs at the FC region. For calculating the electron kinetic energies (KE), specific shifts were introduced to the KE of electrons from D₀, D₁, and D₂, since CASSCF overestimates the neutral excited state energies and underestimates the cationic state energies near the FC region. Shifts of 3.0 eV, 2.8 eV, and 2.3 eV were introduced to the KEs due to photoelectrons coming from ionization to D₀, D₁, and D₂, respectively. These shifts were based on a comparison of the neutral S₁ energy to the experimental absorption peak and the cationic state energies to those calculated at the XMS-CASPT2(5,6)/cc-pVDZ level. D₃ and D₄ are quite high in energy at the CASSCF(3,3) level and, hence, are not accessible and were neglected for the calculation of the photoelectron spectrum. The energies of the photoelectrons at each delay step for each trajectory were convolved with a Gaussian function having a 0.5 eV width in order to account for the finite energy resolution of our measurements. Along the pump–probe delay axis, the calculated spectra were also convolved with the IRF from the experimental measurements.

In order to check the validity of the CASSCF TSH dynamics on which we base the calculation of the TRPES, we also carried out TSH calculations at the XMS-CASPT2 level. The details of these calculations are given in Appendix C 2.

In Fig. 2, we show both the measured and calculated time-resolved photoelectron spectra for cc-COD. As the measurements contained contributions from both UV and VUV driven dynamics near zero time delay given our finite pulse durations, we first performed a global 2D fitting analysis and subtracted off the portion of the signal due to VUV driven dynamics. For the calculated photoelectron spectrum, we calculated the time-dependent Dyson norms for ionic states probed by our photon energy. The details of the 2D fitting analysis as well as the Dyson norm calculations are included in Appendices B 1 and C 1, respectively.

Figure. 2(a) shows the TRPES as a function of pump–probe delay and KE. As one can see, the TRPES shows two main peaks near time-zero. One broad peak has high KE between 1.2 eV and 4 eV, and a narrower peak has lower KE below ∼1 eV. As discussed in more detail below, these two peaks can be interpreted in terms of ionization to two different states of the cation. In order to follow the shifts of the two peaks, we extracted the peak positions for each delay and plotted these in Fig. 2(b). It is clear that the higher energy peak shifts systematically from about 2.5 eV to about 1.7 eV with an increasing pump–probe delay. This is in contrast to the lower energy peak, which does not shift significantly with the pump–probe delay, staying around 0.7 eV. Appendix B 2 contains a comparison of different approaches to quantifying the time-dependent shift in the peak positions.

In Fig. 2(c), we plot the calculated TRPES for each pump–probe delay. Fig. 2(d) shows the total ionization yield (energy integrated TRPES) for both the measurements and the calculations as a function of pump–probe delay, as well as the ground (S₀) and excited state (S₁ + S₂) populations. We interpret the measured and calculated spectra with the theoretical calculations in Sec. V.

In order to validate the statistical significance of features in the experimental measurements, a standard bootstrapping analysis was employed to estimate the uncertainties. One standard deviation (STD) is treated as an error bar, and this is indicated by the green and yellow shadings around the data points in Fig. 2(b). According to this analysis, the shift in the position of the high energy peak is about three times the standard deviation. Details of the bootstrapping analysis are included in Appendix B 2.

Our electronic structure and trajectory surface hopping calculations find that ultrafast decay of the excited state occurs through radiationless transitions to the ground state facilitated by conical intersections.⁵⁶ In general, there are three types of conical intersections between the ground and excited states, characterized by structural deformations around the two double bonds. Twisting and pyramidalization of a single double bond lead to two types of conical intersections, depending on whether the carbon adjacent to the other double bond or the carbon adjacent to the single bonds is pyramidalized. A third type of conical intersection involves twisting of both double bonds. The role of these conical intersections in the dynamics has been discussed in a separate theoretical study.⁵⁶ Here, we discuss how the Dyson norms along pathways connecting them to the initial FC point vary, as they determine the photoelectron yield as the wave packet evolves on the excited state.

Figure 3 shows which ionic states are produced when ionization from S₁ occurs based on Koopmans’ theorem: D₀ and D₃. These two states have a large Dyson norm with S₁, which can be explained by the fact that removing a single electron from either the HOMO or the LUMO leads to the dominant configuration in these two states. As also shown in the figure, ionization to D₃ and D₀ is energetically allowed, and thus, one expects ionization from S₁ to lead to both D₃ and D₀, at least near the FC point. In order to see how the energies of photoelectrons associated with ionization to these Dyson correlated states proceeds as the wave packet evolves on S₁, we calculated the energies of the neutral and ionic states between the FC point and the three different groups of CIs noted above.

Figure 4 shows the calculated electronic energies of the lowest lying neutral and ionic states at several points interpolated between the FC point and the three different CIs. Three important states are highlighted—the optically bright first excited state of the neutral, S₁ (red), the ground cationic state, D₀ (green), and the third excited state of the cation at FC, D₃ (cyan). We note that since D₃ crosses a number of ionic states en route to the third CI, it is of mixed character and is, therefore, labeled D_mix. In Fig. 4, the x-axis is the fraction from FC to the CI region.

It is clear from the figure that while the cyan and red lines are roughly parallel for all three panels as one moves away from the FC point, the green and red lines diverge. This means that one expects the low energy peak to remain roughly in the same place, while the high energy peak should shift to lower energies as the wave packet moves away from the FC point. This is consistent with our observation of a shifting high energy peak and a steady low energy peak in the photoelectron spectrum.

Figure 5 shows calculated photoelectron spectra along the LIICs by evaluating the Dyson norms between the neutral state S₁ and cation states along the paths to different CIs, as shown in Fig. 4. In the calculation, the photoelectron kinetic energy is obtained by subtracting the energy difference between D₀/D₃ and S₁ from the VUV probe photon energy in the experiment, which is 7.92 eV. The yield as a function of energy is then given by the norm of the Dyson orbital calculated by projecting each ionic state onto the neutral. In Fig. 5, one can see that the calculated photoelectron spectrum shows two bands from the FC point toward all three CIs, with the higher energy band decreasing in KE as a function of the fraction from FC to the CI. In contrast, the lower KE band maintains a relatively constant energy around 1 eV. The behavior of these two bands with the position along the LIIC is consistent with the calculated and measured time evolution of our photoelectron spectrum, showing two main peaks—one at higher energy, which shifts with position/delay, and the other at lower energy, which does not shift with position/delay.

The calculated energy for ionization to D₀ differs from the measured peak in the photoelectron spectrum for two main reasons. One is the error/uncertainty in the calculations, which is about 0.5 eV. The other is the fact that in the experiment, the pump laser excites the molecules on the red side of the absorption spectrum, meaning that the excitation is from the edges or tail of the ground state vibrational wavefunction on S₀ to lower vibrational levels on S₁. Since S₁ and D₀ diverge as one moves away from FC, the photoelectron spectrum is shifted to lower energies than one would calculate at FC. This leads to a lower measured photoelectron energy for ionization to D₀ than the calculations, although we note that the high energy shoulder of the measured high energy peak extends to roughly 5 eV (the calculated value), as one would expect based upon the explanation given above. In addition, we note that the measurements and calculations roughly agree on the location of the low energy peak (for ionization to D₃), consistent with the fact that D₃ is roughly parallel to S₁ near the Franck–Condon point, in contrast to the divergence of D₀ and S₁. A more detailed discussion on the calculated Dyson norms and what they reveal about the electronic structure of the neutral and cationic states is given in Appendix C 1.

The fact that the positions of the two peaks in the photoelectron spectrum vary differently with time delay, while the amplitudes of the peaks vary similarly with delay, is consistent with the fact that they both arise from lifting the same neutral wave packet on S₁ to different cationic states.

The measured and calculated photoelectron spectra shown in Figs. 2(a) and 2(c) agree on a number of features, but also have significant differences. They both contain two peaks in the spectrum near zero time delay (at ∼2.5 eV and ∼0.5 eV for the measurements and ∼4 eV and ∼1.5 eV for the calculations), with the high energy one shifting with delay and the low energy peak not shifting with delay. Both peaks decay on a 50 fs–100 fs time scale. However, the measurements and calculations disagree on the relative weightings of the two peaks and the exact time scale for the decay in the yield. In addition, the shift of the high energy peak in the calculations is such that the two peaks merge and appear as one for delays of 50 fs and greater. The shift of the high energy peak in the calculations is obscured to some extent by the convolutions with the experimental response function. The unconvolved calculation results show a clearer shift of the high energy peak, which are consistent with the results shown in Fig. 5. This level of agreement between calculations is consistent with other TRPES studies which compare measurements and calculations.^77–79 The disagreement on the relative positions of the peaks as a function of time is likely because we have introduced a constant shift for the gap between neutral and ionic states, assuming that the error in this gap is constant across the potential energy surfaces. This, however, is not correct, as we have shown by detailed comparisons in our previous work,⁵⁶ with the shifts becoming smaller as a function of time. It is, however, not possible to introduce a variable energy shift in our calculations.

As noted above, while the calculations and measurements agree on the qualitative features in the TRPES, they disagree on the decay times for the peaks. The experimental decay time agrees well with previous studies,^45,47 but it is about a factor of two faster than the theoretical decay time. Similar discrepancies have been observed before in ethylene.⁸⁰ In the earlier work, the discrepancy between the calculations and measurements was determined to be a result of windowing, since the probe pulse was not energetic enough to ionize the wave packet from everywhere on the excited state. While our probe photon energy is much larger than the ionization potential to D₀ everywhere on the excited state, ionization to higher lying states may be restricted as the wave packet evolves on S₁. This is explored in more detail in Appendix C 2. Another possibility that we explored is that the underlying electronic structure theory is not accurate enough. There are many points where the electronic structure theory can lead to slower dynamics. These include inaccurate energy gaps, slopes, and non-adiabatic couplings (NACs). In order to assess the influence of the level of theory on the decay dynamics, we compared our original CASSCF calculations with calculations at the XMS-CASPT2 level. A comparison of the state populations at these two levels of theory is also shown in Appendix C 2. The two methods give very similar results, so dynamical correlation does not seem to be the issue. It is possible that there are other effects, such as the basis set or the inclusion of Rydberg states that are missing and play a key role. Finally, we note that our calculations of the TRPES do not account for any of the probe photon energy going into vibrations. This has been accounted for in previous work^75,76 and would have the effect of lowering the energies of the peaks in the spectrum, which could impact the possible windowing effects discussed in more detail in  Appendices A–C. Before concluding, we note that calculations that we have carried out for a similar molecule, 1,3-cyclohexadiene, are in much better agreement with measurements reported in the literature⁷⁸ and measurements that we have performed. This suggests that the issues which lead to the differences between calculations and measurements for cc-COD are not systematic, but specific to cc-COD.

In conclusion, we have performed UV pump VUV probe measurements of cc-COD, using TRPES. The measurements are interpreted with the help of electronic structure and trajectory surface hopping dynamics calculations. The calculations allow us to assign the features in the TRPES and understand their behavior in terms of both the neutral and ionic state variations along the reaction coordinate. The calculations predict a slower decay than that measured by the experiment. Future work aims at addressing this discrepancy.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors are grateful to Chuan Cheng and Brian Kaufman for insightful discussions and technical help in the experiments described in this manuscript. The authors gratefully acknowledge support from the U.S. Department of Energy under Award Nos. DE-FG02-08ER15984 and DE-FG02-08ER15983. Most of the computational work was performed using the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation grant (No. ACI-1548562).

An important aspect of the pump–probe measurements is determining time-zero and the instrument response function (IRF). We made use of pump–probe measurements with ethylene in order to locate the time-zero and determine the IRF of the experiment given the very short decay time associated with VUV excitation. The ethylene pump–probe scans were conducted with the same running conditions as those for cc-COD. Figure 6 shows the ethylene TRPES signal and fitting results. Our previous measurements are consistent with earlier measurements in the literature and indicate that ethylene undergoes rapid internal conversion after VUV excitation, with a decay constant of 20 fs–40 fs.^12,60

Figure 6(a) shows the photoelectron spectrum for ethylene. Positive delay corresponds to the VUV pump and UV probe. One can see that the spectrum is slightly asymmetric due to the VUV induced relaxation dynamics. The energy integrated photoelectron yield was fit with a Gaussian convoluted with a single exponential function, as shown in Fig. 6(b). The center of the fit Gaussian function defines the time-zero for the experiment, with the full-width at half-max (FWHM) representing the IRF, i.e., the cross correlation of the UV and VUV pulses. Our fitting results yield a Gaussian FWHM ∼80 fs and an exponential decay of ∼26 fs.

In this section, we describe the detailed analysis of the measured photoelectron spectra. Figure 7(a) shows the photoelectron spectrum as a function of pump–probe delay with a linear intensity scale, including both negative and positive delays. One can see from the figure that there are dynamics on both sides of time-zero, which can complicate the interpretation for small positive delays, given the finite pump and probe pulse durations. The photoelectron yield is higher for negative delays than that for positive. This is a common feature for these pump and probe wavelengths, since many molecules have stronger absorption in the VUV (at ∼160 nm) than that in the deep UV (at ∼260 nm). In earlier work,²¹ we established that with an accurate determination of time-zero, we were able to extract information for the positive delays despite the larger signal for negative delays which can bleed into the positive delay signal because of the finite IRF.

In order to separate the VUV and UV driven dynamics, we carried out a 2D bilateral global least squares fit to the data.^35,81 Here, we define Y(t; ε_k) to be the differential photoelectron yield as a function of pump–probe (or probe–pump) delay and electron kinetic energy for the UV (or VUV) induced dynamics. In general, Y(t; ε_k) is globally fitted to

where the S_i(ε_k) are the time-independent decay related spectra and G(t) represents the Gaussian cross correlation function associated with the IRF. The energy-resolved amplitudes of the fitting components have decay constants τ_i, assuming that the population of the excited states follows an exponential decay. In Eq. (B1), a given S_i(ε_k) is associated with the energy-related photoionization cross section σ_i(ε_k) of the ith excited state. If the underlying dynamics involve photoionization from more than one excited state (i.e., i = 1, 2, …), one can establish kinetic models in the global fitting where more than one decay constant could be extracted from the fitting. The total electron yield is then convolved with G(t), the Gaussian cross correlation function that represents the IRF. In the case of cc-COD, we fit to only one exponential decay for the UV driven dynamics signals, since we are focusing on dynamics involving a single excited state. Thus, Eq. (B1) can be simplified as

This type of analysis of the 2D global fitting is based upon the assumption that the spectrum S₁(ε_k) is time-independent, and only the overall amplitude varies with time. In other words, the different portions of the photoelectron spectrum are fitted with the same decay constant, which means that spectral components associated with the dynamics of the relevant state remain frozen at the Franck–Condon point. Under this assumption, the molecule does not experience large amplitude motion away from the Franck–Condon point before the decay of the excited state, and such an analysis provides a reasonable estimation of the excited state lifetime.^82,83

In order to eliminate the VUV driven dynamics, we subtracted the fitted negative time signal from the total TRPES measurements.^35,81 In this way, we obtained the UV driven signal alone. This is important for the interpretation around zero time delay, where the pump and probe pulses overlap. We assumed a single exponential model for the VUV driven dynamics (negative time signal), leading to the following fitting function, which includes both positive and negative delays,

Figures 7(b) and 7(c) show the UV driven (pump–probe) and VUV driven (probe–pump) dynamics contributions to the TRPES, respectively, obtained by subtracting the fitted VUV and UV driven spectra from the measured spectra. We focus on the UV driven dynamics (pump–probe signal). With the subtraction of the fitted VUV driven dynamics (probe–pump signal), we are able to interpret the pump–probe signal for positive delays without the influence of the VUV excitation. As discussed earlier in connection with Figs. 2(a) and 2(b), the photoelectron yield of the higher energy peak (ε₂) displays a systematic shift toward lower energy for increasing positive pump–probe delays, whereas the lower energy peak (ε₁) does not shift. This feature is more obvious in Fig. 7(d) with a logarithmic intensity scale. In Fig. 7(d), we plot the time-dependent peak positions of ε₁ and ε₂ on top of the spectra of each delay step along the positive direction. Figure 7(e) shows the energy integrated electron yield from the measured and fitted spectra. It is clear that the signal from the VUV driven dynamics is stronger than the signal from the UV driven dynamics and contributes even for positive time delays out to ∼50 fs.

In addition to allowing for the subtraction of the VUV driven dynamics, performing a global 2D fit can serve as a test of whether or not there is a significant motion on the excited state before internal conversion. If there is no significant structure in the residuals associated with the fit, then this suggests that the excited state wave packet does not experience much displacement before internal conversion. On the other hand, the systematic structure in the residuals can serve as an indication or corroboration of motion on the excited state, which is viewed via the changes in the photoelectron spectrum as the wave packet moves away from the Franck–Condon point.

In the case of cc-COD, one of the features we observe for positive delays in the TRPES is the shift of the higher KE peak toward lower energy as a function of pump–probe delays, which we interpret in terms of the fact that S₁ and D₀ have the opposite slopes near the Franck–Condon point so that the wave packet motion along S₁ leads to a decreasing photoelectron energy with ionization to D₀. Thus, we expect to see systematic residuals in the 2D global fit. Figure 7(f) shows the residuals for the 2D global fitting results. The residuals show a systematic structure in multiple regions for positive delays. We focus our discussion on two key regions. One is a positive peak for early times over a broad range of energies (from ∼10 fs to 40 fs and from 0 eV to 5 eV). We interpret this feature in terms of an imperfect subtraction of the contributions from VUV driven dynamics, which can contribute while there is still some overlap of pump and probe pulses. The other is a narrow positive peak just below 2 eV between 50 fs and 150 fs, together with a broad negative peak between 2 eV and 4 eV over the same range of delays. This is consistent with a shift of the high energy peak in the spectrum to lower energy over these delays. The positive residuals for energies above 4 eV that do not vary significantly with time delay are likely due to stray electrons, which may also contribute to the positive residuals over all energies after 150 fs. Based upon the structures in the residuals, we conclude that the global 2D fit is consistent with the peak shift analysis above, indicating that our TRPES measurements provide evidence of a significant motion along S₁ before internal conversion.

In the TRPES measurements, the photoelectron spectra show two bands corresponding to ionization from the initial neutral state to different states of the cation. The higher kinetic energy peak shows a systematic shift toward lower energy, while the position of the lower kinetic energy peak stays relatively constant. We have carried out two approaches to follow the evolution of the peak locations. These approaches are based on finding the maximum for each peak and finding the center of mass (CoM)—from both the raw data and fitted spectra. We fit each of the two peaks shown in Fig. 2(b) to a single Gaussian function, which does not capture each peak perfectly, but allows for a simple analysis of the time-dependent location for each peak.

Figure 8 shows the time-dependent peak locations calculated based on finding the peak maxima and CoM from the raw and fitted data. Figures 8(a) and 8(b) show the time-dependent peak locations and CoM values from the raw spectra, while Figs. 8(c) and 8(d) plot the results from the fitted spectra. We first discuss the results from finding the locations of peak maxima. As shown in Figs. 8(a) and 8(c), the higher kinetic energy peak shifts toward lower energy monotonically, with a shift of ∼0.9 eV over 100 fs. In contrast, the lower kinetic energy peak oscillates and shifts less than 0.4 eV. Figures 8(b) and 8(d) plot the CoM values from both the raw and fitted spectra. The CoM of the higher energy peak also shows a monotonic shift toward the lower energy, although the shift is only about 0.4 eV. On the other hand, the CoM of the lower energy peak is roughly constant. The results indicate that there is a significant difference in the behavior of the lower and higher energy peaks, independent of the approach to quantifying the peak shifts.

A bootstrapping analysis was carried out in order to check the statistical uncertainty for the time-dependent kinetic energy changes. For each pump–probe delay, the photoelectron kinetic energy distribution was calculated from the averaging of 350 velocity-mapped images of the measurement. This dataset was re-sampled 100 times, with 350 images randomly selected each time, using the standard bootstrapping method. Each bootstrapped photoelectron kinetic energy distribution was analyzed using the full data analysis route, and the standard deviation from these 100 results was taken as the standard error. Specifically, we extracted the peak positions of the lower and higher energy peaks from the bootstrapped dataset and calculated the standard deviations as the error bar. This is indicated by the green and yellow shadings in Fig. 2(b).

A more careful examination of the calculated Dyson norms shows information about the changes in the character of both the neutral and cationic states. The neutral S₁ state initially corresponds mainly to a HOMO → LUMO excitation, but the character becomes more mixed along the relaxation pathways. In particular, mixing with the character of S₂ can occur, which corresponds to two different configurations, as illustrated in Fig. 3. Furthermore, the cationic states can become mixed and non-adiabatically change character along the pathways. A signature of these mixings is present in the third pathway involving the non-local CI. This is apparent when calculating the Dyson norms based only on the initial characters, as shown in Fig. 3. Using the simple correlations based on Koopmans’ theorem shown in that figure, it is predicted that the HOMO → LUMO excitation will lead to either a hole in the HOMO or a hole in the LUMO. The Dyson norms to those states are only shown in the bottom panels of Fig. 9. Comparing the top panels of that figure, which include all Dyson norms, to the bottom panels, which only include the two Koopmans’ correlated Dyson norms, we see that the first two paths generally agree, but there is a larger discrepancy for the third path, especially in the second half. This is a signature of increased mixing. The additional features come from ionization to the D₂ state, which has a mixed character, including the configurations (HOMO − 1)¹(HOMO)¹(LUMO)¹, (HOMO − 1)¹(HOMO)², and (HOMO − 1)¹(LUMO)². The first configuration can be produced by ionization from HOMO → LUMO (initial S₁ state), while the other two can be produced by ionizing from S₂. The non-local pathway distinguishes itself from the other two in more than one way. It has a stronger signature of mixing, as discussed above. This is similar to dynamics in butadiene, where the non-local CI is described as having more doubly excited character.³⁶ In addition, both the pathways, as shown in Fig. 4, and the calculated photoelectron spectra, shown in Fig. 5(c), are somewhat different from the other pathways. A barrier on the S₁ surface predicts that this pathway may be less dominant in the dynamics, and the trajectory surface calculations in our theoretical study⁵⁶ indicate that, indeed, it only plays a role in about 16% of the trajectories.

In addition to the TSH simulations at the CASSCF level, we also carried out TSH calculations at the XMS-CASPT2(4,4)/cc-pVDZ level to investigate if the level of theory is responsible for the discrepancy between the measured ionization yield and the calculated excited state population (and ionization yield). The same initial conditions (those used in the CASSCF TSH calculations) were used to simulate the absorption spectrum (with a Lorentzian line shape and a phenomenological broadening of 0.3 eV) and propagate the dynamics from the bright S₁ state. The XMS-CASPT2 absorption spectrum was plotted using an in-house code SArCASM,⁸⁴ whilst the dynamics was performed using Newton-X.⁸⁵ On-the-fly energies, gradients, and non-adiabatic couplings were generated using the Bagel package for the XMS-CASPT2 level.^86–88 Since both conformers showed similar dynamics and similar time scales at the CASSCF level, TSH at the XMS-CASPT2 level was performed for only the lowest energy conformer of cc-COD. The trajectories at the XMS-CASPT2 level were propagated starting from the S₁ state as it is the bright state in the FC region with the active space chosen for this study.⁵⁶ The fewest switches surface hopping (FSSH)⁸⁹ algorithm was employed to take into account non-adiabatic coupling (NAC) between the S₂, S₁, and S₀ states. Decoherence corrections were taken into account using the approach of non-linear decay of mixing by Granucci and Persico⁹⁰ with the recommended value of the empirical parameter, α = 0.1 Hartree.⁹¹ The velocity Verlet algorithm was used to integrate Newton’s equations of motion with a time step of 0.5 fs. The semiclassical time-dependent Schrödinger equation was integrated using fifth-order Butcher’s algorithm with a time step of 0.005 fs. The simulations were performed for 400 fs using XSEDE’s computational resources.⁹² 

Figure 10 compares the excited and ground state populations as functions of time for the TSH calculations at the CASSCF and XMS-CASPT2 levels of theory. It is clear from the figure that the two calculations produce very similar results. The figure also shows the population decay using only the trajectories corresponding to the excitation window used in the TRPES measurements. Again, it is clear that the population decay is similar to that of the trajectories corresponding to the whole absorption spectrum, suggesting that the dynamics are not very sensitive to exactly where the pump pulse is centered within the absorption spectrum.

In order to test the sensitivity of the calculations to the energies of the neutral and ionic states, and the extent of windowing effects that limit the experiments but are not captured by the calculations, we calculated the total ionization yield for different limits of integration in the photoelectron spectrum. This is shown in Fig. 11, which shows the calculated photoionization yield, integrating the calculated photoelectron spectrum between different energetic limits, together with the measured photoelectron yield. It is clear from the figure that one can produce calculated results that decay on time scales longer than or comparable to the experimental measurements, depending on what range of energies one integrates the photoelectron yield over. This highlights the sensitivity of the agreement to errors in the energies of neutral and ionic states as well as the ability to determine how much energy goes into the photoelectron vs vibrations during the photoionization process.