Attosecond electronic and nuclear quantum photodynamics of ozone monitored with time and angle resolved photoelectron spectra

Recently we reported a series of numerical simulations proving that it is possible in principle to create an electronic wave packet and subsequent electronic motion in a neutral molecule photoexcited by a UV pump pulse within a few femtoseconds. We considered the ozone molecule: for this system the electronic wave packet leads to a dissociation process. In the present work, we investigate more specifically the time-resolved photoelectron angular distribution of the ozone molecule that provides a much more detailed description of the evolution of the electronic wave packet. We thus show that this experimental technique should be able to give access to observing in real time the creation of an electronic wave packet in a neutral molecule and its impact on a chemical process.


Theoretical Background
A molecule such as ozone can be viewed as a collection of N nuclei and n electrons.
denote the position vectors of the nuclei and the electrons, respectively. Using a semi-classical approach with respect to the external electromagnetic field and the so-called dipole approximation, the non-relativistic Coulomb molecular Hamiltonian operator for the system interacting with a time-dependent external electric field, E t ( ), reads with Ψ r R t ( , , ) the wave packet of the molecule. The adiabatic electronic basis functions, Φ r R ( ; ) i , satisfy for each R Scientific RepoRts | 6:36613 | DOI: 10.1038/srep36613 where R are to be viewed as parameters and E R ( ) i el play the role of potential energy surfaces for the nuclei. Here, we consider only a pair of adiabatic electronic states for ozone: X( 1 A 1 ), the ground state, and B( 1 B 2 ), the Hartley excited state. The total wave function of the molecule can be expanded as In the following, we assume the Born-Oppenheimer approximation to be valid and thus neglect the non-adiabatic couplings between the two electronic states stemming from the nuclear kinetic energy operator. The only coupling between X and B is induced by the external field through the term µ , where the transition dipole is defined as . We also neglect the diagonal terms involving µ R ( ) XX and µ R ( ) BB since E t ( ) is an external field resonant between X and B with respect to the central wavelength of the spectrum of the pulse.
Thus, the evolution of Ψ R t ( , ) is governed by a set of two coupled equations involving only To solve this set of equations, i.e. to solve the Schrödinger equation for the nuclei, we use the MCTDH method [16][17][18][19][20][21]29 . The nuclear wave functions are expanded in a basis set of timedependent functions, the so-called single-particle functions (SPFs), where f denotes the number of nuclear degrees of freedom (Q κ are single coordinates or groups of coordinates involved in R). There are n κ SPFs for the κth nuclear degree of freedom. The equations of motion [16][17][18][19][20][21] for the A-coefficients and the SPFs are derived from a variational principle that ensures optimal convergence.
The parameters defining E t ( ), the laser pump pulse (see Fig. 1) are: central wavelength at 260 nm, intensity of 10 13 W/cm 2 , Gaussian envelope with a full duration at half maximum (FDHM) equal to 3 fs. Note that, due to the C 2v symmetry of the ozone molecule at the Franck-Condon (FC) point (R 1 = R 2 = 1.275 Å; α = 116.9°), the y-component (B 2 ) of the transition dipole between X and B is the only non-vanisihing one at the FC point and is thus primarily responsible for the light-induced electronic transitions. Consequently, the effective polarization axis of the electric field is y.
Further details regarding our calculations -the (time-independent) primitive basis sets, the parameters for the complex absorbing potentials, the refitting of the potential energy and transition dipole surfaces in a form adapted to MCTDH, and the number of SPFs -can be found in previous work, for instance in Sec. 3 of ref. 31.
Starting from the vibrational ground state in the electronic ground state X, MCTDH calculations will generate Ψ R t ( , ) X and Ψ R t ( , ) B at any subsequent time. Assuming that only the B electronic state is populated by the laser pulse (see Fig. 1), the total molecular wave packet (see Eq. 4) can be constructed provided the corresponding adiabatic electronic wave functions are known.
Thus, with this approach, we can obtain in principle the full electronic and vibrational wave packet (note again that we only consider the case where the total angular momentum is equal to 0). However, this quantity cannot  be observed directly in actual experiments and we need a time-resolved property that will characterize the time evolution of the system: the TRPES for instance, which can be measured and compared to calculations. The procedure that we used to compute this quantity is explained below. As a first approximation, we can consider that the early stages of the process will be dominated by the behavior of the wave packet at the FC point, R FC . The corresponding renormalized density matrix of the molecule at the FC point (see Sec. II B of ref. 9 for further details) reads, for i, i′ = X, B, Note that such local populations of X and B are not classical quantities but extracted from the actual quantum wave packets. Assuming a "stationary" picture, the approximate photoelectron spectra from either X 35 or B at the FC point appear as stick spectra, where  is the kinetic energy (KE) of the ejected electron, i = X or B, and k is used to label the various cation states.  ik are the corresponding peaks appearing in the spectra. They satisfy where E photon denotes the energy of the probe photon, 95 eV here. E i are the energies of the X and B states at the FC geometry, E k the energies of the cation that can be populated by the photon at the same geometry, and IP ik are the relative ionization potentials. Our calculations show that 19 cation states can be populated (up to about 20 eV above the X state) 11 . For the calculation of the peak intensities, I ik , we adopt an approach based on Dyson orbitals 10 . The latter are defined as where Φ i el are the electronic functions of the neutral molecule as defined above and Φ k cat the electronic functions of the cation. We calculated Dyson norms at the FC point (see ref. 9) at the CASSCF(17,12)/aug-cc-pVQZ (no state average) level of theory for the cation wave functions and CASSCF(18, 12)/aug-cc-pVQZ (no state average) for the neutral wave functions with the MOLPRO quantum chemistry package 36 . The energies of the neutral and the cation were further refined with MRCI-SD(Q) calculations, including Davidson corrections, and based on the previous CASSCF references.
If a sudden approximation is assumed, the squares of the Dyson norms, φ φ i k , are proportional to the relative ionization probabilities I ik . Ionization potentials and are reported in 1. The corresponding stick spectrum is displayed in Fig. 2. To obtain the energy resolved spectra we convoluted the stick spectra with a Gaussian envelope function G(ε) to mimic the bandwidth of the XUV probe pulse, Here σ is the standard deviation of the intensity: σ = 1.5 eV for a probe pulse of FDHM equal to 500 as. Electron kinetic energy / eV  Table 1) are labeled according to the order given in ref. 35 Let us now consider the full photoionization dynamics. Assuming a randomly oriented molecular sample, the differential cross section in the laboratory frame (LF) coordinate system is given by the following expression: jk jk jk jk jk jk 2 is the second order Legendre polynomials and θ is the angle between the direction of the electron momentum and the polarization of the electric field. Ω is the angle relative to electron emission momentum in the LF system and the two energy dependent parameters are σ jk (partial cross section) and β jk (asymmetry parameter). (The LF system defines the experiment i.e. the direction of the polarization and propagation of light as well as the direction of electron detection. The reference system is the molecular frame (MF) system in which the molecule is fixed and the electronic structure, transition dipole moment etc. calculations are performed.) Calculation of σ and β parameters require an explicit description of the continuum wave function for the final state. Neglecting interchannel coupling effects, generally very small far from thresholds, a single channel approximation of the form is generally quite accurate. Here ϕ κ −  ( ) describes an electron with asymptotic momentum κ  (and incoming wave boundary conditions, appropriate for photoionization), and A describes antisymmetrization and proper symmetry couplings. Actually it is computationally easier to work in an angular momentum basis, employing eigenstates 1 1 where the continuum wavefunctions ϕ εlm are characterized by suitable asymptotic conditions, in our case K-matrix boundary conditions, defined as lm l m l l l m m l l m lm l m , which has the advantage of working with real wave functions. Here f l and g l are regular and irregular coulomb functions. The ϕ εlm so obtained can be transformed to incoming wave boundary conditions and then to linear asymptotic momentum by standard transformation 37 The same transformation can be directly applied to the transition dipole moments. The many-particle transition dipole moment ik lm k cat lm i el ; reduces to the single particle moment involving the Dyson orbital (9) plus an additional term (conjugate term) which is generally small and is usually neglected 38 . Here γ is the Cartesian component of the dipole, D and d are the many-particle and the single particle dipole operators. From dipole moments (and the K-matrix) σ jk (ε) and β jk (ε), as well as any angular distribution from oriented molecules, can be computed according to well known formulas 37 .
In our formulation, the continuum wave function (13) is computed as an eigenfunction of the Kohn-Sham Hamiltonian defined by the initial state electron density ρ

KS lm lm
where V eN is the nuclear attraction potential, V C the coulomb potential and V XC the exchange-correlation potential defined in terms of the ground state density ρ. The latter is obtained from a conventional LCAO SCF calculation, employing the ADF program with a DZP basis 39,40 . A special basis is employed for the continuum solutions of (19). Primitive basis functions are products of a B-spline radial function 41,42 times a real spherical harmonic ilm i l m The full basis comprises a large one-center expansion on a common origin, with long range R max0 , and large maximum angular momentum, L max0 . This is supplemented by additional functions centered on the nuclei, of very short range, R maxp , and small angular momenta L maxp . A short range is necessary to avoid almost linear dependence of the basis, which spoils the numerical stability of the approach. Despite the very limited number of LCAO functions these choices ensure a very fast convergence of the calculated quantities. The basis is then fully symmetry adapted. The calculation of continuum eigenvectors is performed at any selected electron kinetic energy by the Galerkin approach originally proposed in ref. 43 and the generalized to the multichannel case 44,45 . From the energy independent Hamiltonian H and overlap S matrices continuum vectors are obtained as eigenvectors of the energy dependent matrix A(E) = H − ES with eigenvalues very close to zero. These give the correct number of independent open channel solutions, and are efficiently obtained by block inverse iteration, since they are separated by large gaps from the rest of the spectrum. Actually the more stable form A + A is currently employed 46 . Final normalization to K-matrix boundary conditions is obtained by fitting the solutions to the analytical asymptotic form at the outer boundary R max0 .
In the present calculation the LB94 V XC potential 47 was employed, due to the correct asymptotic behavior, important in photoionization. Parameters were L max0 = 12, R max0 = 25.0 a.u., with 135 B-splines of order 10, L maxp = 2, R maxp = 1.50 a.u. for the O atoms, for a total of 23013 basis functions. Such an approach, called static-DFT proves in general remarkably accurate for the description of cross sections and angular distributions 41,48,49 . In conjunction with the Dyson orbital formulation it is able to describe ionization involving multiconfigurational initial and final cationic states 38,50 . We refer to previous work for details of the implementation 41,51 . σ jk and β jk are obtained on a dense electron KE ε jk grid, so that the value at any KE dictated by the given photon energy can be accurately obtained by interpolation. With these the angularly resolved photoelectron intensity becomes: Applying the same convolution procedure as in Eq. 9 of ref. 11 we arrive to the appropriate formula of the angle resolved photoelectron spectrum: k kk k Here the ρ kk (τ) comes from eq. 6 and from now on the above expression (eq. 23) will serve as our working formula in the forthcoming part of the paper. Figure 3 displays the intensity of the ejected electrons as a function of energy and time delay between the pump and probe pulses for three different fixed values of the orientation angle, θ. It can be seen that the ionization probability is larger for smaller angles. For θ > 45° it is drastically reduced. At early times, when t delay < − 2 fs, ionization can only take place from the ground state, X. Here, two clearly distinct high intensity bands are observed within the 75-78 eV and the 80-85 eV energy intervals. These are consistent with the large Dyson norms calculated between the X state of the neutral and some of the states of the cation (see Table 1). In particular, large Dyson norms are found between X and the 1st (0.  Fig. 4, the electron emission orientation is given against the energy of the ejected electrons at several consecutive times. We observe that, up to t delay = − 1 fs, only two energy regions, (75-78) eV and (81-84) eV, exhibit significant intensity. They correspond to ionization taking place from X only. Ionization occurring from B, once t delay > − 2 fs, is characterized by the third band that appears around 88 eV and disappears slowly beyond t delay > 4 fs. Within the t delay = 1-2 fs time interval, the strengthening of the middle band reflects the combined impact of ionization occurring from both states together. Again, one clearly sees that, as a general trend, the intensity decreases monotonically as the angle between the ejected electrons and the direction of the polarization increases.

Results and Discussion
In Fig. 5, the electron emission orientation is plotted as a function of the time delay for several fixed electron energy values. Again, one observes large intensities in the (75-77) eV energy region and t delay < 0 fs time interval for small orientation angles. The latter corresponds to the lack of population of the B state resulting in ionization taking place only from the X state. For t delay > 0 fs, the decrease of the intensity indicates depletion of the X state. For > 80  eV, a joint effect of ionizations from X and B is observed, more substantially from X. Again, the shape and the structure of the band for  > 85 eV and t delay = (− 2)-6 fs is typical of ionization occurring from B.
From Fig. 5 it also appears that the angular distribution is strongly peaked along the probe field polarization, which is consistent with a high β value, close to two, for all ionizations. This is not surprising because of the high photon energy of the probe, 95 eV, which implies high kinetic energy of the outer valence ionized electrons, typically characterized by high β values, similar for all ionizations.
Finally the oscillatory patterns appearing in Figs 3 and 5 are clear fingerprints of the time dependence of the external electric field. Specifically, the pump pulse is a few-cycle pulse of width 3 fs and period 0.87 fs, centered around 260 nm (4.8 eV) in the deep UV (UV-C) domain and therefore its oscillation is faster than the nuclear motion.
In summary, the most representative signal is perhaps the upper-right panel in Fig. 3 (intensity against electron kinetic energy at different time delays for θ = 0°). It is clear that the largest temporal change in the spectrum is associated with the highest kinetic energies, from 86 to 89 eV, which are exclusively emitted from the B state, where the intensity increases significantly just after the pump pulse. Correspondingly, the decrease of the intensity after the pump is most evident in the low kinetic energy region, from 75 to 78 eV, due to the depleting of the X state, which is the dominant contribution in this energy window. Conclusions A numerical simulation protocol has been developed for describing the electron dynamics of the ozone molecule in the Franck-Condon region involving only the ground (X) and Hartley (B) electronic states in the dynamics. Assuming isotropic initial distribution for the molecular ensemble, angle resolved photoelectron spectra have been calculated for various time delays between the pump that creates the wave packet (coherent superposition of X and B) and the probe that ionizes from either X or B. This physical quantity can be measured in actual experiments and compared to our calculations.
The present results are very encouraging and call for further improvements concerning the accuracy of the dynamics simulations. Therefore, our future aim is to perform more realistic simulations upon going beyond the presently assumed limiting hypotheses: isotropic initial distribution and populations extracted at the FC geometry only. This will be manifested by two significant changes in the numerical protocol: i) after the pump pulse is off alignment of the molecular ensemble will be assumed; ii) instead of performing calculations at a single FC geometry, several other nuclear geometries will be involved in the FC region where the nuclear density has significant value too.
We stress again that given the dipole matrix elements and K-matrix, all photoionization observables can be computed, like photoionization from fixed-in-space molecules (MFPADS) or partially oriented molecules, as well as suitable averages over final detector energy and angle resolution 49 , to accurately describe any specific experimental setup. Actually the 95 eV pulse employed in the present study was suggested by an experimental colleague. With hindsight angular distribution from unoriented molecules turn out not to be very informative, given the β values close to 2 for all final states at this relatively large photon energy. Working at lower energies would produce larger anisotropies. Moreover working with oriented molecules, which is a goal actively pursued in such studies, would further much enhance anisotropies, different for each initial and final state.
The present numerical simulations clearly indicate that angle and time resolved photoelectron spectra can be used in molecular attophysics to characterize the creation of an electronic wave packet in a neutral molecule on the subfemtosecond time scale. We expect our computational study to be followed by experiments showing similar results.
As the number of experimental choices is quite large, we found it important to set up a fully ab-initio general formulation that will accommodate any specific experimental setup. We look forward to upcoming experiments to validate the theoretical framework provided here.