Role of electronic correlations in photoionization of NO 2 in the vicinity of the 2 A 1 / 2 B 2 conical intersection

We present the first ab initio multi-channel photoionization calculations for NO2 in the vicinity of the A1/ B2 conical intersection, for a range of nuclear geometries, using our newly developed set of tools based on the ab initio multichannel R-matrix method. Electronic correlation is included in both the neutral and the scattering states of the molecule via configuration interaction. Configuration mixing is especially important around conical intersections and avoided crossings, both pertinent for NO2, and manifests itself via significant variations in photoelectron angular distributions. The method allows for a balanced and accurate description of the photoionization/photorecombination for a number of different ionic channels in a wide range of photoelectron energies up to 100 eV. Proper account of electron correlations is crucial for interpreting time-resolved signals in photoelectron spectroscopy and high harmonic generation (HHG) from polyatomic molecules.


Introduction
Understanding molecular photochemical reactions is a challenging task due to both the large number of excited states that usually participate in the reaction and the various intra-molecular radiationless processes, which redistribute the charge and vibrational energy of the molecule. These include, e.g., internal conversion, isomerization, and proton or electron transfer. [1][2][3] Greater insight into these photochemical reactions, and the corresponding underlying dynamical processes, can be obtained by performing experiments in the time domain. Typically, this involves pumpprobe experiments where: (i) a pump laser pulse excites the molecule, creating a wavepacket consisting of a coherent superposition of the molecular eigenstates; (ii) the wavepacket propagates in the energetically accessible regions of the reaction coordinates; (iii) and finally, the wavepacket is probed by the most suitable technique (e.g. laser induced fluorescence, 4 resonant multiphoton ionization, 5,6 photoionization, 7 high harmonic generation [8][9][10] and transient absorption 11,12 ). Timeresolved information on the photochemical reaction is obtained by using a variable delay between the pump and the probe.
Interpreting the information recorded in the pump-probe signal is, however, challenging. In particular, the dynamics induced by the pump step is often entangled with the complexity of the final state, or a superposition of such states, accessed by the probe. Here, we focus on how this challenge manifests itself when photoelectron spectroscopy or high harmonic generation is used as a probe step.
Time-resolved photoelectron spectroscopy (TRPES) is a powerful and wide-spread technique for investigating photochemical reactions. 13 There are several reasons for using photoelectron spectroscopy as a probe in time-resolved experiments. First, XUV light can be used to photoionize valence electrons, thus providing valuable insight into the chemically active molecular orbitals of the target. Second, photoionization is always an allowed process, so it can be used to probe, in principle, all the reaction coordinates. Third, low intensity XUV radiation can be treated using first-order perturbation theory, thus greatly simplifying the theoretical treatment of the ionization process.
Moreover, new FEL facilities and HHG sources are now able to deliver broadband XUV radiation with B10 fs duration, giving access to the fastest molecular processes with unprecedented temporal resolution. 7,14 These advances have stimulated a growing number of experimental and theoretical efforts for investigating several aspects of molecular dynamics, e.g., to map the electronic character of molecules along the reaction path, [15][16][17] to shed light on the deactivation pathways of complex molecules, [18][19][20] and most relevant for this work, to elucidate the non-adiabatic dynamics in polyatomic molecules. [21][22][23][24][25][26] High harmonic generation (HHG) provides complementary information to photoelectron spectroscopy. Indeed, the radiative photorecombination step in HHG, which is responsible for the emission of a high-frequency photon, is an inverse of the photoionization step in photo-electron spectroscopy. Notably, when several final states of the ion participate in photoionization, the photoelectron spectra contain incoherent superposition of the contributions correlated to these final ionic states. In contrast, high harmonic generation encodes the interference of these ionic channels and hence the relative phases between them, which determine the dynamics of the hole created by ionization and can be reconstructed from the measurements. [27][28][29][30][31][32][33] The interferometric property of HHG places particular premium on a balanced and accurate description of the photorecombination step from the number of different participating ionic channels, therefore demanding a high-level modelling of correlation and polarization. Additionally, for complex molecules the use of longer-wavelength (mid-infrared) pulses of the driving field is necessary. 34,35 This in turn leads to larger energies of the recombining electron. Therefore, accurate theoretical modelling of HHG requires the extension of the energy region available in photoionization calculations, which is a challenging task for current state-of-art high-quality multichannel calculations. The data presented in this work satisfy these demanding requirements and open the way towards the calculations of HHG from NO 2 including the effects of nuclear motion.
To describe the photoionization step, the existing state-ofthe-art wavefunction-based theoretical studies of TRPES 21,36 use a high-quality configuration interaction (CI) description of the bound molecular states, but a single Slater determinant description of the residual ionic states. This approach neglects the electronic correlation in the ions, an assumption that is particularly crude for treating highly excited ionic states, or the region of the reaction coordinates where non-adiabatic dynamics may affect the ion. We address this limitation here using the multichannel R-matrix method. The calculations are performed using the UKRmol + package, [37][38][39] which allows us to describe the electronic correlation in several neutral and ionic molecular states via configuration interaction. In this paper, we concentrate exclusively on the photoionization process in TRPES. The coupling of the photoionization calculations performed here with nuclear dynamics for TRPES or time-resolved high harmonic spectroscopy is the subject of a future work.
In this work, we extend our previous results for the photoionization of NO 2 in the equilibrium geometry, 40 and demonstrate the first application of the R-matrix method for the photoionization of polyatomic molecules for several nuclear geometries. Here, we restrict our analysis to photoionization from the ground and the first excited neutral states of NO 2 , leaving the residual ion in the lowest energy singlet and triplet ionic states for an asymmetric stretching varying from 0 to 0.2 Å and bending angles from 801 to 1601. We are interested in studying photoionization transitions close to the conical intersection between the ground and the first excited neutral states of NO 2 , 21 which occurs at a bending angle of E1071. Previous quantum chemical calculations 41 have shown the existence of an avoided crossing between the two lowest energy singlet ionic states of NO 2 also at E1071. Consequently, a precise evaluation of the photoionization matrix elements for such a transition requires a careful description of electronic correlation for both neutral and the singlet ionic states. This paper is organized as follows: in Section 2, we briefly describe the principles of the R-matrix method. In Section 3, we present the details of the calculations and in Section 4, we report the orientationally averaged results and the photoelectron angular distributions (PADs) for the photoionization of the NO 2 molecule in a range of nuclear geometries. Finally, in Section 5, we present the conclusions and future directions of this work.

Photoionization calculations using the R-matrix method
The molecular-frame photoelectron angular distribution for photoionization from the F N i bound state and leaving the ion in the f NÀ1 f state is given by where a is the fine structure constant, N is the number of electrons in the molecule, a 0 is the Bohr radius, o is the photon energy in atomic units, d is the dipole operator in the length gauge, and k f is the momentum of the ejected electron. C (À) N f are continuum wavefunctions which satisfy incoming wave boundary conditions. In the case of a randomly oriented ensemble of molecules and linearly polarized laser fields, eqn (1) becomes where b is the asymmetry parameter, s fi is the partial photoionization cross-section, P 2 is the second-order Legendre polynomial and y is the electron ejection angle measured wrt the polarization of the incident photon in the laboratory frame. The R-matrix method divides the configurational space into an inner (r r a) and an outer (r 4 a) region. The R-matrix radius, a, is chosen in such a way that only the continuum electron reaches the R-matrix boundary. In the inner region, both the continuum and the bound-state wavefunctions are given in terms of the basis functions, c N k , as where A (À) kf (k f ) and B ik are the coefficients of the expansion. The basis functions c N k are given in terms of the close coupling expansion c N k x 1 ; :: where Z ij are the continuum orbitals orthogonalized with respect to the target orbitals, A is the antisymmetrisation operator, b kp and a kij are coefficients determined by diagonalization of the Hamiltonian. In the second summation, w (N) p are the so-called L 2 configurations, which improve the description of electron correlation, polarization and help in converging the partial wave expansion by allowing the continuum electron to enter high angular momentum states close to the nuclei. They are formed by allowing the continuum electron to occupy target molecular orbitals. The reason they need to be explicitly added is that the one electron continuum functions are orthogonalized to the target molecular orbitals. Restricting the continuum electron to the continuum orbitals would neglect an important part of the Hilbert space.
Following the assumption that only the continuum electron reaches the R-matrix boundary, the outer region solutions are found by matching with known asymptotic solutions. The coefficients of the continuum and the bound state wavefunctions (A (À) kf (k f ) and B ik , respectively) are then found after matching the solutions of the inner and the outer region at the R-matrix boundary (see ref. 39 for further details).

Calculations
As illustrated in Fig. 1, NO 2 has three vibrational modes: symmetric stretch, scissors and asymmetric stretch. The latter is defined as the simultaneous shortening of one of the N-O bond lengths by r a /2, and the stretching of the other N-O bond length by r a /2. The ground state equilibrium geometry has C 2v symmetry with a bond distance of r s E 1.25 Å and an O-N-O angle of g E 1331. In this paper, the molecule is fixed in the zx plane, with the z-axis coinciding with the molecular axis (see Fig. 1). We consider here variations of the asymmetric stretch (r a ) and bending angle (g). For this reason, the calculations are performed in the C s point symmetry group, since a non-zero asymmetric stretch breaks the symmetry of the molecule with respect to reflection in the zy plane. Consequently, the ground and first excited neutral states of NO 2 , which belong to the A 1 and B 1 irreducible representations in C 2v (only for r a = 0), respectively, are both mapped into A 0 symmetry in C s .

Molecular orbitals and target models
In C s symmetry, the ground electronic configuration of NO 2 at the equilibrium geometry is given by with 1-3a 0 orbitals being the core orbitals and the valence space extending up to the 12a 0 , 3a 00 molecular orbitals.
In our calculations, we generate using MOLPRO 43 the molecular orbitals on a grid of geometries using a state averaged MCSCF procedure using the 6-311G** atomic basis. 44 Hartree-Fock orbitals at the equilibrium geometry were used as an initial guess for the state averaged MCSCF procedure. We included in the state averaging the (1-2) 1 3 A 00 ionic and the (1) 2 A 0 and (2) 2 A 0 neutral states. Both neutral states had a weight of 30% in the averaging. Following our previous work on the photoionization of NO 2 at the equilibrium geometry, 40 we defined an active space consisting of the 6-10a 0 and 1-3a 00 orbitals, and tested its two different variants: model-1, which includes single excitations to the 11-14a 0 orbitals, and model-2 that extends the single excitations up to the 11-15a 0 , 4a 00 orbitals (see Table 1). Note that in our approach, the target configuration state functions (CSFs) are coupled to the continuum electron via the close-coupling equation, see eqn (5). Consequently, it is important to use an optimal number of target CSFs, in order to keep the size of the N-electron Hamiltonian manageable and also to obtain an accurate description of the continuum and bound wavefunctions.

Inner and outer regions
The continuum orbitals were generated using Gaussian-type orbitals (GTOs), with angular momentum up to l = 6 and optimized to represent the continuum up to 90 eV. The exponents of the continuum GTOs centered at the center of mass of the molecule are listed in Table 2. The molecular integrals in   the inner region were computed in quadruple precision. This allowed us to avoid numerical linear dependency problems and retain all the continuum functions in the basis. For each geometry, the 16 lowest energy virtual orbitals were included in the calculations for improving the description of polarization in the inner region. Convergence for all models, with respect to the number of states included in the close coupling expansion, was achieved for calculations with at least 80 target states. Due to computational limitations, we only performed convergence tests at the photoionization level for a few critical molecular geometries at the extremes of the angular grid. As indicated in Fig. 3 and 4, there is a good agreement on the asymmetry parameter and partial cross-sections between calculations using model-1 and model-2. Since model-1 is computationally cheaper than model-2, the former was chosen for the photoionization calculations over the whole grid of nuclear geometries.
We used the Perl scripts 45 to prepare the input files and to drive and monitor the complex sequence of R-matrix programs which have to be executed to complete the calculation for multiple geometries. The scripts were originally developed for the calculations of electron-molecule collisions and subsequently adapted to enable photoionization and execution on large computer clusters. The structure of the scripts comprises definitions, for each program, of the input generation, program execution and output analysis, followed by a specification of the calling sequence of all programs. This system is therefore very flexible with respect to the inclusion of new programs in the suite and is very helpful in organizing large amounts of output data.

Orientationally averaged results
4.1.1 (1) 1 A 0 residual ionic state. In Fig. 5(a) we show the crosssections for photoionization from both neutral states leaving the ion in the (1) 1 A 0 state. The most interesting feature in Fig. 5(a) is the sudden variation of the cross-section at g E 1071 for the (1) There is a strong correlation between the partial photoionization cross-sections and the dominant electronic configurations Fig. 6 Partial photoionization cross-sections (in megabarns) in logscale (a) and asymmetry parameters (b) for 801 r g r 1601 for the triplet final ionic state. In each subplot, the left column (right column) corresponds to the photoionization from the (1) 2 A 0 ((2) 2 A 0 ) neutral state. In each subplot, the first row corresponds to r a = 0 Å, the second to r a = 0.1 Å and the third to r a = 0.2 Å. of the neutral (initial) and cationic (final) states. Changes in molecular geometry can lead to changes in the dominant electronic configurations, thus leading to variations in the partial photoionization cross-section. This is observed in Fig. 5(a), where far from the conical intersection, the (1) 1 A 0 ionic and (1)  Note that far from the conical intersection, the leading electronic configuration of the (1) 1 A 0 ionic state is always a 1h state with respect to the (1) 2 A 0 neutral state. Consequently, we observe a large magnitude of the partial photoionization crosssection for the (1) 2 A 0 -(1) 1 A 0 photoionization transition for g o 1001 and g 4 1151 (see the left column of Fig. 5(a)). On the other hand, far from the conical intersection, the leading electronic configuration of the (1) 1 A 0 ion is always a satellite state (2h-1p) with respect to the (2) 2 A 0 neutral. The photoionization transitions are then mediated through at least one nondominant electronic configuration, thus leading to a low magnitude of the partial photoionization cross-section for the (2) 2 A 0 -(1) 1 A 0 transition (see the right column of Fig. 5(a)).
Close to the conical intersection, the neutral states are strongly coupled, thus, we observe that the partial photoionization cross-section decreases for the (1) 2 A 0 -(1) 1 A 0 photoionization transition, and increases for the (2) 2 A 0 -(1) 1 A 0 photoionization transition. We note that the avoided crossing between the (1) 2 A 0 and the (2) 1 A 0 ionic states at g E 1071 also leads to a strong coupling between these states. Consequently, photoionization calculations that do not describe correctly the effect of electronic correlation in the (1) 1 A 0 state fail to accurately evaluate the (1) 2 A 0 -(1) 1 A 0 and (2) 2 A 0 -(1) 1 A 0 photoionization transitions at g E 1081.
In Fig. 5(b) we plot the asymmetry parameter for photoionization from both neutral states leaving the ion in the (1) 1 A 0 state. At photon energies below 50 eV, we notice a pronounced change of the asymmetry parameter as g passes by the conical intersection. The increased sensitivity at low energies goes hand in hand with the variation of the photoionization cross-sections at low energies. On the other hand, at photon energies above 50 eV there is almost no signature of the change in electronic configurations of the neutral or the ionic states for 801 o g o 1601.
4.1.2 (1) 3 A 0 residual ionic state. As shown in Fig. 6(a), the partial photoionization cross-sections are approximately homogeneous for the (1) 2 A 0 -(1) 3 A 0 photoionization transition. This can be understood by the same analysis reported in the last section: the (1) 3 A 0 state has a . . .(9a 0 ) 1 (10a 0 ) 1 leading electronic configuration in the whole geometric grid considered here. Thus, it is a 1h state with respect to the leading electronic configuration of both neutral states for all geometries considered here.
On the other hand, the asymmetry parameter for the (1) 2 A 0 -(1) 3 A 0 and (2) 2 A 0 -(1) 3 A 0 transitions shown in Fig. 6(b) exhibits a sudden change in amplitude when crossing the conical intersection at low photon energies. This indicates that even though the dipole transition strength is approximately constant for all geometries considered, the angular distribution of the ejected electrons varies considerably as the molecule crosses the conical intersection. This is confirmed in Fig. 7, where we plot the Dyson orbitals in the zx plane for the (1) 2 A 0 -(1) 1 A 0 and (2) 2 A 0 -(1) 1 A 0 photoionization transitions in the vicinity of the conical intersection. We can see that as the system goes through the conical intersection, the Dyson orbitals for the (1) 2 A 0 -(1) 1 A 0 and (2) 2 A 0 -(1) 1 A 0 transitions are swapped. The result is the appearance and disappearance of the nodal plane, respectively, along the molecular axis in the two Dyson orbitals as g varies from 1081 to 1031. This shows that even for a non-oriented ensemble of molecules, one could use the asymmetry parameter to identify changes in the electronic configuration of the neutral states.

Photoelectron angular distributions
In this work, we calculated the PADs for parallel and perpendicular alignments for transitions from the (1) 2 A 0 and (2) 2 A 0 neutral states to the (1) 1 A 0 and (1) 3 A 0 ionic states in the range of molecular geometries of 801 r g r 1601, 0 Å r r a r 0.4 Å. However, we will restrict ourselves to the discussion of the PADs for the (1) 2 A 0 -(1) 1 A 0 photoionization transition near the conical intersection and only examine emission in the xz plane, which is the preferred plane for photoemission 8 in NO 2 .
In Fig. 8, we show the PADs for the (1) 2 A 0 -(1) 1 A 0 transition for parallel and perpendicular alignments at bending angles of 1151, 1081, and 1031 and asymmetric stretchings of 0 Å, 0.1 Å and 0.2 Å. The emission angle is measured from the molecular axis in the zx plane. At r a = 0 Å, we observe a node along the molecular axis in the PADs for g o 1071 in the case of parallel alignment ( Fig. 8(a)) and for g 4 1071 in the case of perpendicular alignment (Fig. 8(b)). The origin of this node can be explained by analyzing the photoionization transition in the C 2v point group which, at r a = 0, includes the invariance of the molecule upon reflection in the zy plane. In C 2v , the (1) 2 A 0 neutral state maps into the (1) 2 A 1 state for g 4 1071, and into the (1) 2 B 1 state for g o 1071. Since, in the case of parallel alignment, the dipole operator has A 1 symmetry, the scattering states must necessarily have A 1 symmetry for g 4 1071 and B 1 symmetry for g o 1071. Consequently, since the singlet state (1) 1 A 0 maps into A 1 in C 2v , the continuum electron has to belong to A 1 symmetry for g 4 1071 and to B 1 symmetry for g o 1071. As the states with B 1 symmetry are antisymmetric under reflection in the zy plane, the continuum electron exhibits a node at the molecular axis for g o 1071 (see Fig. 7  and 8(a)). A similar analysis shows that the reverse occurs for perpendicular alignment, i.e., the node at the molecular axis occurs for g 4 1071 (see Fig. 8(b)).
The PADs also can be used to identify the main partial wave channels that contribute to the photoelectron signal. Taking again as an example the PAD for parallel alignment and r a = 0, the continuum electron necessarily belongs to B 1 symmetry for g o 1071 (see the PADs corresponding to g o 1071 and r a = 0 in Fig. 8(a)). In B 1 symmetry, the partial waves can have any value of l and m must be odd. However, an examination of the Legendre polynomials for this set of quantum numbers demonstrates that only partial waves with an odd l value have a maximum at an angle of 901 with respect to the molecular axis. Consequently, we conclude that the partial waves with l = 3 and l = 5 dominate the PADs for photon energies from 30 to 50 eV.

Conclusions
We performed NO 2 photoionization calculations for a range of geometries close to the conical intersection between the (1) 2 A 0 and (2) 2 A 0 neutral states. We analyzed partial cross-sections, asymmetry parameters, Dyson orbitals and PADs for photoionization leaving the ion in the (1) 1 A 0 and (1) 3 A 0 states. At r a = 0 Å, we identify a peak in the partial cross-section for photoionization from the (2) 2 A 0 state leaving the ion in the (1) 1 A 0 state, which appears due to the conical intersection in the neutral and the avoided crossing between the (1) 1 A 0 and (2) 1 A 0 ionic states at 1021 o g o 1071. Away from this region the signal is low. This is due to the character of the ionic state changing from a satellite to a main line state when passing through the conical intersection (g E 1071), and then back to a satellite state when passing through the avoided crossing (g E 1021). The opposite happens with ionization from the (1) 2 A 0 state with a dip showing up at 1021 o g o 1071. The cross-section for photoionization leaving the ion in the (1) 3 A 0 state is relatively homogeneous throughout the range of geometries investigated since this ionic state is always a 1h state with respect to the initial neutral state. In both cases, the asymmetry parameter exhibited a clear signature of the molecular orbitals from which the photoelectron was ejected, i.e., a significant change in magnitude as the molecule goes through the conical intersection, especially at energies below 40 eV.
To the best of our knowledge, we presented in this paper the first photoionization calculations over a large grid of nuclear geometries using robust CI methods for describing both the neutral and the ionic states. We find that the effects of configuration-interaction are especially important to describe accurately photoionization leaving the ion in the (1) 1 A 0 state (see Fig. 5(a)) and should be taken into account when interpreting the time-resolved photoelectron spectra, which map the dynamics taking place on multiple potential energy surfaces of the molecule and their interplay via the conical intersection, 21 since the range of geometries for which the configuration mixing is important is relatively large.
Our approach allows us to perform calculations of TRPES using broad-band XUV probe laser pulses, which leave the ion in highly excited states after ionization, thus providing rich information regarding the electron dynamics of polyatomic molecules in non-adiabatic regions of the reaction coordinates. In addition, we highlight that our results are also relevant for high harmonic spectroscopy, 8 since the recombination step in HHG is strongly linked to the one-photon ionization described here (see ref. 33 and 46-48).