Natural ionization orbitals for interpreting electron detachment processes

Theoretical simulation of electron detachment processes is widely used to interpret and understand experiment in a diverse range of chemistry and physics applications.^1–13 In order to understand how the wavefunction changes upon ionization, two families of methods exist—those based on Koopmans’ approach¹⁴ in which the hole is described by an orbital in the N electron calculation and those based on a delta self-consistent field (ΔSCF) between an N and N − 1 electron calculation. Formally the ionization energy is equal to the negative of the orbital energy of the ionized electron in both a restricted (open-shell) Hartree-Fock (HF) and Density Functional Theory (DFT) calculations,¹⁵ although in DFT calculations this value is highly dependent on the functional used and often significantly incorrect.¹⁶ ΔSCF methods typically calculate the ionization potential by removing the ionized electron excited into a continuum state and explicitly calculating the N − 1 electron system. Such simulations usually employ DFT, owing to the excellent balance of simplicity, cost, and accuracy. However, the generality of the approach enables any ground state (GS) method to be employed including post-self-consistent field (SCF) methods. While analysis of the wavefunction response to ionization in the Koopmans method and models derived from it is simple, in studies that employ ΔSCF inclusion of relaxation effects accompanying ionization mean that the resulting canonical orbitals are usually difficult to interpret as they are often spatially diffuse and may not easily map directly to the ionization event.

In order to better interpret and understand the nature of electronic transitions, unitary transformations of the orbital space to which the energy is invariant can be applied. Prior approaches based on such transformations have been reported for analyzing simulations of electron transitions. In the limiting case of electronic excitation to a non-continuum state, Martin’s Natural Transition Orbitals (NTO)¹⁷ and the attachment-detachment density analysis of Head-Gordon et al.¹⁸ have been proposed. Ziegler et al. have developed the natural orbitals for chemical valence (NOCV) scheme for analysis of chemical bonding which, in a similar approach to this work, utilizes a natural orbital (NO) orbital transformation of a difference density, Δ_P.¹⁹ The NOCV approach constructs Δ_P from the difference of total and fragment (promolecular) densities while this study uses a Δ_P constructed from the density of initial and final states of an electronic transition.

Methodologies with specific application to ionization calculations in which the size of the occupied-virtual (ov) orbital space changes include the ΔSCF method of Martin and Davidson.²⁰ This approach, which is directly analogous to the NTO scheme, constructs an N_elec × (N_elect − 1) overlap matrix between molecular orbitals (MOs) in the ground and ionized states, where N_elec is the number of electrons in the ground state. Singular value decomposition of this matrix accomplishes a corresponding orbital transformation²¹ with eigenvalues identifying the overlap between MOs in the ground and excited states. The ground state orbital without a corresponding excited state occupied orbital can be labeled the orbital from which the ionized electron originated.

The study of excitations beyond single reference wavefunctions can be achieved through calculation of Dyson orbitals.^22–25 Dyson orbitals are calculated from the overlap between the N electron wavefunction and the N − 1 electron wavefunction of the ionized state and generalize Koopmans’ theorem to correlated wavefunctions obtained through post-SCF methods. While the Dyson orbital approaches are incredibly powerful, they exist within a wavefunction framework and cannot be easily extended to the widely used DFT family of methods. Furthermore, the requirement for a correlated wavefunction limits the size of system that can be tackled with the Dyson orbital approach. Approaches for approximating Dyson orbitals from DFT calculations to reduce computational scaling have been proposed by Dauth et al.²³ However, beyond wavefunction analysis, Dyson orbitals have had particular application in electron propagator theories (EPTs), the formally equivalent one-electron Green’s function, and equation-of-motion (EOM) approaches to obtain accurate ionization energies. Most closely related in scope to this work are the developments of EOM-coupled cluster (CC) by Krylov et al. to the study of ionization processes^26–31 and the electron propagator theory (EPT) work of Ortiz et al.^32–38 

Experimentally, (e,2e) or electron momentum spectroscopy (EMS) provides the ability to probe the electron density of individual orbitals during ionization events.^39–47 Although EMS can access the complete valence shell, it does this with lower resolution than more widely used photoelectron spectroscopy (PES). However, the EMS cross section is proportional to the square of the Dyson orbital associated with the transition and in this regard has the potential to image the electron/hole density of individual orbitals.

Here we describe an alternative approach that is entirely general to the methodology used for computing ground and ionized states in electron detachment calculations. Using an NO based approach for ionization events, we propose a scheme we refer to as Natural Ionization Orbitals (NIO)s. More specifically the approach uses a NO analysis scheme to find a compact representation of electron detachment and indicates the key orbitals involved in the process through calculated occupation change eigenvalues. Pulay and Hamilton were likely the first to suggest using a single set of NOs obtained from an unrestricted Hartree-Fock (UHF) calculation,⁴⁸ although the one-electron density matrix was used to yield the correlation space for Complete Active Space (CAS) calculations. We emphasize that this differs from the approach detailed in this article, in which a density difference approach is employed for analysis of simulated electron excitation processes, although we do posit that NIOs may have application to multireference approaches based on density difference formalisms such as difference dedicated configuration interaction (DDCI).^49,50 We show that the method is able to differentiate between different types of ionization processes and provide key mechanistic information that would be difficult to determine by visualization of canonical orbitals or employing population analysis schemes. The method is then compared to existing approaches and the strengths and weaknesses contrasted.

In the case of a single electron detachment process the density difference matrix, Δ_P, can be computed for the vertical detachment processes by taking the difference of one-particle density matrices (1PDMs) computed in the ground, P_g, and ionized, P_i, states,

The change in the number of electrons Δ_elec can then be computed from

where S is the atomic orbital (AO) overlap matrix and the angled brackets denote the matrix trace. Utilizing the fact that the trace of Δ_PS is invariant to unitary transformations we compute a set of orbitals ϕ_p which we term NIOs from the canonical orbitals, {ψ_p},

The matrix U can be determined by solving the eigenvalue equation

where δ_elec is a diagonal matrix containing the occupation change number for each NIO. The sum of these individual occupation change numbers gives the total change in the number of electrons Δ_elec. The eigenvectors U can then be back transformed to the AO basis to give the NIO coefficients V.

Extension of this methodology to unrestricted calculations in which α and β electrons are treated separately is carried out by using the α 1PDM difference matrix, $ΔPα$⁠, and β 1PDM difference matrix, $ΔPβ$⁠, to compute two sets of NIOs for each spin type, V^α and V^β.

Conceptually the NIOs can be considered analogously to the standard NO scheme of Löwdin.⁵¹ Transformation to the set of NOs has general applicability by improving the convergence behavior of full configuration interaction (FCI) calculations. The set of NOs yields an associated set of occupation numbers as eigenvalues and the natural expansion of the configuration interaction space is over NOs that differ in their occupation number from zero. As such, the NOs represent the fastest converging set of orbitals for the FCI expansion and in the limiting case of occupation numbers either zero or one collapse to the single Slater determinant. We propose that NIOs can be considered similarly as the natural basis for describing transitions from a bound state to the continuum state. The NIOs with the largest occupation change numbers give the most compact basis for describing the transition.

To ascertain the value of NIO analysis, we compare NIO orbitals to those obtained from both the Koopmans approach (i.e., the Highest Occupied Molecular Orbital (HOMO)) and calculated the Dyson orbitals. The EPT methods of Ortiz et al.^52–56 are used to generate Dyson orbitals D which can be expanded in the spin-orbital basis as D = UC. The C are the MO coefficients, while the U_nr are related to the pole strength P_n, which determines the likelihood of transition from orbital n, according to

Within the diagonal self-energy approximation,⁵⁷ the Dyson orbital n can thus can be calculated from the relevant HF orbital according to $Dn=PnCn$⁠. The Dyson orbital associated with each ionization energy is normally closely related to the canonical orbitals, particularly for ionization from the ground state. In all the systems modeled here the pole strength was greater than 0.9 indicating the validity of Koopmans approach and the diagonal self-energy approximation.

In order to illustrate the utility of the NIO procedure, we demonstrate its application in three different electron detachment analyses using an implementation within a local development version of Gaussian⁵⁸ which interfaces with a separate standalone code. NIO analysis is carried out in both unrestricted and restricted (open-shell) frameworks using HF and various DFT functionals with different basis sets as indicated throughout the text. The valence orbital ionization calculations display various degrees of orbital relaxation and correlation effects, and we compare them to Dyson and Koopmans (HOMO) orbitals to demonstrate the ways different methods present the wavefunction response to ionization.

The first system that we examine is electron detachment of the ¹A₁ ground state (GS) of formaldehyde (CH₂O). Figure 1 shows the relevant orbitals for the discussion within this section. The top row of Fig. 1 shows the HOMO of HF calculations obtained with both orbital spatial and spin symmetry constraints (left) and unrestricted symmetry (right) using the 6-311++G basis set. The HOMO orbitals in the HF calculations correspond to the electron hole within Koopmans’ theorem. Strictly, Koopmans’ theorem is only valid with the context of restricted HF calculations. However, both UHF and unrestricted DFT B3LYP calculations reveal that no symmetry breaking occurs in the ¹A₁ state of formaldehyde. In both restricted and unrestricted cases Koopmans’ theory predicts the electron to be ejected from a b₂ symmetry orbital with mixed σ(H) and π(O) character. However, while Koopmans’ theorem is valid in the limit that the orbitals of the ionized state are identical to those of the ground state, it assumes no orbital relaxation in the ionized state and has an incomplete description of electron correlation owing to the underlying HF calculation.

In order to properly account for correlation effects, we compute Dyson orbitals from the EPT approach within the diagonal self-energy approximation (Fig. 1, bottom). The EPT calculations provide the correlation correction to the HF reference orbitals. The calculated pole strengths P_n in these calculations were found to be 0.926 80 using both closed-shell and open-shell models, yielding a Dyson orbital corresponding to scaling of the HOMO of 0.962 70. The Dyson orbitals thus strongly resemble the HOMO molecular orbitals demonstrating that in this case Koopmans’ theorem provides a reasonable picture of electron detachment.

Unlike the Koopman and Dyson orbital approaches, NIO analysis requires the calculation of both the GS and the ionized wavefunction and thus it accounts for orbital relaxation effects accompanying ionization in a ΔSCF picture. The resulting NIOs are, therefore, dependent on the correct states being obtained in both calculations. The central row of Fig. 1 shows the hole in the NIO analysis computed using both restricted and unrestricted B3LYP using 6-311++G basis set. Calculation of the ionized state using restricted open-shell B3LYP was found to converge to the ²B₁ state. Attempts to converge to the ²B₂ state using restricted open-shell HF were unsuccessful, and in Fig. 1 we show the hole of the ¹A₁ → ²B₁ + e⁻ transition computed with NIO analysis. The ²B₁ state is determined to result from electronic detachment from the CO π orbital. The ²B₂ state was obtained without issue using unrestricted B3LYP and the NIO hole shown in Fig. 1 resembles both the HOMO and the Dyson orbital of the ¹A₁ → ²B₂ + e⁻ transition. The validity of the Koopman picture is supported by, as well as the high pole strength in the EPT calculations, the small δ_elec eigenvalues for electron relaxation (paired) orbitals, which are all smaller than 0.15 indicating a small level of orbital relaxation accompanying ionization.

While formaldehyde provides a simple test system in which Koopmans’ model is valid and orbital relaxation effects are small, in more complicated systems the effect of correlation and orbital relaxation would be expected to significantly degrade the qualitative picture given by the HOMO. We thus turn to the ²A″ → ¹A′ + e⁻ transition of the nitromethyl radical which exhibits more significant orbital relaxation accompanying the ionization. The HOMO orbitals of restricted and unrestricted HF/6-31G(3d,3p) calculations are shown in Fig. 2, top row. Unlike the formaldehyde example there is a significant difference between the restricted and unrestricted HOMO owing to symmetry breaking in the unrestricted calculation as demonstrated by the total spin squared expectation value $〈Sˆ2〉$ of 1.17 (for a spin pure doublet this value should be 0.75). The main difference between the restricted and unrestricted HOMO orbitals is the sign of the p atomic orbital on the backwards-facing oxygen atom, which changes the position of the node from between the CH₂ and NO₂ groups in the Restricted Open-Shell Hartree-Fock (ROHF) calculation, to between the H₂CNO and O groups in the UHF calculation. Thus in the ROHF calculation the ejected electron mainly resides in a p-type orbital on the CH₂ group, while in the unrestricted calculation the ionized electron is located within the π orbital delocalized over the H₂CNO moiety.

The Dyson orbital computed using the unrestricted HF reference is shown in Fig. 2, bottom right. The calculated pole strength P_n in this calculation is 0.911 42 and thus the resulting scale factor of the HOMO is 0.9547. The resulting orbital strongly resembles that predicted from Koopmans’ theorem and implies that Koopmans’ theorem is qualitatively valid, a conclusion that is further supported by the relatively modest change in ionization energy between EPT and UHF calculations of 0.516 eV. The Dyson orbital itself supports the idea that correlation effects result in an ionized electron that is more localized on the H₂CNO group. The EPT model is not implemented for restricted open-shell models and thus we are unable to compute Dyson orbitals using the ROHF reference.

The NIO hole calculated from restricted open-shell HF with a 6-31G(3d,3p) basis set on both the ²A″ and ¹A′ states is shown in Fig. 2, bottom left. This calculation indicates that there is very little orbital relaxation accompanying the ionization as all other NIO pairs have occupation change numbers δ_elec of less than 0.08. The NIO hole strongly resembles the HOMO orbital in this case. Calculations of NIOs using a restricted (open-shell) formalism effectively force spatially equivalent NIO pairs in the α and β sets with equal occupation change numbers δ_elec except from one pair in each spin set that corresponds to the NIO hole state and its partner in the opposite spin set. Conjectural evidence based on the few systems we have studied points to the fact that NIO analysis based on a restricted (open-shell) calculation tends to indicate much smaller amounts of electronic relaxation upon ionization than calculations based on unrestricted calculations. This presumably results from enforced spatial pairing of α and β spin orbitals limiting the amount of spin polarization in response to ionization. Although the restricted open-shell framework accounts for relaxation effects absent in the Koopmans picture, electron correlation effects are ignored in both pictures. In contrast, NIOs based on spin-unrestricted calculations do account for both orbital relaxation and some electron correlation contributions (though often at the cost of introducing spin contamination).

Fig. 3 shows the NIOs with significant occupation change numbers δ_elec computed using UHF/6-31G(3d,3p). The NIO hole (1α) is entirely localized in the non-bonding π MO over the H₂CNO moiety. This differs from the UEPT Dyson orbital and the UHF HOMO orbital as there is no contribution from the p orbital of the backwards-facing oxygen atom. The origin of the difference is that the NIO approach includes relaxation effects that are absent from the Koopmans model and the correlation corrected Dyson orbitals. Separation of the electron hole and simultaneous relaxation of the wavefunction in response to ionization obtained through NIO analysis limits the hole to the non-bonding π orbital which allows a more intuitive understanding of where the ionized electron originates. The contribution of the backwards-facing oxygen p orbital is found in the NIOs that describe relaxation of the wavefunction. Relaxation in the β spin set is shown in the 1β – 3β orbital pair (Fig. 3, right), indicating that β electron density moves from the backwards-facing oxygen p orbital to the CH₂ moiety. Relaxation in the α spin set defined by the 2α − 3α NIOs (Fig. 3, left) demonstrates that electron density moves from the π bonding orbital over the H₂CNO moiety to a NO π antibonding orbital, with a greater contribution from the nitrogen p orbital. Based on the few systems we have calculated, unrestricted NIOs tend to localize all relaxation into a single NIO pair for each electron spin, producing five relevant NIOs including the hole NIO, demonstrating the highly convergent nature of NIOs.

In the corresponding ionization orbital (CIO) model of Martin and Davidson,²⁰ the overlap between orbital pairs can be related to the probability that an ionization event to the continuum will be accompanied by additional finite excitations. Such additional finite excitations define shakeup and shakeoff states and, owing to the fact that the eigenvalues resulting from the singular value decomposition of the overlap matrix are physically related to the NIOs, the eigenvalues δ_elec similarly imply the occurrence of shakeup and shakeoff states. In the case of nitrodimethyl radical the δ_elec eigenvalues associated with relaxation are sufficiently small that we consider this transition to be of a one-electron nature; however, it is possible to find orbital pairs with very large δ_elec values approaching unity.

The final example of the use of NIOs that we describe in this report is that of the ³A′ → ²A″ + e⁻ electron detachment in $MoNbO4−$⁠. Figure 5 illustrates the significant NIOs as determined by δ_elec > 0.1. This was computed using UB3PW91 with the SDDPlusTZ basis set consisting of SDD pseudopotentials on molybdenum and vanadium, an augmented SDD valence basis set on the two transition metals, and aug-cc-pVTZ on the oxygen atoms.⁵ This case motivated the development of NIOs as it was unclear whether the ionized broken symmetry state that was found could be formed from a single electron transition or whether it was formed as a result of a “shake-up” transition in which a one electron ionization is simultaneously accompanied by a finite transition of a second electron. Interpretation of the canonical MOs is difficult due to significant electronic rearrangement; it was not possible in our initial analysis to identify a one-to-one correspondence between anion and neutral states. Thus it was not even possible to identify a particular MO from which the electron detachment originated, let alone determine if any accompanying transitions had taken place. Calculated Mulliken spin densities appeared to imply that all electrons on molybdenum were paired and two electrons on vanadium were unpaired in the anion while in the neutral there were two unpaired α electrons on molybdenum and an unpaired β electron on vanadium. The Mulliken spin densities therefore appeared to imply that some form of shake-up transition was involved, but spin densities are unable to provide information on the mechanism of any such electron transfer and indeed the reliability of such values is not certain when using diffuse basis functions owing to the limitations of Mulliken partitioning.

The HOMO orbitals computed using UHF and UB3PW91 are shown in Fig. 4 along with the Dyson orbital computed from the UHF reference. These calculations were carried out using UB3PW91 with the SDDPlusTZ basis set consisting of SDD pseudopotentials on molybdenum and vanadium, an augmented SDD valence basis set on the two transition metals, and aug-cc-pVTZ on the oxygen atoms.⁵ The UB3PW91 HOMO shows that α electron density is removed from 3d_z² orbitals on both vanadium and molybdenum (Fig. 4, left). The UHF HOMO shows a similar picture with greater density taken from molybdenum 3d_z² and bridging oxygen 2p orbitals (Fig. 4, right). Unrestricted EPT calculation of the Dyson orbital yielded a computed pole strength of 0.9091 (Fig. 4, center). The ionization energy is significantly different with UHF predicting 5.23 eV while the EPT corrected ionization energy is 2.22 eV. However, the Koopmans and Dyson orbitals describe the electron hole in one-electron transitions. Broken symmetry calculations can potentially give multiple solutions depending on the ionized state which correspond to different transitions. For example in $MoV O4−$ two different ²A″ solutions were obtained and it was not clear how they differed. Using NIO analysis it was determined that the higher energy solution is a one-electron transition, while the lower energy solution involves electronic excitation of an electron from the molybdenum to the vanadium d-shell. This example clearly shows the interpretive value of NIOs. Indeed, our own difficulties interpreting the Koopmans orbitals for this system served as our initial motivation for exploring new interpretive solutions.

Figure 5 illustrates the significant NIOs for the lowest energy transition as indicated by δ_elec > 0.1. From the NIOs in Fig. 5 it is clear that electron detachment involves the loss of an α electron from the vanadium 3d_z² orbital as indicated by δ_elec = − 1.00 for orbital 1α. The remaining occupation change eigenvalues indicate rearrangement of the electron density in response to the ionization and are produced as matching hole-particle pairs with δ_elec of equal size but opposite sign. Thus the two important hole-particle pairs identified in the ³A′ → ²A″ + e⁻ transition are those we label 2α → 3α and 1β → 3β. It is apparent that the 1β → 3β transition corresponds to promotion of a β electron from the molybdenum 3d_z² orbital to the 3d_z² orbital. The vanadium 3d_z² orbital is the same as the one from which the electron was ionized. The 2α → 3α hole-particle pair represents the compensating effect of the electron density partially filling the hole left by promotion of the β orbital. Orbital 2α suggests that superexchange between the d_yz orbitals on vanadium and molybdenum through p_z orbitals on bridging oxygen is disrupted and that α density is diminished on vanadium and enhanced on molybdenum. The large δ_elec eigenvalue associated with the 1β − 3β transition implies that more than just orbital relaxation is occurring, but that a whole electron is transferred from the molybdenum to the vanadium center. Thus the ³A′ → ²A″ + e⁻ transition might be recognized as a shake-up state, in which the ²A″ state differs from the ³A′ not only by an ionized electron but also by a swapped occupied-virtual orbital pair.

In this article we have described a natural orbital analysis of the electron density difference between ground and ionized electronic states that provides an improved representation for the interpretation of ionization processes at mean-field cost. This approach is equally valid whether making use of restricted (open-shell) or unrestricted DFT and HF approaches. Depending on the spin restriction framework in which the NIO analysis is carried out (unrestricted vs. restricted open-shell, for instance), electron correlation effects may be included to different extents and lead to different descriptions of the ionization process. Using the occupation change numbers it is possible to determine not only from which orbital the detached electron originated but also from the electronic transitions which describe relaxation of the density after ionization. This analysis therefore provides an enhanced understanding of electronic processes occurring in a number of applications involving ionization and provides further connections between experiment and simulation.

We gratefully acknowledge Professor J. Vincent Ortiz (Auburn University) for helpful discussions, financial support from UC Merced, and computing time on the Multi-Environment Computer for Exploration and Discovery (MERCED) cluster which is supported by National Science Foundation Grant No. ACI-1429783.