Ultrafast photodissociation dynamics of 2-ethylpyrrole: adding insight to experiment with ab initio multiple cloning

Abstract

The ultrafast photodissociation dynamics of 2-ethylpyrrole (2-EP) is simulated in a fully quantum manner on the S₁ and S₂ πσ* states by the ab initio multiple cloning (AIMC) method. AIMC treats electrons with accurate electronic structure methods “on the fly”, and nuclear dynamics with wavefunction propagation via a basis set of Ehrenfest trajectory guided Gaussian wavepackets. Total kinetic energy release (TKER) spectra are produced, as well as velocity map images and N–H dissociation times. These are compared to results from time-resolved velocity map imaging studies, and the AIMC method is able to provide quantitative reproduction of experimental data, including dissociation times of 50–80 fs. Novel insight into the dissociation mechanism is then obtained, with the experimentally obtained time constant shown to be composed of two components. Firstly, there is a contribution in <50 fs from 2-EP molecules that have sufficient energy in the N–H stretch coordinate to dissociate almost immediately over the barrier, and this is followed by a second slower contribution from 2-EP molecules that must sample the potential energy surface before finding a way around the barrier to dissociate. This two component mechanism is not observed experimentally due to the temporal widths of the laser pulses obscuring the dynamics in the <50 fs window, and is shown for the first time via theory. Calculations are also performed on selectively deuterated 2-EP, demonstrating that AIMC is able to produce a kinetic isotope effect for the dissociation time constant, and correctly predict a shift to lower energy in the TKER spectrum. The S₂ πσ* state is also shown to be unstable with respect to the S₁ πσ* state, with the N–H dissociation proceeding along S₁ when initially excited to S₂. This work demonstrates that the combination of state of the art theory and experiments can provide unprecedented novel insight into the N–H dissociation mechanism, with the tantalising prospect of providing insight into more general heteroatom hydride bond dissociation.

---

1 Introduction

A number of fundamental processes in chemistry and biology involve ultrafast excited state dynamics following photo-absorption, including light harvesting in plants,¹ UV photodamage in DNA,² and photoprotection in melanin pigments.³ Nitrogen containing aromatic heterocycles are found in molecules responsible for the above processes, and in particular pyrrole is a component of chlorophyll as well as being present in chromophores of other important natural compounds such as: vitamin B12, heme, bilirubin, biliverdin, and tryptophan. Due to the prevalence of pyrrole, in recent years there has been significant interest in its excited state dynamics both experimentally^4–20 and theoretically.^19–40

Whilst a large component of effort has been directed to understanding the excited state dynamics of pyrrole, much less well studied are its derivatives. This is also an important area to explore how modification alters its function. Recently, 2-ethylpyrrole (2-EP) has been studied experimentally by H (Rydberg) atom photofragment translational spectroscopy (HRA-PTS),³⁷ time-resolved velocity map imaging (TR-VMI) and time-resolved ion yield (TR-IY) mass spectrometry⁴¹ to examine the effects of ring-substitution on pyrrole. The only theoretical investigations on 2-EP have consisted of electronic structure calculations,³⁷ however dynamics calculations are desirable to further elucidate the photodissociation process. Fully quantum dynamics calculations still remain a huge challenge, however they are possible with ab initio multiple cloning (AIMC),^42,43 and in this work we use the method to offer insight to experiment on 2-EP.⁴¹ The effect of selective deuteration at the N–H bond has also been studied experimentally,⁴¹ so we consider modelling this as well to illustrate the capability of AIMC to reproduce experimental kinetic isotope effects.

Initially, it is worthwhile to briefly summarise the current understanding of the ultrafast excited state dynamics of pyrrole, and to contrast and compare to 2-EP. Early electronic structure work on pyrrole^21–24 helped to characterise the lowest energy excited states, with the seminal work by Sobolewski et al.^22,23 identifying the importance of the low lying πσ* states with regards to biological photoprotection. These states are dissociative along the N–H stretch coordinate, and provide an ultrafast radiationless transfer route through a conical intersection with the ground state, deactivating potentially reactive excited species. There are two low lying πσ* states in pyrrole that have spin symmetry labels 1¹A₂ and 1¹B₁ due to its C_2v molecular symmetry. Transitions from the X¹A₁ ground state to the lower energy of these two, the 1¹A₂(πσ*) state, are formally electric dipole forbidden. However, it can become directly populated through vibronic mixing with nearby higher lying ππ* states 1¹A₁ and 1¹B₂ (albeit with low transition cross-sections), or indirectly via internal conversion from the 1¹B₂(ππ*) state at shorter wavelengths. Electronic structure calculations for 2-EP revealed an increase in oscillator strength for the ground to first excited πσ* state transition compared to pyrrole,³⁷ which is to be expected due to the reduction in molecular symmetry to C_s in the anti conformation and C₁ in the lowest energy gauche conformer.

Early experimental studies on pyrrole observed two H atom dissociation channels in total kinetic energy release (TKER) spectra: one with a sharp high kinetic energy distribution of H atoms emitted perpendicular to the transition dipole moment; and one with a broader lower kinetic energy distribution, and a more isotropic emission.^4–9 It was posited that the former high kinetic energy dissociation channel was due to rapid N–H dissociation along the 1¹A₂(πσ*) state, followed by conical intersection with the ground state to produce a pyrrolyl radical in the 1²A₂ ground state. The latter lower kinetic energy dissociation channel was thought to be due to internal conversion to a vibrationally “hot” ground state, followed by dissociation. The high kinetic energy channel is believed to be dominant at longer pump wavelengths, with Wei et al.⁵ observing approximately 76% following this route at λ = 243.1 nm via the integration of the kinetic energy spectrum. At shorter pump wavelengths (λ < 218 nm) the lower kinetic energy channel is dominant,⁷ with 1¹B₂(ππ*) populated initially before passing through conical intersections and dissociating, with the pyrrolyl radical produced in the ground state. The addition of a σ donating ethyl group in 2-EP causes the onset of H atom dissociation to appear at longer wavelengths: 267 nm in 2-EP³⁷ compared to 254 nm in pyrrole.⁷ Excitation in the range 248 ≤ λ ≤ 263 nm leads to dissociation from the lowest energy πσ* state.³⁷

Time-resolved studies on pyrrole provided quantitative measures of the H atom appearance lifetimes in each of these dissociation channels.^15–20 Lippert et al. conducted the first of these studies by using TR-IY to obtain time-constants of τ₁ = 110 ± 80 fs for the high kinetic energy channel and τ₂ = 1.1 ± 0.5 ps for the lower kinetic energy channel at λ = 250 nm.¹⁵ Subsequently, Roberts et al. used both TR-VMI and TR-IY to observe a single time constant of 126 ± 28 fs for the high kinetic energy channel at λ = 250 nm.¹⁸ It was postulated that tunnelling out of the quasi-bound region of the 1¹A₂(πσ*) state plays a role at this excitation wavelength, as a time constant of 1.4 ± 0.3 ps for selectively deuterated pyrrole-d₁ dissociation gave a kinetic isotope effect (KIE) of ∼11. Further pump wavelengths were used by Roberts et al. in this study, with excitation at λ = 238 nm yielding a time constant of 46 ± 22 fs for undeuterated pyrrole, and 136 ± 38 fs for pyrrole-d₁. Due to the smaller KIE (∼3), this was attributed to faster, over the barrier dissociation with negligible contribution from tunnelling. Finally, at λ = 200 nm a time constant of 52 ± 12 fs was observed, attributed to initial population of the 1¹B₂(ππ*) state and subsequent rapid internal conversion to the 1¹A₂(πσ*) state and dissociation. A low kinetic energy feature was observed at this wavelength, however it had a time-constant of 1.0 ± 0.4 ns rather than the picosecond timescale observed by Lippert et al.¹⁵ It was suggested that the picosecond timescale observed by Lippert et al. could be due to undesired multiphoton dissociative ionisation events, and the nanosecond time-constant was due to C–H dissociation following internal conversion to the ground state, with the analogous process for N–H having an even longer timescale.¹⁸

TR-VMI and TR-IY studies on 2-EP in the 248 ≤ λ ≤ 265 nm range found high kinetic energy H atom appearance lifetimes on the order of 50–80 fs, and selectively deuterated 2-EP-d₁ at λ = 257 nm had an appearance lifetime of 140 ± 20 fs giving a KIE of ∼2, suggesting that tunnelling does not play a significant role in the dissociation mechanism. At excitation wavelengths shorter than 248 nm, an additional low kinetic energy feature appeared in the H atom TKER spectrum, similar to pyrrole. However, the appearance lifetime was ∼1.5 ps, intimating it was not due to dissociation following internal conversion to the ground state, and the explanation of a second dissociation channel involving population of higher-lying singlet states was offered once undesired multiphoton dissociation events and clustering in the molecular beam were ruled out.⁴¹

Previous quantum dynamics calculations on pyrrole have used reduced dimensionality potential energy surfaces,^25–27 surface hopping,^29–33 and multiconfigurational time-dependent Hartree (MCTDH) with parameterised potential energy surfaces.^{19,20,34–36} Our own ab initio multiconfigurational Ehrenfest (AI-MCE)⁴⁴ and ab initio multiple cloning^42,43 methods have also previously been used for dynamics studies on pyrrole.^38–40 In the former, using the AI-MCE method, conical intersections of the 1²A₂ and 1²B₁ radical states following H dissociation were observed, helping to explain the lack of experimentally observed pyrrolyl radicals in the 1²B₁ state.³⁸ In the latter, using the AIMC method, it was suggested that some low kinetic energy H atoms are formed in an ultrafast manner as a result of dissociation where the radical does not transfer to the ground state immediately.^39,40 Both methods use frozen/fixed width Gaussian basis functions (or equivalently coherent states) guided by Ehrenfest trajectories to describe the nuclear wavepacket, and calculate electronic potential energy surfaces at the centre of these Gaussian basis functions as the simulation is running (“on the fly”). This provides a fully dimensional and accurate description of nonadiabatic quantum dynamics.

Fully quantum simulations of multiatomic species still remain a huge challenge in chemical dynamics, and in the present work we seek to provide a full description of the ultrafast photodissociation dynamics of 2-EP using the AIMC method, for rigorous comparison to previously published experimental results.⁴¹ The computational details for the AIMC calculations carried out in this work are described in Section 3, and it is used to study the ultrafast dynamics of 2-EP in Section 4, alongside its selectively deuterated form. Initially however, a brief overview of the experimental method used to compare to the theoretical results is presented in the following section.

2 Experimental details

The experiments for 2-EP used TR-VMI apparatus with temporally delayed femtosecond pump and probe laser pulses to provide ultrafast time-resolved photochemical measurements. A full description of the apparatus used may be found in ref. 41.

The pump laser pulse was tuned over the wavelength range 238 ≤ λ_pump ≤ 265 nm and the probe pulse was produced at a wavelength of λ_probe = 243.1 nm to facilitate the 2 + 1 resonance enhanced multiphoton ionisation (REMPI) of H/D atom photofragments via the two photon allowed 2s ← 1s transition, and subsequent one photon ionisation. A cross-correlation/instrument response function at the temporal overlap of pump and probe pulses (Δt = 0) was obtained via the TR-IY of Xe⁺ (by 2 + 1′ ionisation). This measurement was fitted to a Gaussian distribution, which yielded a Gaussian instrument response function G_IRF(t) with width parameter σ_XC = 38 fs. The laser beams were then focused to intersect a molecular beam containing 2-EP, and the resulting cationic photoproducts accelerated down a time-of-flight tube and impacted onto a detector selectively gated for the exclusive detection of H⁺ (or D⁺) ions. The design and electric potentials of the ion optics and detector projected the three-dimensional photofragment velocity distribution onto the two-dimensional imaging detector, resulting in a velocity map image. From the velocity map image, the original three-dimensional photofragment velocity distribution was reconstructed and integrated using the polar onion peeling method,⁴⁵ resulting in one-dimensional photofragment velocity and angular distributions. After calibration to the known dissociation energetics of HBr at 200 nm,⁴⁶ the velocity distribution was transformed into a one-dimensional TKER spectrum. Multiple TKER spectra were recorded at various time delays Δt to produce H/D transients, which were used to obtain appearance lifetimes.

3 Computational details

As mentioned in the introduction, the AIMC method uses frozen Gaussian basis functions guided by Ehrenfest trajectories to describe the nuclear wavepacket, and calculates electronic potential energy surfaces at the centre of these Gaussian basis functions “on the fly” to provide a fully quantum description of nonadiabatic dynamics. The full numerical details of the method has been described elsewhere,^42,43 and brief overview of the working equations is given in the ESI.† For the discussion in the main text, it suffices to say that the Ehrenfest trajectories that guide the Gaussian basis incorporate components from all electronic states involved in the calculation, and hence include nonadiabatic effects at all times. In regions of strong nonadiabatic coupling when an Ehrenfest trajectory is comprised of an unphysical average of multiple electronic states, a “cloning” procedure is applied to expand the basis, mimicking wavepacket splitting. This usually occurs following a nonadiabatic transition with incomplete population transfer, when the electronic states with significant population also have differing gradients.

As with previous AIMC studies,^39,40,42 dynamics were simulated using a modified version of AIMS-MOLPRO^47–49 that incorporates Ehrenfest trajectories. Electronic structure calculations were performed using the complete active space self-consistent field (CASSCF) method. The electronic basis set used was Dunning's cc-PVDZ set,⁵⁰ with one additional diffuse s function, one additional set of p functions, and one additional set of d functions added to the nitrogen atom; as well as one additional diffuse s function, and one additional set of p functions added to the dissociative hydrogen atom. The active space used has 8 electrons in 7 orbitals: three ring π orbitals and two corresponding π* orbitals, and the N–H σ and corresponding σ* orbital. State averaging was performed over three states and dynamics were performed on three states: the ground and two lowest excited singlet states. The width of the Gaussian basis functions γ was taken to be 4.7 Bohr⁻² for hydrogen, 6.6 Bohr⁻² for deuterium, 22.7 Bohr⁻² for carbon and 19.0 Bohr⁻² for nitrogen, as suggested previously.⁵¹

Initial positions and momenta for the nuclei were sampled from the ground state vibrational Wigner distribution in the harmonic approximation using vibrational frequencies and normal modes calculated at the same level of CASSCF theory as above. As in previous AIMC works,^39,40,42 excitation from ground to excited state is approximated by projecting the wavefunction onto the desired excited state, which amounts to setting the normalised electronic state amplitude for the desired excited state to unity for each trajectory, whilst the other electronic state amplitudes are set to zero. It is noted that the finer details of initial photoexcitation are not accounted for completely by this approximation, and there may be some small energetic sampling deficiencies that are discussed in the following section. The temporal widths of the pump and probe laser pulses from the photoexcitation are taken into account when calculating time-constants for H dissociation however, the details of which shall also be explained in the following section.

The initial number of trajectories, starting adiabatic electronic state, and number of cloning events for each calculation carried out in this work are shown in Table 1. Undeuterated 2-EP calculations are performed with dynamics starting on the S₁ and S₂ states, and deuterated 2-EP calculations are performed with dynamics starting on only the S₁ state. Hereafter, they are referred to with the starting adiabatic electronic state in parenthesis afterwards. The parameters used to determine a cloning event are the same for each system, and are explained in more detail in the ESI.† The smaller number of cloning events for 2-EP(S₂) compared to 2-EP(S₁) is due to more complete population transfer at conical intersections. The smaller number of cloning events for 2-EP-d₁(S₁) compared to 2-EP(S₁) is due to the lower frequency of the N–D vibration compared to N–H, meaning fewer conical intersections are encountered along this coordinate. Furthermore, the cloning procedure incorporates a description of wavepacket splitting which is less likely to occur for D than H as it is heavier. Calculations were run for 350 fs, or until the N–H bond exceeded 4 Å, which was defined as the point of dissociation. All calculations used a timestep of ∼0.6 fs (2.5 a.u.).

Table 1 Molecule labels, initial number of trajectories, starting electronic state, number of cloning events for each calculation carried out in this work, calculation time, and percentage of trajectories dissociated by the end of the calculation time

Molecule label	Initial trajectories	Initial adiabatic electronic state	Cloning events	Calculation time (fs)	Percentage dissociated
2-EP(S1)	600	S1	142	350	62%
2-EP(S2)	600	S2	51	350	81%
2-EP-d1(S1)	600	S1	18	350	55%

---

4 Results

In the following section, calculated H/D atom TKER spectra are presented and compared to experiment. The dissociation energies obtained from calculation are smoothed with Gaussian functions (with width σ = 200 cm⁻¹) to create a curve for the TKER spectrum, as opposed to a stick spectrum of delta functions. This allows better comparison to experiment, and this procedure has been used before to obtain simulated TKER spectra.^39,40 An experimental velocity map image relative to the electric component of the laser field is reproduced by the 2-EP(S₁) calculation using a calculated velocity distribution relative to the molecular axes.

Calculated dissociation times are also presented, in terms of raw dissociation times from the trajectory data, and following a smoothing procedure to take into account experimental laser pump and probe temporal widths (this procedure has also been carried out recently⁵²). The smoothing procedure is required because the simulation has well defined start and end points, however, experimentally these points are “blurred” by the temporal widths of the laser pulses. The process for this smoothing procedure is as follows: the pump and probe laser pulses are assumed to be Gaussian in time, with widths obtained from the experimental G_IRF(t) measurement (which has width σ_XC, equal to the cross-correlation/convolution of pump and probe pulses). The pump and probe laser pulses are both produced in the same manner, so it is expected that their temporal widths will be similar. However, as the probe consists of a 2 + 1 REMPI of H atom photofragments, this simultaneous absorption of 3 photons will reduce the width of the Gaussian that models it by a factor of so that . Therefore we have σ_pump = 33 fs and fs, which reproduce well the experimentally obtained σ_XC = 38 fs.

The raw dissociation times from the trajectory data are then represented as Gaussian probability distributions, centred around each raw dissociation time, with width parameter σ_pump. For each point in this distribution, the dissociating trajectory will have a probability of being “probed” by another Gaussian probability distribution with width σ_probe. These Gaussian probability distributions are convoluted and summed for each trajectory with associated weights from calculated amplitudes, enabling a smoothed transient to be produced that can be compared to experimental data. This smoothed transient is then fit using the kinetic model employed in experiment^18,41 to obtain time constants, and these fits are analysed relative to the raw dissociation times to add insight to experimental results. The kinetic model is given by an exponential rise with H/D-atom appearance lifetime τ and time zero correction t₀, convoluted with the Gaussian instrument response function G_IRF(t), multiplied by amplitude A and Heaviside unit step function u(t) (which is equal to 1 for t ≥ 0 and 0 otherwise)

S(t) = A·G_IRF(t) × [(1 − e^{−(t−t₀)/τ})u(t)]. (1)

4.1 TKER spectra

Calculated TKER spectra are shown in Fig. 1, with experimental TKER spectra shown in the insets,⁴¹ and intensity of the main peaks maxima set to 1. The theoretical spectra show good qualitative agreement to those obtained by experiment, with one main peak and some smaller features at lower kinetic energies. Comparing the spectra from S₁ and S₂ dynamics in Fig. 1a, a broader profile of the main peak from 2-EP(S₂) than 2-EP(S₁) is observed. This is also seen experimentally upon decreasing the pump wavelength, where the effect of going from a 262 to 252 nm pump pulse is illustrated in the inset. This indicates that excitation at shorter wavelengths may involve some population of the S₂ state initially, rather than just the S₁ state. Quantitatively, the energy of the main peak is 3000–4000 cm⁻¹ larger from simulation than experiment, due to inaccuracies in the CASSCF potential energy surface. To attempt to quantify these inaccuracies, we performed multireference perturbation theory calculations with multistate corrections (MS-CASPT2) at the Franck–Condon region of the S₁ and S₂ potential energy surfaces, and the asymptotic region of the S₀ state (when the N–H distance is equal to 4 Å). We calculated the difference in potential energy between these two regions, and noted that CASSCF overestimated the difference by 3750 cm⁻¹ for S₁ and 2228 cm⁻¹ for S₂ compared to MS-CASPT2. This causes the shift in peak in the TKER spectra, and if it were computationally feasible to use a MS-CASPT2 potential energy surface for the dynamics, then the energies would be much more comparable to experiment. A similar effect has been noted previously with pyrrole.^38,39 There are also slightly larger shoulders on the high energy side than experiment that can be ascribed to the sampling procedure producing some artificially high energy trajectories. This issue may be remedied by a new sampling procedure, recently developed, that simulates the initial photoexcitation by taking into account the laser pump wavelength and pulse shape. Initial results using this procedure appear to confirm this hypothesis.⁵³

Fig. 1 (a) Calculated H atom TKER spectra for 2-EP(S₁) and 2-EP(S₂), with experimental results at a pump wavelengths of λ_pump = 262 nm and λ_pump = 252 nm in the inset.⁴¹ (b) Calculated H atom TKER spectrum for 2-EP-d₁(S₁) overlaid with that of 2-EP(S₁), and experimental results at a pump wavelength of λ_pump = 257 nm in the inset.⁴¹

The TKER spectrum for 2-EP-d₁(S₁) is shown in Fig. 1b, with the 2-EP(S₁) TKER spectrum overlaid for comparison. An experimental spectrum for 2-EP-d₁ at a pump wavelength of 257 nm is shown in the inset, with undeuterated 2-EP at the same pump wavelength overlaid.⁴¹ It can be seen that AIMC correctly reproduces the shift of the main peak to lower kinetic energies, due to the difference in zero point energies of the deuterated and undeuterated forms of 2-EP. Both parts of Fig. 1 therefore demonstrate the ability of AIMC to reproduce spectral features observed experimentally.

When comparing the 2-EP(S₁) TKER spectrum to that of pyrrole previously calculated in an analogous manner by AIMC,⁴⁰ we observe a shift in the main peak to lower energy for 2-EP compared to pyrrole (∼10 000 cm⁻¹ for 2-EP(S₁) and ∼11 000 cm⁻¹ for pyrrole). This is also seen experimentally (∼6600 cm⁻¹ for 2-EP and ∼7200 cm⁻¹ for pyrrole¹⁸), where the overestimation of the simulated peaks is due to inaccuracies of the CASSCF potential energy surface as noted above.

4.2 Velocity map images

Fig. 2a shows an experimental velocity map image at λ_pump = 262 nm, obtained via reconstruction of a slice through the centre of the original three-dimensional photofragment velocity distribution using the polar onion peeling method.⁴⁵ It demonstrates the ejection of H atoms relative to the electric component of the laser field ε, which is mapped onto the molecular axes via the transition dipole moment. It can be seen that the high kinetic energy H atoms are predominantly emitted perpendicularly to the electric field, and that the lowest kinetic energy H atoms are emitted relatively isotropically, albeit this is a less intense feature of the velocity map image, as with the TKER spectrum.

Fig. 2 (a) Experimental velocity map image at λ_pump = 262 nm. (b) Simulated velocity map image for 2-EP(S₁) with respect to laser pulse polarisation, using calculated transition dipole moments. (c) Definition of molecular axes for 2-EP with an illustration of the S₀ to S₁ transition dipole moment. (d) Three-dimensional velocity distribution of ejected H atoms from calculation, with points coloured according to their density. A projection of this point density on the xy plane is shown via a contour map.

A simulated velocity map image may be produced to further validate the theory compared to experiment by averaging the calculated velocity distribution of hydrogen atoms relative to the axes of 2-EP over all possible orientations of the molecule. Integrating over Euler angles and using a narrow Gaussian function for the velocity vector |v| as opposed to a δ function, we obtain^39,40

where θ is the angle between |v| and the transition dipole of 2-EP, and σ = 0.001 Å fs⁻¹ is the width of the Gaussian function. The definition of the molecular axes for 2-EP is shown in Fig. 2c, along with an illustration of the S₀ to S₁ transition dipole moment at the optimised ground state geometry calculated by CASSCF at the same level of theory as the dynamics calculations. The transition dipole moment is sensitive to geometry, therefore it was calculated by CASSCF for the initial geometry associated with each trajectory for use in eqn (2) to obtain the simulated velocity map image, shown in Fig. 2b. One can see that the direction of the main anisotropic high kinetic energy component of the experimental velocity map image is reproduced well by the simulations, albeit the energetics are overestimated as discussed in the previous section. The low kinetic energy components are less well resolved due to lower statistics, however there is a small component of isotropic low kinetic energy hydrogen atoms seen in the calculated velocity map image.

The calculated velocity distribution used to obtain the simulated velocity map image is shown in Fig. 2d as a three-dimensional scatter plot of the velocities of ejected H atoms for each trajectory. The points are coloured according to the density of dissociation velocities, and a two-dimensional projection of this onto the xy plane is shown. This permits straightforward analysis of the dissociation direction relative to the molecular axes of 2-EP. It can be seen that the H atoms have velocities with significant z component, i.e. along the N–H stretching coordinate, and from the projection of the density of velocities onto the xy plane is can be seen that there is some preference for the negative x and y directions, i.e. away from the ethyl group. The small contribution of trajectories to the low kinetic energy isotropic component of Fig. 2b may also be seen from those with low velocity in the z-direction. Analysis of these individual trajectories reveals that they undergo cloning events and retain some character on the excited states as dissociation occurs. This effect has also been noted to be responsible for the similar weak intensity low kinetic energy portion of the pyrrole TKER spectrum,³⁹ and the radical formed does eventually relax to the ground state.³⁸ This individual analysis of trajectories is possible because they are run separately and combined at the end of the calculation, see Section S1 in the ESI† for further details.

A greater isotropic emission of low kinetic energy H atoms appears experimentally when a second N–H dissociation channel appears at λ_pump < 248 nm, with a longer (∼1.5 ps) time constant.⁴¹ This feature begins to appear predominantly in the TKER spectrum after 200 fs, and is the main source of ejected H atoms at λ_pump = 238 nm. The calculations showed no large increase in the low kinetic energy feature in the 200–350 fs range, and is unlikely to appear if the calculations were continued beyond 350 fs, as the majority of trajectories have dissociated by this point. This, combined with the fact that the absorption edge of the higher lying ππ* states begins to appear at λ_pump < 248 nm,⁴¹ suggests that initial excitation to the S₃ or S₄ ππ* states occurs at these pump wavelengths rather than S₁ (or S₂). From that point, the mechanism of the low kinetic energy H dissociation is unclear, as currently dynamics calculations involving additional excited states (hence greater active space for the CASSCF calculation) for longer times are prohibitively expensive for our AIMC method. However, it is possible that the mechanism will involve 2-EP molecules that retain some S₃/S₄ excited state character as they dissociate, similar to the much less intense low kinetic energy feature seen in 2-EP(S₁) dynamics, where the trajectories retain some S₁/S₂ excited state character as they dissociate.

4.3 Dissociation kinetics

Turning to the kinetics of dissociation, a cumulative sum of raw dissociation times, alongside smoothed H/D atom appearance transients with associated fits, and experimental data are presented in Fig. 3 for 2-EP(S₁), 2-EP(S₂), and 2-EP-d₁(S₁). The smoothing procedure was performed as outlined in the introduction to this section, taking into account the temporal widths of the pump and probe laser pulses. The smoothed transients obtained from these probability distributions were then fit using the kinetic model in eqn (1). A time zero correction was applied to the experimental data as determined in ref. 41, whilst it is not necessary (t₀ = 0) for the AIMC data. The experimental data show non-zero baseline at negative times due to “reverse dynamics”. In these reverse dynamics, the probe pulse acts as a pump to photoexcite 2-EP, and then subsequently provides two photons in a 2 + 1′ REMPI scheme for the dissociated H atom, with the third photon provided by the “pump” pulse. This phenomenon has been explained previously in ref. 15 and 18.

Fig. 3 Raw cumulative sum of dissociation times from trajectories, alongside smoothed H atom appearance transients with associated fits and experimental data for (a) 2-EP(S₁) with experimental data at a pump wavelength of λ_pump = 262 nm, (b) 2-EP(S₂) with experimental data at a pump wavelength of λ_pump = 252 nm, (c) 2-EP-d₁(S₁) with experimental data at a pump wavelength of λ_pump = 257 nm.

Visually, it can be seen that the smoothed transients compare well to experimental data, albeit with a slight lifetime shift for 2-EP-d₁(S₁). The lifetimes obtained of 68.2 ± 0.5 fs, 72.6 ± 0.4 fs, and 104.6 ± 1.1 fs for 2-EP(S₁), 2-EP(S₂), and 2-EP-d₁(S₁), respectively, also compare well to the experimentally obtained 55 ± 13 fs and 70 ± 20 fs at λ_pump = 262 nm and 252 nm for 2-EP, and 140 ± 20 fs at λ_pump = 257 nm for 2-EP-d₁. The KIE from simulation of ∼1.5 is slightly less than the experimentally obtained KIE of ∼2, however if the calculation was extended for a longer period of time so that more trajectories dissociate it is expected that this value will increase as undissociated trajectories are still in the S₁ state. Furthermore, there may be some tunnelling effects that are not accounted for by the calculation which would give rise to the larger experimental KIE. The effect of tunnelling with regards to dissociation energetics has been explored by calculation previously with pyrrole,⁴⁰ however it has not been considered with regards to dissociation kinetics, and this could be a subject of a future avenue of research to determine the influence on the calculated KIE, with the expectation that the KIE would increase to more closely match experiment.

From the lifetimes obtained from the fit, it would appear that the kinetics of dissociation are similar from the S₁ and S₂ surfaces. However, when considering the proportion of total trajectories that are dissociated by 350 fs (shown in Table 1) it is observed that a greater proportion have dissociated when dynamics are started on the S₂ surface (81%) than the S₁ surface (62%). The reason for both of these effects may be seen from the averaged electronic state populations for 2-EP(S₂) trajectories in the first 50 fs, illustrated in Fig. 4. This figure shows that there is an immediate and significant transfer of population from S₂ to S₁ as the calculation begins, indicating that the S₂ state is unstable with respect to S₁. Around 10 fs later some trajectories reach the S₁–S₀ conical intersection, there is steady population transfer, and trajectories begin to dissociate. The immediate S₂ to S₁ transfer produces trajectories in the S₁ state with high energy, and is the reason why a larger percentage of 2-EP(S₂) trajectories dissociate compared to 2-EP(S₁). Furthermore, this may also explain the broader profile of the main peak in the 2-EP(S₂) TKER spectrum compared to 2-EP(S₁) in Fig. 1a, as the higher energy trajectories result in a broader distribution of kinetic energies for the emitted H atoms.

Fig. 4 Averaged electronic state populations for 2-EP(S₂) trajectories.

Returning to the transients in Fig. 3, further insight into the dissociation kinetics may be obtained by considering the raw cumulative sum of dissociation times from the trajectories. Examining these a few things may be noticed, particularly in the sub 50 fs regime, that are masked by the temporal widths of the laser pulses experimentally. Firstly, no trajectories dissociate until 14.6 fs for 2-EP(S₁), 12.3 fs for 2-EP(S₂), and 24.8 fs for 2-EP-d₁(S₁). This is merely a consequence of defining the point of dissociation as 4 Å, with the delay the time taken for the N–H/D bond to stretch to this distance. Experimentally this will also occur, as the N–H bond extends over the πσ* surface following the pump laser pulse. However, lack of temporal resolution does not permit this to be observed in the experimental transient, and the exact point of dissociation is less easy to define. The longer time for 2-EP-d₁(S₁) trajectories to begin to dissociate compared to 2-EP(S₁) is due to the lower vibrational frequency of the N–D bond compared to N–H.

More interestingly, from the initial dissociation point onwards for the next ∼40 fs a rapid increase in the number of trajectories dissociating is observed, as those prepared in geometries with the correct orientation and sufficient energy in the N–H stretching coordinate to dissociate over the barrier do so either immediately or following a vibrational period. Fig. S1 in the ESI† illustrates the tendency for trajectories with a large amount of energy in the N–H stretching coordinate to dissociate rapidly.

The majority of trajectories that dissociate within the calculation time do so by this mechanism in the first 50 fs for 2-EP(S₁) and 2-EP(S₂), and in the first 70 fs for 2-EP-d₁(S₁). This rapid increase of dissociating molecules is not seen to the same extent experimentally, instead there is a much smoother rise due to the “blurring” effect of the laser pulses. Following this initial rapid rise, the rate of dissociation slows for the rest of the calculation as the remaining trajectories do not have enough energy in the stretching coordinate to immediately dissociate, and must first sample more of the potential energy surface to find a way around the barrier. See Fig. S2 (ESI†) for an example dynamic potential energy of a trajectory that (a) dissociates via the fast mechanism, and (b) dissociates via the slow mechanism. It should be noted that the initial rapid dissociation process may be slightly overestimated by the sampling technique producing too many high energy trajectories, as noted in Section 4.1 with the additional high energy shoulder on the calculated TKER spectrum relative to experiment, and this is a further possibility of the cause the time shift of the 2-EP-d₁(S₁) transient relative to experimental results. The recently developed sampling procedure also mentioned in Section 4.1 may remedy this.⁵³

5 Conclusions

We have used AIMC to simulate the ultrafast photodissociation of 2-EP in a fully quantum manner with dynamics starting on the S₁ and S₂ πσ* states, and deuterated 2-EP with dynamics starting on the S₁ state. TKER spectra and H/D-atom appearance lifetimes from the N–H/D dissociation have been produced and compared to time-resolved velocity map imaging experimental results.⁴¹ The TKER spectra reproduce the structure of the main peak and less intense low kinetic energy features from experiment. Comparative features observed experimentally have also been reproduced in each of the spectra, such as the shift of the main peak to lower energies due to deuteration, the lower energy of the peak for 2-EP compared to previous AIMC calculations on pyrrole,⁴⁰ and the broadening of the main peak due to shorter pump wavelengths. The latter effect was observed in the 2-EP(S₂) calculations compared to 2-EP(S₁), and was explained by considering the averaged electronic state populations of 2-EP(S₂) trajectories. These populations demonstrated that the S₂ state is unstable with respect to S₁, leading to rapid transfer from S₂ to S₁ that produced trajectories in a high energy S₁ state and led to a broad distribution of kinetic energies for the dissociated H atoms.

Velocity distributions were also calculated for 2-EP(S₁), illustrating predominant dissociation in the direction of the N–H bond, but with a slight preference for the H atom to be ejected away from the ethyl group. There was also a small contribution of low kinetic energy isotropically distributed H atoms, due to trajectories that retain excited state character upon dissociation, before the radical eventually relaxes to the ground state.

The H/D atom appearance lifetimes were obtained by applying a smoothing procedure to the raw dissociation times obtained from calculation. This smoothing procedure took into account the experimental pump and probe laser pulse temporal widths to allow a direct comparison between simulation and experiment. Visually, these smoothed transients compared extremely well to the experimental results, and the time constants obtained from fitting to the kinetic model used in the experiments were also in very good agreement. Insight into the kinetic model and sub-50 fs dynamics was achieved by considering the raw dissociation time data, and how it compared to the smoothed transients and experimental data. The raw dissociation time data showed that the dissociation was essentially a two-step process: firstly, 2-EP molecules with the correct geometry to dissociate over the barrier do so either immediately or following a vibrational period; and secondly, this is followed by a slower rate of dissociation for the 2-EP molecules that must sample more of the potential energy surface before finding a way around the barrier. This two-step process is “blurred” by the temporal width of the pump and probe laser pulses, leading to a single lifetime for the smoothed transients and experimental data, and it is observed for the first time by the AIMC method.

Performing quantum simulations of multiatomic species still remains a significant challenge, and in this work AIMC has been applied to its largest system yet, demonstrating its ability to reproduce experimental quantities and add insight into them. Further calculations into molecules of photochemical interest are either currently underway or planned; including pyrazole, imidazole, and aniline. The method is not limited to photochemistry however, and can be broadly applicable to a great many small molecule systems, with additional insight into fundamental ultrafast processes on the horizon.

Conflicts of interest

There are no conflicts of interest to declare.

Acknowledgements

The work has been supported by the Leverhulme trust RPG-2015-190. J. A. G. is currently supported by EPSRC grant EP/N007549/1 and previously has been supported by the University Research Scholarship from the University of Leeds, as well as funding from the School of Chemistry, University of Leeds. D. M. is supported by EPSRC grant EP/P021123/1 and also would like to acknowledge the previous support by the EPSRC grant EP/N007549/1. N. C. C. F. thanks the Leverhulme Trust for postdoctoral funding. V. G. S. thanks the EPSRC for equipment grants (EP/J007153 and EP/N010825) and the Royal Society and Leverhulme Trust for a Royal Society Leverhulme Trust Senior Research Fellowship.

References

1. Y. C. Cheng and G. R. Fleming, Annu. Rev. Phys. Chem., 2009, 60, 241.
2. R. Improta, F. Santoro and L. Blancafort, Chem. Rev., 2016, 116, 3540.
3. L. A. Baker, B. Marchetti, T. N. V. Karsili, V. G. Stavros and M. N. R. Ashfold, Chem. Soc. Rev., 2017, 46, 3770.
4. D. A. Blank, S. W. North and Y. T. Lee, Chem. Phys., 1994, 187, 35.
5. J. Wei, A. Kuczmann, J. Riedel, F. Renth and F. Temps, Phys. Chem. Chem. Phys., 2003, 5, 315.
6. J. Wei, J. Riedel, A. Kuczmann, F. Renth and F. Temps, Faraday Discuss., 2004, 127, 267.
7. B. Cronin, M. G. D. Nix, R. H. Qadiri and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2004, 6, 5031.
8. B. Cronin, A. L. Devine, M. G. D. Nix and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2006, 8, 3440.
9. L. Rubio-Lago, D. Zaouris, Y. Sakellariou, D. Sofikitis, T. N. Kitsopoulos, F. Wang, X. Yang, B. Cronin, A. L. Devine, G. A. King, M. G. D. Nix, M. N. R. Ashfold and S. S. Xantheas, J. Chem. Phys., 2007, 127, 064306.
10. M. N. R. Ashfold, B. Cronin, A. L. Devine, R. N. Dixon and M. G. D. Nix, Science, 2006, 312, 1637.
11. M. N. R. Ashfold, G. A. King, D. Murdock, M. G. D. Nix, T. A. A. Oliver and A. G. Sage, Phys. Chem. Chem. Phys., 2010, 12, 1218.
12. G. M. Roberts and V. G. Stavros, Chem. Sci., 2014, 5, 1698.
13. S. W. Crane, M. M. Zawadzki, J. O. F. Thompson, N. Kotsina, O. Ghafur and D. Townsend, J. Chem. Phys., 2016, 145, 234304.
14. S. L. Horton, Y. Liu, P. Chakraborty, S. Matsika and T. Weinacht, J. Chem. Phys., 2017, 146, 064306.
15. H. Lippert, H.-H. Ritze, I. V. Hertel and W. Radloff, ChemPhysChem, 2004, 5, 1423.
16. R. Montero, A. P. Conde, V. Ovejas, M. Fernández-Fernández, F. Castaño, J. R. Vázquez de Aldana and A. Longarte, J. Chem. Phys., 2012, 137, 064317.
17. R. Montero, V. Ovejas, M. Fernández-Fernández, A. P. Conde and A. Longarte, J. Chem. Phys., 2014, 141, 014303.
18. G. M. Roberts, C. A. Williams, H. Yu, A. S. Chatterley, J. D. Young, S. Ullrich and V. G. Stavros, Faraday Discuss., 2013, 163, 95.
19. G. Wu, S. P. Neville, O. Schalk, T. Sekikawa, M. N. R. Ashfold, G. A. Worth and A. Stolow, J. Chem. Phys., 2015, 142, 074302.
20. O. M. Kirkby, M. A. Parkes, S. P. Neville, G. A. Worth and H. H. Fielding, Chem. Phys. Lett., 2017, 683, 179.
21. M. H. Palmer, I. C. Walker and M. F. Guest, Chem. Phys., 1998, 238, 179.
22. A. L. Sobolewski and W. Domcke, Chem. Phys., 2000, 259, 181.
23. A. L. Sobolewski, W. Domcke, C. Dedonder-Lardeux and C. Jouvet, Phys. Chem. Chem. Phys., 2002, 4, 1093.
24. B. O. Roos, P.-Å. Malmqvist, V. Molina, L. Serrano-Andrés and M. Merchán, J. Chem. Phys., 2002, 116, 7526.
25. V. Vallet, Z. Lan, S. Mahapatra, A. L. Sobolewski and W. Domcke, Faraday Discuss., 2004, 127, 283.
26. V. Vallet, Z. Lan, S. Mahapatra, A. L. Sobolewski and W. Domcke, J. Chem. Phys., 2005, 123, 144307.
27. Z. Lan, A. Dupays, V. Vallet, S. Mahapatra and W. Domcke, J. Photochem. Photobiol., A, 2007, 190, 177.
28. M. Barbatti, M. Vazdar, A. J. A. Aquino, M. Eckert-Maksić and H. Lischka, J. Chem. Phys., 2006, 125, 164323.
29. M. Vazdar, M. Eckert-Maksić, M. Barbatti and H. Lischka, Mol. Phys., 2009, 107, 845.
30. B. Sellner, M. Barbatti and H. Lischka, J. Chem. Phys., 2009, 131, 024312.
31. M. Barbatti, J. Pittner, M. Pederzoli, U. Werner, R. Mitrić, V. Bonačić-Koutecký and H. Lischka, Chem. Phys., 2010, 375, 26.
32. M. Barbatti and K. Sen, Int. J. Quantum Chem., 2016, 116, 762.
33. M. Sapunar, A. Ponzi, S. Chaiwongwattana, M. Malis, A. Prlj, P. Decleva and N. Doslic, Phys. Chem. Chem. Phys., 2015, 17, 19012.
34. S. P. Neville and G. A. Worth, J. Chem. Phys., 2014, 140, 034317.
35. H. Köppel, E. Gromov and A. Trofimov, Chem. Phys., 2004, 304, 35.
36. S. Faraji, M. Vazdar, V. S. Reddy, M. Eckert-Maksic, H. Lischka and H. Köppel, J. Chem. Phys., 2011, 135, 154310.
37. T. N. V. Karsili, B. Marchetti, R. Moca and M. N. R. Ashfold, J. Phys. Chem. A, 2013, 117, 12067.
38. K. Saita, M. G. D. Nix and D. V. Shalashilin, Phys. Chem. Chem. Phys., 2013, 15, 16227.
39. D. V. Makhov, K. Saita, T. J. Martínez and D. V. Shalashilin, Phys. Chem. Chem. Phys., 2015, 17, 3316.
40. D. V. Makhov, T. J. Martínez and D. V. Shalashilin, Faraday Discuss., 2016, 194, 81.
41. N. C. Cole-Filipiak, M. Staniforth, N. d. N. Rodrigues, Y. Peperstraete and V. G. Stavros, J. Phys. Chem. A, 2017, 121, 969.
42. D. V. Makhov, W. J. Glover, T. J. Martínez and D. V. Shalashilin, J. Chem. Phys., 2014, 141, 054110.
43. D. V. Makhov, C. Symonds, S. Fernandez-Alberti and D. V. Shalashilin, Chem. Phys., 2017, 493, 200.
44. K. Saita and D. V. Shalashilin, J. Chem. Phys., 2012, 137, 22A506.
45. G. M. Roberts, J. L. Nixon, J. Lecointre, E. Wrede and J. R. R. Verlet, Rev. Sci. Instrum., 2009, 80, 053104.
46. P. M. Regan, S. R. Langford, A. J. Orr-Ewing and M. N. R. Ashfold, J. Chem. Phys., 1999, 110, 281.
47. B. G. Levine, J. D. Coe, A. M. Virshup and T. J. Martínez, Chem. Phys., 2008, 347, 3.
48. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, WIREs Comput. Mol. Sci., 2012, 2, 242.
49. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, et al., Molpro, version 2010.1, a package of ab initio programs, 2010, see http://www.molpro.net.
50. T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007.
51. A. L. Thompson, C. Punwong and T. J. Martínez, Chem. Phys., 2010, 370, 70.
52. N. d. N. Rodrigues, N. C. Cole-Filipiak, K. N. Blodgett, C. Abeysekera, T. S. Zwier and V. G. Stavros, Nat. Commun., 2018, 9, 5188.
53. D. V. Makhov and D. V. Shalashilin, Chem. Phys., 2018, 515, 46.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c8cp06359a

---