Exploring Ultrafast Dynamics of Pyrazine by Time ... - ACS Publications

Gaia Tomasello†, Alexander Humeniuk†, and Roland Mitrić‡. † Department of Physics, Freie Universität Berlin, Arnimallee 14, 14195 Berlin, Ge...
0 downloads 0 Views 6MB Size
Article pubs.acs.org/JPCA

Exploring Ultrafast Dynamics of Pyrazine by Time-Resolved Photoelectron Imaging Gaia Tomasello,† Alexander Humeniuk,† and Roland Mitrić*,‡ †

Department of Physics, Freie Universität Berlin, Arnimallee 14, 14195 Berlin, Germany Institut für Physikalische und Theoretische Chemie, Julius-Maximilians Universität Würzburg, Emil-Fischer-Straße 42, 97074 Würzburg, Germany



ABSTRACT: We present the simulation of time-resolved photoelectron imaging spectra of pyrazine in the gas phase. The approach we have adopted is based on the combination of the ab initio nonadiabatic molecular dynamics ”on the fly” with an approximate treatment of the photoionization process using Dyson orbitals and Coulomb functions to describe the bound and ionized states of the photoelectron. The method has been implemented (Humeniuk, A.; et al. J. Chem. Phys 2013, 139, 134104) in the framework of the time-dependent density functional theory and has been applied here to interrogate the ultrafast internal conversion between the S2 and S1 states in pyrazine. Conventional time-resolved photoelectron spectra without angular resolution fail to locate the S2 → S1 internal conversion, because the ionization potentials relevant for the photoionization channels S2 → D1 (π−1) and S1 → D0 (n−1) are almost identical. Introducing the angular resolution in the photoelectron spectra reveals evidence of such ultrafast internal conversion and provides a more detailed picture of the overall dynamics. The simulated time- and energy-dependent anisotropy map obtained within the Dyson/time-dependent density functional theory approach is in good agreement with its experimental counterpart provided by Horio et al. (Horio, T.; et al. J. Am. Chem. Soc. 2009, 131, 10932). Our theoretical approach represents a general tool for mapping the time- and angle-resolved photoelectron spectra in complex systems and thus can be used to investigate the ultrafast relaxation processes occurring in isolated molecules.



photodynamics.5,6 However, the comparison with theory is indispensable for the interpretation of complex experimental TRPEI spectra and for gaining a deeper insight into the nonadiabatic dynamics of polyatomic molecules. Photoionization is a complicated scattering process especially as the probe pulse becomes shorter and energy-richer and a lot of theoretical work has gone into determining which approximations are good enough to reproduce the main signatures of radiationless transitions as probed by detecting photoelectrons.7 Different methodologies have been developed to track nuclear dynamics through conical intersections by calculating time-resolved PADs. If a single reaction coordinate can be identified, the ionization potentials and photoangular distributions can be mapped along this coordinate using high-quality quantum chemistry methods and scattering calculations. This was done by the authors of ref 8 to study the isomerization of acetylene to vinylidene via a conical intersection and follow the evolution of the PADs in the molecular frame. For very small molecules (essentially only diatomics and triatomics) fully quantum-mechanical nuclear dynamics on coupled diabatic states can be performed as done by the authors of ref 9 to

INTRODUCTION One of the most powerful experimental techniques for studying nonadiabatic phenomena is the time-resolved photoelectron spectroscopy (TRPES) due to its intrinsic sensitivity both to the character of the electronic wave function and to the nuclear dynamics.1−3 Briefly, a TRPES experiment proceeds as follows: A pump laser pulse creates a wave packet on an excited state, which is probed after a variable time delay by a second pulse that leads to ionization. The intensity of the ejected electrons is measured as a function of their kinetic energy, providing a “picture” of how the system evolves.4 The time-resolved photoelectron imaging (TRPEI) technique1 is even more powerful since it provides, in addition to photoelectron energies, also information about their angular distribution, which can be directly linked to the character of the electronic state. The photoelectron angular distribution (PAD) therefore represents a more direct probe of the initial and target electronic wave functions and is able to provide a wealth of information about the nature of the electronic state that is ionized as well as about the molecular geometry and orientation as the electron departs. The photoelectron anisotropy parameters β2 are the basic quantities extracted from the PAD signal and depending on the orbital symmetry assume values between −0.5 and +2. The parameter β2 captures the information on the nature of the excited states involved in the © 2014 American Chemical Society

Special Issue: A. W. Castleman, Jr. Festschrift Received: February 18, 2014 Revised: June 13, 2014 Published: June 13, 2014 8437

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

⃗ wave functions, where the integration is scattering (|k⟩) performed over N − 1 electrons according to

investigate the signature in the TRPEI spectra of a wave packet passing through a conical intersection in the NO2 molecule. For larger molecules these approaches become less applicable as the TRPEI spectrum contains signals from different decay or reaction channels which might overlap and are difficult to disentangle. Computing photoangular distributions at a high level of theory for each potentially important nuclear configuration becomes prohibitively time-consuming. Our approach combines “on the fly” trajectory-based molecular dynamics10−12 accounting for nonadiabatic effects in the framework of the surface hopping method13,14 with an approximate treatment of the photoionization process using Dyson orbitals as initial and Coulomb waves as the final electron states.15−17 For a detailed description of the approach we refer the reader to our previous publication.24 We believe that the shortcomings of using fast time-dependent density functional theory vs wave function based methods to describe the electronic structure and our simplistic description of photoionization vs an exact solution of the scattering problem are offset by the ability to directly extract complete TRPEI spectra from a nonadiabatic dynamics simulation in a reasonable time, as if it were an experiment performed in the computer. In this work, we apply our methodology to simulate the TRPEI spectra of pyrazine, which has been used as a model for the investigation of the internal conversion (IC) through a conical intersection in a number of experimental18,19 and theoretical studies.20,21 Previous experimental studies have demonstrated that the photoelectron kinetic energy (PKE) distribution does not significantly change upon the S2[1B2u (π,π*)] → S1[1B3u (n,π*)] internal conversion since the energy gaps for the dominant S2(π,π*) → D1(π−1) and S1(n,π*) → D0(n−1) photoionization channels are almost the same.5 However, the time dependence of the experimental photoelectron angular anisotropy parameter β2 observed for photoionization of pyrazine after photoexcitation to the S2 state, reported by Suzuki et al.,5,22 clearly exhibited an abrupt change at around 20 fs reflecting the change of the electronic character from π,π* to n,π* upon S2 → S1 internal conversion through a conical intersection. These studies have thus demonstrated that angular distributions of photoelectrons reveal the underlying nonadiabatic dynamics quite clearly.

N

⟨k |⃗ ⊗ ⟨Φ N − 1|∑ E ⃗ · ri |⃗ Φ⟩ = ⟨k |⃗ E ⃗ · r |⃗ ΦD⟩ i=1

Since the integration only extends over the coordinates of N − 1 electrons, the Dyson orbital still has a dependence on one electron coordinate, thus taking into account the electronic reconfiguration between the neutral molecule and the ion. Notice that since we calculate the Dyson orbitals using relaxed cationic states (calculated at the linear response−time-dependent density functional theory (LR-TDDFT) level), the norm of the Dyson orbital can be smaller than one. The Dyson orbital will look like the highest occupied molecular orbital (HOMO) or lowest unoccupied molecular orbital (LUMO) and will have a norm equal to 1, only if the electronic structure of the cation contains the same molecular orbitals as the neutral molecule. In general this is not the case and the norm of the Dyson orbital for the ionization from a neutral to a cationic state can be taken as a measure of the importance of this channel which is automatically included in our approach. Since a full description of the photoionization process requires the solution of a scattering problem, which is computationally extremely demanding for polyatomic molecules, we have assumed that the continuum states of the photoelectron are approximately described by Coulomb waves, which is reasonable given the cationic character of the remaining molecule after ionization. In our previous study24 we have shown by comparison with experimental data that the Dyson orbitals obtained from the TDDFT can be used to capture the leading features of the experimental photoelectron angular distributions. However, it should be noted that the form of Dyson orbitals depends on the density functional used due to the self-interaction error, and range-separated functionals should generally provide better performance. As noted in ref 16, the description of the continuum electron can be improved considerably by orthogonalizing the Coulomb waves with respect to the Dyson orbital. This effectively imprints knowledge of the molecular potential onto the continuum wave function. This is particularly important at low photokinetic energies where the scattering wave function of the leaving electron is strongly perturbed by the molecular potential of the remaining cation. In order to test this approximation, we have in our previous studies on furan24 numerically investigated the description of scattering states at selected geometries by expanding the electrostatic potential of the cation into spherical harmonics and integrating the coupled-channel equations up to high angular momentum orders. Although, in detail, the scattering states looked different from Coulomb waves in the region around the molecule, we did not observe a qualitative change in the cross-section and the anisotropy parameters. We point out that the present approach is computationally much more efficient and in contrast to coupled-channel simulations can be used for a large number of geometries along the nonadiabatic trajectories. The Coulomb wave could be orthogonalized by the Gram− Schmidt procedure:



METHODS Theoretical Description of the Photoionization Process. The process of photoionization from a polyatomic molecule is simulated by a model where the initial and final states are multielectron states consisting initially of an electronically excited molecule with a bound electron that interacts with a photon to produce an ion and a free electron. By noting that the molecule interacts with light through coupling to the dipole operator, which is a single-particle operator, the problem can be reduced to a single-particle process. Within this approximation the Dyson orbital1,17 arises as a single-particle wave function defined as the overlap between the neutral state with N electrons and the cation state with N − 1 electrons: |φD⟩ =

N ⟨Φ N − 1|Φ N ⟩

(2)

|k⟩̃ ∝ |k⟩ − ⟨ψ D|k⟩|ψ D⟩

(1)

The probability of an electron being ejected in a certain direction (generating the PAD signal) is given by the transition dipole integral, which connects the initial (|ϕD⟩) and final

(3)

However, since we are only interested in the matrix element of the dipole operator between the Dyson orbital and the continuum state, we shift the Dyson orbital in space, so that the 8438

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

of cos(θk) but these are usually small. To be able to use analytical expressions for σ and β2,23 the dipole interaction with a polarized electromagnetic field E⃗ mp (mp = 0, −1, 1 stand for linearly, left, and right polarized probe pulse, respectively) and the initial and final states ψD and ψ(−) are expanded into spherical harmonics with respect to r ⃗ and k.⃗ Then the formulas (13) and (14) from ref 24 are used to compute σ and β2 from the partial waves of the Dyson orbital and the continuum orbital. For details about the orientation averaging, see refs 15, 16, 23, and 24. Because Dyson orbitals depend parametrically on a nonadiabatic nuclear trajectory R⃗ (t) and since the ionization takes place on a much shorter time scale than the movement of the nuclei, we calculate the ground and excited states of the neutral and ionized molecule at the same structure, corresponding to extracted snapshot geometries along the dynamics. For the calculation of Dyson orbitals along the trajectories we employ the procedure previously introduced in ref 24. Computational Methods. The geometry of pyrazine has been optimized on the ground state using density functional theory employing the B3LYP25,26 functional and the TZVP basis set27 as implemented in the TURBOMOLE28 package. The nonadiabatic dynamics simulations have been carried out in the frame of the time-dependent density functional theory29 combined with Tully’s surface hopping method.13,14,30,31 For the TDDFT calculations, the same exchange-correlation functional and basis set as for the geometry optimization have been used. The 100 initial coordinates and momenta have been generated by propagating a trajectory in the electronic ground state at constant temperature (298 K) for 10 ps and then sampling it each 100 fs. For the propagation of the nuclei, the Newtonian equations of motion have been integrated using the Verlet algorithm32 with a time step of 0.1 fs. Along each trajectory, the electronic Schrödinger equation has been integrated using the fourth-order Runge−Kutta method with a time step of 10−5 fs. In our simulation, we have performed nonadiabatic surface hopping molecular dynamics in the manifold of the ground state and the lowest five excited singlet states, starting directly from the bright state S2. We have propagated the trajectories for 200 fs with a time step of 0.1 fs. Dyson orbitals and the anisotropy parameters have been determined by calculating the electronic structure of the neutral molecule (N electrons) and the ionized molecule (N − 1 electrons) for each geometry extracted at intervals of 1 fs along each trajectory. We have considered the photoionization processes involving the current state Sn and the two lowest cation states D1 and D0, which are the primary processes in agreement with the experimental data.22 In order to calculate orientation-averaged transition dipole matrix elements, the Dyson orbitals were shifted by their dipole moments and expanded into spherical harmonics on a Lebedev grid.33 Spherical wave components for Coulomb waves were computed on the same radial grid for a range of PKE values. For each trajectory i, maps for the total cross-section σi(PKE,t) and the anisotropy parameters β2,i(PKE,t) have been calculated using the Dyson orbitals at time t as initial states and Coulomb waves for an electron with momentum k = (2PKE)1/2 as final states. In order to calculate the anisotropy parameter along one nonadiabatic trajectory, we obtained the transition dipole moment d⃗0→I = (dx,dy,dz) between the ground state and the initial excited state I = 3 (S2) in the first time step, which coincides with the time of the pump pulse. The photoionization

expectation value of the dipole operator in the Dyson orbital becomes zero. The matrix element will then be the same as for the Gram−Schmidt orthogonalized states:16

Moreover it has to be noticed that the electromagnetic field E⃗ and electron momentum k⃗ are vectors defined in the laboratory frame, while the Dyson orbital and the position vector are defined in the molecular frame. Thus, one frame has to be rotated into the other by a rotation matrix R in order to calculate the matrix element in eq 2: ⟨Ψ(R(−)α , β , γ )k ⃗|(R(α , β , γ )E ⃗) · r |⃗ ΨD⟩

(5)

α, β, and γ are the Euler angles which parametrize the relative orientation of the molecular to the laboratory frame. The continuum orbital |ψ(−) k⃗ ⟩ is a scattering state, which asymptotically resembles a free electron with linear momentum in the direction of k̂ in the laboratory frame. The detector is positioned so as to detect electrons ejected into this direction. |ψ(−) Rk⃗ ⟩ represents the same state as that seen from the molecular frame, where the vector k,̂ which points toward the detector, now depends on the relative orientation of the molecule. Among the isotropically oriented molecules in the ground state, a linearly polarized pump pulse preferentially excites those ⃗ (where I is the molecules whose transition dipole moment d0→I initially pumped state) is aligned with the polarization vector E⃗ . The transition probability T0→I depends therefore on the angle Θ = ∠(E⃗ ,d⃗0→I) between the two vectors: T0 → I ∝ ⟨Ψ|I E ⃗ · r |⃗ Ψ0⟩||2 ∝ cos2(∠(E ⃗ , d0⃗ → I))

(6)

This leads to an excited state population with aligned dipole moments. If the time delays between pump and probe pulse are shorter than the rotational period, the molecules cannot reorient. The probe pulse therefore encounters an aligned ensemble in the excited state that can be characterized by a distribution of orientations f(R) which is a function of the Θ: faligned (R) =

3 2 cos (β(R) − Θ) 2

(7)

where β(R) − Θ indicates how much the Euler angle of the molecular orientation matrix R deviates from the perfectly aligned transition dipole moment. In order to calculate the photoionization cross-section, averaging over molecular orientations weighted by the distribution f(R) needs to be performed: dσ ∝ dΩ

∫ dR |⟨Ψ(R−k)⃗ |E⃗·(R−1r ⃗)|ΨD⟩|2 f (R)

(8)

The averaged scattering cross-section depends only on the angle between the polarization vector (or the z-axis for light circularly polarized in the xy-plane) and the electron momentum and can be characterized for nonchiral molecules by two quantities, the total scattering cross-section σ and the anisotropy parameter β2 (ranging between −0.5 and 2.0): dσ σ = [1 + β2P2(cos(θ ))] dΩ 4π

(9) 6

where P2(x) is the second Legendre polynomial. A TRPEI experiment is fully described by maps of these two parameters in the time and PKE domain. For prealigned molecules, the averaged cross-section can also depend on even higher powers 8439

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

Figure 1. Electronic absorption spectrum of pyrazine: (a) theoretical and (b) experimental. The vertical excitations have been calculated at the B3LYP/TZVP level on 100 structures from a thermal ensemble. The individual transitions (black sticks) have been convoluted with a Lorentzian.

Figure 2. Potential energy diagram for electronic and ionization transition energies calculated (a) on the optimized structure on the ground state (B3LYP/TZVP). Panel b shows the character of the relevant electronic excited states.



cross-section was then averaged over molecular orientations with the (3/2)cos 2 (β − Θ) weight factor, where Θ = arccos(dz/|d⃗|). Finally the maps of all trajectories were combined. Each trajectory produces one curve in (PKE,t) space, which is defined by energy conservation hν = PKE + IE(t) (ν is the frequency of the probe pulse and IE(t) the energy difference between the current neutral state Sn and the cation states D1 and D0 along the trajectory). After running sufficiently many trajectories with different initial conditions, the relevant part of the (PKE,t) space is filled with data points. The contributions of different trajectories i to the same point in the anisotropy map are averaged using the following formula:

The experimental absorption spectrum of pyrazine in the gas phase is characterized mainly by two bands corresponding to the two lowest bright singlet states S1 and S2 at the energies ∼ 3.8 and ∼4.8 eV. The simulated absorption spectrum, calculated at the B3LYP/TZVP level on a thermal ensemble, is able to reproduce the general experimental shape (see Figure 1). At the B3LYP/TZVP level we have located a first excited state S1(1B3u) characterized by an n → π* transition and a minor oscillator strength. The next higher excited state 1Au is dark and also of n → π* character. The excited state S2(1B2u) features a π → π* (HOMO−1 → LUMO) transition and is the bright state. Notice that the state numbering was performed according to the experimental convention. Between the bright

β(PKE, t ) = (∑ σi(PKE, t ) βi (PKE, t ))/(∑ σi(PKE, t )) i

RESULTS. TRPEI ON PYRAZINE

i

8440

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

simulation. In both trajectories the ultrafast S2 → S1 internal conversion is completed within ∼30 fs. Then the dynamics proceeds on the S1 potential energy surface (trajectories of type I) and just in a few cases (5%, trajectories of type II) the kinetic energy re-distribution leads to the excitation of the ringpuckering mode, responsible for the decay to the ground state via the low-lying S1/S0 conical intersection. Thus, while the trajectories of the type II decay to the ground state, the type I trajectories remain in the excited state within the simulation time. To provide a detailed explanation of the nonadiabatic dynamics of pyrazine and the consequent ionization process, we analyze first the selected nonadiabatic trajectory of the type I presented in Figure 4. The dominant electronic configuration

states denoted by Sn several dark singlet states are additionally present (cf. Figure 2). Higher in energy is the third optically active state S5(1B1u), which has π → π* character consisting of a HOMO−1 → LUMO+1 electronic excitation. The order and the character of the states is summarized in Figure 2, where we show the vertical excitations calculated on the ground state optimized structure. The calculated ionization potentials S0 → D0(n−1) (9.41 eV) and S0 → D1(π−1) (10.21 eV) are also in agreement with the experimental data (9.22 and 10.18 eV).34 In order to simulate time-resolved photoelectron imaging spectra of pyrazine, we first performed the nonadiabatic dynamics propagating 100 trajectories for 200 fs. The timedependent populations on the excited states have been determined by monitoring the number of trajectories in each state. Starting from the bright state S2(1B2u), the population first relaxes to the lower dark state 1Au and then to S1(1B3u). The population on the initial state decays with a time constant of 41 fs (see Figure 3), in quite good agreement with the

Figure 3. Time-dependent excited states populations obtained by averaging 100 nonadiabatic trajectories. Nonadiabatic transitions between the lowest five excited singlet states were included. Gradients and surface hopping probabilities were computed at the TDDFT/ B3LYP level of theory.

experimental measurements18,34 and previous calculations.20,21 In accord with previous findings, the ground state is not populated significantly in the first 200 fs. Indeed, from the entire ensemble of the propagated trajectories, only 5% decay to S0 through a nonradiative mechanism involving a ringpuckering coordinate. The time-resolved photoelectron spectra have been calculated including the photoionization processes to both the cationic states D0 and D1. The calculated Dyson-orbital norm and total cross-section reveal which photoionization processes occur during the nonadiabatic dynamics. According to Koopmans’ theorem and in agreement with the experiments,5,22 within the first 40−50 fs, the S2(π,π*) → D1(π−1) photonization process contributes primarily, being replaced by the S1(n,π*) → D0(n−1) photonization process for later times, which becomes progressively dominant until the end of simulation. In order to illustrate the relaxation mechanism, we have selected two examples of nonadiabatic trajectories, representing the typical reaction mechanisms observed during the

Figure 4. Trajectory of type I. (a) Relative energies of the lowest four neutral electronic states and the lowest two ionized states (D1 and D0) along a selected nonadiabatic trajectory. The actual state in which the trajectory resides is labeled by the orange dashed line. Dyson orbitals corresponding to the photoionization processes S2 → D1(π−1) and S1 → D0(n−1) are plotted each 60 fs with the associated β2 values for a probe pulse energy of 6.26 eV. (b) Character of the relevant orbitals in the wave function of the current state along the trajectory.

at the beginning of the dynamics, at t = 0 fs, corresponds to the bright state S2. As can be seen from the HOMO and LUMO orbitals dominantly involved in the transition (cf. panel b in Figure 4), the character of the state changes from π → π* to n → π* during the dynamics within 60 fs. The time-dependent energies of the four excited states of the trajectory in Figure 4a are well separated from the ground state within the simulation period, which is in agreement with the absence of population transfer to the ground state. As the time proceeds, the trajectory remains on the S1 potential energy surface and the n → π* character is conserved until the end of the simulation. 8441

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

orbital (see Figure 4) and do not change significantly during the propagation time. In Figure 6a we present first the simulated β2 parameters for the Franck−Condon geometry for both photoionization

Although only a small number of nonadiabatic trajectories belongs to the type II, which decays to the ground state, it is worthwhile to briefly comment on a representative trajectory of this type shown in Figure 5. Here, although at t = 0 fs the bright

Figure 6. Calculated anisotropy coefficients β2 at the same Franck− Condon geometry for both photoionizations S2 → D1(π−1) and S2 → D0(n−1).

channels S2 → D1(π−1) and S2 → D0(n−1) as a function of the photoelectron kinetic energy. As can be seen, the PKE dependence for both ionization channels is very similar, decreasing from an initial value of 0.6 to 0.0 over the PKE range between 0.0 and 3.0 eV (cf. Figure 6). A similar trend has been previously noticed in the anisotropies calculated by Suzuki and Suzuki using the FOCI/CMSXα method35 for the same photoionization processes. During the course of the nonadiabatics dynamics, most trajectories relax from a π → π* to an n → π* state but conserve the same LUMO character. Therefore, the Dyson orbitals and their related anisotropies do not change substantially. The anisotropy increases only because the ejected electron has less kinetic energy at its disposal, so that the anisotropy needs to be evaluated at a lower PKE, where its value is larger. In Figure 7 we analyze the anisotropy rise in a selected trajectory by comparing its values at t = 0 fs and t = 100 fs for the probe pulse energies of 6.26 and 7.84 eV. For both photoionization channels, the anisotropy parameters evaluated at different time steps along the same trajectory increase as the PKE values decrease. At t = 0 fs the photoionization occurs mainly from the S2 bright state to the first cationic state D1, while at t = 100 fs it occurs from the S1 to the cationic ground state D0. Note that, whereas the transitions S2(ππ*) → D1(π−1) and D0(nπ*) → D0(n−1) are single-electron processes, ionization from S2(ππ*) to D0(n−1) is not, since the configuration of D0(π−1) cannot be obtained from S2(ππ*) by simply removing an electron. Therefore, the norm of the Dyson orbital for the transition S2(ππ*) to D0(n−1) is very low under the assumption that the laser intensity is low enough to exclude multiphoton processes. The anisotropy parameter β2, being a function of PKE, reflects of course also the dependence on the energy of the probe pulse and is shifted to lower values as the energy of the probe pulse increases. Therefore, the overall anisotropy increase can be attributed to the rise of the ionization potentials as the dynamics proceeds. The comparison between experimental and simulated timeand energy-dependent anisotropy maps for the first 200 fs is depicted in Figure 8. We compare the theoretical anisotropy

Figure 5. Trajectory of type II. (a) Relative energies of the lowest four neutral electronic states and the lowest two ionized states (D1 and D0) along a selected nonadiabatic trajectory decaying on the ground state. The actual state in which the trajectory resides is labeled by the orange dashed line. Dyson orbitals corresponding to the photoionization processes S2 → D1(π−1) and S1 → D0(n−1) are plotted at t = 0 fs, t = 80 fs, and t = 112 fs with the associated anisotropy values for a probe pulse energy of 6.26 eV. (b) Character of the relevant orbitals in the wave function of the current state along the trajectory.

state is again of π → π* nature, the internal conversion leads to an energetically accessible deactivation channel to the electronic ground state. The radiationless decay to the ground state is indeed observed at 112 fs, due to a low-lying S1/S0 conical intersection mediated by a ring-puckering coordinate. In the following, we focus the attention on the PAD calculated by averaging over the whole propagated ensemble. As already discussed, both the D1−D0 and S2−S1 energy gaps amount to about 0.90 eV. This leads to an almost unchanged PKE distribution upon S2 → S1 internal conversion, since the ionization energies for the S2 → D1(π−1) and S1 → D0(n−1) processes are almost identical. In addition, due to the electronic character of the states S2(π→π*) and S1(n→π*) involved in the internal conversion, the electron is in both cases ejected from the same π* orbital. We therefore expect to obtain similar values of the anisotropy parameter β2 for both channels. This is also supported by the appearance and symmetry of the Dyson orbitals, which mainly reflect the symmetry of the LUMO π* 8442

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

Figure 7. Anisotropy parameters β2 (a) after 0 fs and (b) after 100 fs in a typical trajectory. At t = 0 fs, the dominant process is the photoionization S2 → D1(π−1), while, at t = 100 fs, the photoionization S1 → D0(n−1) has the main contribution. One sees how the anisotropy increases when the PKE decreases and therefore depends on the different probe pulse energies employed.

maps obtained using both the probe pulse energies of 6.26 eV (Figure 8b) and 7.84 eV (Figure 8c) with the experimental map obtained for pump/probe pulse energies of 4.12/6.26 eV, respectively (Figure 8a). A rapid change in the anisotropy values is clearly visible on the time scale of 200 fs. For the 6.26 eV probe pulse the red/green/blue colors correspond to the ejection of an electron from an orbital with anisotropy parameters in the range of 0.0−0.8. These values are typical for an electron ejected from a Dyson orbital with π → π* symmetry, which characterizes the entire S2 → S1 dynamics in pyrazine. The simulated anisotropy maps are in good general agreement with the experiment and reproduce well the relative change in the anisotropy parameters. Notice that while the qualitative features of the signal are well reproduced by the simulation, the PKE range is too broad. This is due to the inaccuracies of TDDFT in describing quantitatively the excited state ionization energies. Moreover, the experimental map exhibits also a narrow region of negative anisotropy parameters. The dark state 11Au(nπ*) is transiently populated in the first 50, which might explain the negative anisotropy in this time window. However, our simulated anisotropy map does not reproduce this feature because the character of the Dyson orbital does not change, neither in the S2(ππ*) → 11Au(nπ*) transition nor in the 11Au(nπ*) → S1(nπ*) transition. In a previous study the first-order configuration interaction method combined with the more accurate multiple scattering Xα approximation has been used to elucidate the origin of this region35 showing that it arises due to the presence of a kag shape resonance within the ionization continuum. The description of the shape resonances requires the full solution of the scattering problem along the nonadiabatic trajectories

Figure 8. (a) Experimental5 time- and energy-resolved anisotropy map obtained with pump/probe pulse energies of 4.70 eV/6.26 eV and a time resolution of 22 fs. (b and c) Simulated anisotropy map β2(PKE,t) from averaging the entire ensemble of 100 nonadiabatic trajectories. The anisotropy parameters have been calculated every 1 fs of the 200 fs long simulation. In panel a the experimental probe pulse energy of 6.26 eV was used, while in panel c the probe pulse energy was increased to 7.84 eV.

which is still computationally out of the reach. The general trend of the lowering of the anisotropy parameter with PKE is also evident if a higher probe pulse energy of 7.84 eV is employed in the simulation. The simulated anisotropy map obtained for a 7.84 eV probe pulse (Figure 8c) leads to a higher PKE and shifts up the bands in the anisotropy map. For this probe energy the red region associated with negative values of β2 is more pronounced and extends over a larger PKE and time interval.



CONCLUSION Theoretical nonadiabatic dynamics simulations combined with experimental time-resolved photoelectron imaging spectroscopy represent a general and most direct approach for the interrogation of coupled nuclear-electron dynamics and thus for the investigation of photochemistry and photophysics in polyatomic molecules. In order to establish the direct link between the theory and experiment, it is highly desirable to develop an efficient theoretical methodology for the simulation 8443

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

(8) Douguet, N.; Rescigno, T.; Orel, A. Time-Resolved MolecularFrame Photoelectron Angular Distributions: Snapshots of AcetyleneVinylidene Cationic Isomerization. Phys. Rev. A 2012, 86, No. 013425. (9) Arasaki, Y.; Takatsuka, K.; Wang, K.; McKoy, V. Time-Resolved Photoelectron Spectroscopy of Wavepackets through a Conical Intersection in NO2. J. Chem. Phys. 2010, 132, No. 124307. (10) Barbatti, M.; Granucci, G.; Persico, M.; Ruckenbauer, M.; Vazdar, M.; Eckert-Maksić, M.; Lischka, H. The On-the-Fly SurfaceHopping Program System NEWTON-X: Application to ab Initio Simulation of the Nonadiabatic Photodynamics of Benchmark Systems. J. Photochem. Photobiol. A:Chem. 2007, 190, 228−240. (11) Mitrić, R.; Bonacčić-Koutecký, V.; Pittner, J.; Lischka, H. Ab Initio Nonadiabatic Dynamics Study of Ultrafast Radiationless Decay over Conical Intersections Illustrated on the Na3F Cluster. J. Chem. Phys. 2006, 125, No. 024303. (12) Lischka, H.; Dallos, M.; Shepard, R. Analytic MRCI Gradient for Excited States: Formalism and Application to the n-pi* Valence- and n-(3s,3p) Rydberg States of Formaldehyde. Mol. Phys. 2002, 100, 1647−1658. (13) Tully, J. Molecular Dynamics with Electronic Transitions. J. Chem. Phys. 1990, 93, 1061−1071. (14) Mitrić, R.; Petersen, J.; Wohlgemuth, M.; Werner, U.; BonacčićKoutecký, V.; Wöste, L.; Jortner, J. Time-Resolved Femtosecond Photoelectron Spectroscopy by Field-Induced Surface Hopping. J. Phys. Chem. A 2011, 115, 3755−765. (15) Oana, C.; Krylov, A. Dyson Orbitals for Ionization from the Ground and Electronically Excited States within Equation-of-Motion Coupled-Cluster Formalism: Theory, Implementation and Examples. J. Chem. Phys. 2007, 127, No. 234106. (16) Oana, C.; Krylov, A. Cross Sections and Photoelectron Angular Distributions in Photodetachment from Negative Ions using Equationof-Motion Coupled-Custer Dyson Orbitals. J. Chem. Phys. 2009, 131, No. 124114. (17) Walter, M.; Haekkinen, H. Photoelectron Spectra from First Principles: from the Many-Body to the Single-Particle Picture. New J. Phys. 2008, 10, No. 043018. (18) Suzuki, Y. I.; Fuji, T.; Horio, T.; Suzuki, T. Time-Resolved Photoelectron Imaging of Ultrafast S2 → S1 Internal Conversion through Conical Intersection in Pyrazine. J. Chem. Phys. 2010, 132, No. 174302. (19) Oku, M.; Hou, Y.; Xing, X.; Reed, B.; Xu, H.; Chang, C.; Ng, C. Y.; Nishizawa, K.; Ohshimo, K.; Suzuki, T. 3s Rydberg and Cationic States of Pyrazine Studied by Photoelectron Spectroscopy. J. Phys. Chem. A 2008, 112, 2293−2310. (20) Werner, U.; Mitrić, R.; Suzuki, T.; Bonačić-Koutecký, V. Nonadiabatic Dynamics within the Time Dependent Density Functional Theory: Ultrafast Photodynamics in Pyrazine. J. Chem. Phys. 2008, 349, 319−324. (21) Werner, U.; Mitrić, R.; Bonačić-Koutecký, V. Simulation of Time Resolved Photoelectron Spectra with Stieltjes Imaging Illustrated on Ultrafast Internal Conversion in Pyrazin. J. Chem. Phys. 2010, 132, 174301−174310. (22) Tang, Y.; Shen, H.; Sekiguchi, K.; Kurahashi, N.; Mizuno, T.; Suzuki, Y.; Suzuki, T. Direct Measurement of Vertical Binding Energy of a Hydrated Electron. Phys. Chem. Chem. Phys. 2010, 12, 3653−3655. (23) Ritchie, B. Theory of the Angular Distribution of Photoelectrons Ejected from Optically Active Molecules and Molecular Negative Ions. Phys. Rev. A 1976, 13, 1411−1415. (24) Humeniuk, A.; Wohlgemuth, M.; Suzuki, T.; Mitric, R. TimeResolved Photoelectron Imaging Spectra from Nonadiabatic Molecular Dynamics Simulations. J. Chem. Phys. 2013, 139, No. 134104. (25) Becke, A. D. Densityfunctional Thermochemistry. III. The Role of Exact Exchange. J. Chem. Phys. 1993, 98, 5648−5652. (26) Becke, A. D. Densityfunctional Thermochemistry. V. Systematic Optimization of Exchange-Correlation Functionals. J. Chem. Phys. 1997, 107, 8554−8560. (27) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully Optimized Contracted Gaussian Basis Sets of Triple Zeta Valence Quality for Atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829−5836.

of the time-resolved photoelectron anisotropy maps which are the observable extracted from the TRPEI experiments. While the theory can in principle provide a full picture of the dynamical processes, only by reproducing the direct experimental observables, confidence in the theoretical simulations can be established and further conclusions about the dynamical processes can be drawn. We have presented here a simulation of the ultrafast nonadiabatic dynamics and time-resolved photoelectron imaging spectra of pyrazine based on our recently introduced Dyson/TDDFT method. Our approach allows us to efficiently calculate the Dyson orbitals along the nonadiabatic surface hopping trajectories. The anisotropy maps are simulated using the Dyson orbitals for the bound electrons and Coulomb waves for the free electrons. We have demonstrated on the example of pyrazine that such an approach can reproduce the main features of the time-resolved anisotropy maps in regions where the influence of shape resonances is small. The agreement between the simulated anisotropy map and its experimental counterpart shows that the TDDFT nonadiabatic dynamics simulation correctly describes the photodynamics of pyrazine. We find that in pyrazine the character of the Dyson orbital does not change significantly as the wave packet switches nonradiatively from S2 to S1. Our theoretical approach should serve to help to analyze and interpret the experimental time- and angle-resolved photoelectron spectra in complex systems. The combination of the approximate theoretical methodology with advanced experimental techniques should serve to advance our understanding of the fundamental aspects of photochemistry and photophysics.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS G.T. acknowledges financial support from the Postdoc fellowship of the A. von Humboldt Foundation (Grant Card No. 7000260646).



REFERENCES

(1) Suzuki, T. Femtosecond Time-Resolved Photoelectron Imaging. Annu. Rev. Phys. Chem. 2006, 57, 555−592. (2) Stolow, A.; Underwood, J. Time-Resolved Photoelectron Spectroscopy of Nonadiabatic Dynamics im Polyatomic Molecules. Advances in Chemical Physics; John Wiley & Sons: New York, 2008; Vol. 139, chapter 6. (3) Suzuki, T. Time-Resolved Photoelectron Spectroscopy of NonAdiabatic Electronic Dynamics in Gas and Liquid Phases. Int. Rev. Phys. Chem. 2012, 31, 265−318. (4) Wu, G.; Hockett, P.; Stolow, A. Time-Resolved Photoelectron Spectroscopy: From Wavepackets to Observables. Phys. Chem. Chem. Phys. 2011, 13, 18447−18467. (5) Horio, T.; Takao, F.; Suzuki, Y.; Suzuki, T. Probing Ultrafast Internal Conversion through Conical Intersection via Time-Energy Map of Photoelectron Angular Anisotropy. J. Am. Chem. Soc. 2009, 131, 10392−10393. (6) Reid, K. Photoelectron Angular Distributions. Annu. Rev. Phys. Chem. 2003, 54, 397−424. (7) Seideman, T. Time-Resolved Photoelectron Angular Distributions as a Probe of Coupled Polyatomic Dynamics. Phys. Rev. A 2001, 64, No. 042504. 8444

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445

The Journal of Physical Chemistry A

Article

(28) Ahlrichs, R.; Bar, M.; Haser, M.; Horn, H.; Kolmel, C. Electronic Structure Calculations on Workstation Computers: The Program System Turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (29) Casida, M. Time-Dependent Density-Functional Response Theory for Molecules. In Recent Advances in Density Functional Methods; Chong, D., Ed.; World Scientific: Singapore, 1995. (30) Mitrić, R.; Werner, U.; Bonačić-Koutecký, V. Nonadiabatic Dynamics and Simulation of Time Resolved Photoelectron Spectra within TDDFT: Ultrafast Photoswitching in Benzylidenaniline. J. Chem. Phys. 2008, 129, No. 164118. (31) Bonačić-Koutecký, V.; Mitrić, R.; Buergel, C.; Nossler, M. Ultrafast Dynamics in Noble Metal Clusters: The Role of Internal Vibrational Redistribution. Chem. Phys. 2008, 350, 111−117. (32) Verlet, L. Computer Experiments on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. Phys. Rev. 1967, 159, 98−103. (33) Lebedev, V. Values of the Nodes and Weights of Quadrature Formulas of Gauss−Markov Type for a Sphere from the Ninth to Seventeenth Order of Accuracy that are Invariant with Respect to an Octahedron Group with Inversion. USSR Comput. Math. Math. Phys. 1975, 15, 48−54. (34) Stert, V.; Farmanara, P.; Radloff, W. Electron Configuration Changes in Excited Pyrazine Molecules Analyzed by Femtosecond Time-Resolved Photoelectron Spectroscopy. J. Chem. Phys. 2000, 112, 4460−4464. (35) Suzuki, Y.-I.; Suzuki, T. Effect of Electron Correlation and Shape Resonance on Photoionization from the S2 and S1 States of Pyrazine. J. Chem. Phys. 2012, 137, No. 194314.

8445

dx.doi.org/10.1021/jp5017328 | J. Phys. Chem. A 2014, 118, 8437−8445