A compact orbital representation of ionization processes is described utilizing the difference of calculated one-particle density matrices. Natural orbital analysis involving this difference density matrix simplifies interpretation of electronic detachment processes and allows differentiation between one-electron transitions and shake-up/shake-off transitions, in which one-electron processes are accompanied by excitation of a second electron into the virtual orbital space.

## I. INTRODUCTION

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^{14} 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,^{15} although in DFT calculations this value is highly dependent on the functional used and often significantly incorrect.^{16} Δ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)^{17} and the attachment-detachment density analysis of Head-Gordon *et al.*^{18} 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}.^{19} 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.^{20} 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^{21} 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.*^{23} 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,^{48} 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.

## II. THEORETICAL METHODS

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 **Δ**_{P}**S** 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, $\Delta P\alpha $, and *β* 1PDM difference matrix, $\Delta P\beta $, 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.^{51} 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*, which determines the likelihood of transition from orbital

_{n}*n*, according to

Within the diagonal self-energy approximation,^{57} 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.

## III. NUMERICAL TESTS

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^{58} 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.

### A. Formaldehyde

The first system that we examine is electron detachment of the ^{1}A_{1} ground state (GS) of formaldehyde (CH_{2}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 ^{1}A_{1} state of formaldehyde. In both restricted and unrestricted cases Koopmans’ theory predicts the electron to be ejected from a b_{2} 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 ^{2}B_{1} state. Attempts to converge to the ^{2}B_{2} state using restricted open-shell HF were unsuccessful, and in Fig. 1 we show the hole of the ^{1}A_{1} → ^{2}B_{1} + e^{−} transition computed with NIO analysis. The ^{2}B_{1} state is determined to result from electronic detachment from the CO *π* orbital. The ^{2}B_{2} 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 ^{1}A_{1} → ^{2}B_{2} + 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.

### B. Nitromethyl radical

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 ^{2}A″ → ^{1}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 $\u3008S\u02c62\u3009$ 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_{2} and NO_{2} groups in the Restricted Open-Shell Hartree-Fock (ROHF) calculation, to between the H_{2}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_{2} group, while in the unrestricted calculation the ionized electron is located within the *π* orbital delocalized over the H_{2}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

_{2}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 ^{2}A″ and ^{1}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_{2}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_{2} 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_{2}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,^{20} 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.

### C. MoV O_{4} anion

The final example of the use of NIOs that we describe in this report is that of the ^{3}A′ → ^{2}A″ + e^{−} electron detachment in $MoNbO4\u2212$. 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.^{5} 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.^{5} The UB3PW91 HOMO shows that *α* electron density is removed from 3*d*_{z2} orbitals on both vanadium and molybdenum (Fig. 4, left). The UHF HOMO shows a similar picture with greater density taken from molybdenum 3*d*_{z2} and bridging oxygen 2*p* 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\u2212$ two different ^{2}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 3*d*_{z2} 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 ^{3}A′ → ^{2}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 3*d*_{z2} orbital to the 3*d*_{z2} orbital. The vanadium 3*d*_{z2} 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*orbitals on bridging oxygen is disrupted and that

_{z}*α*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

^{3}A′ →

^{2}A″ + e

^{−}transition might be recognized as a shake-up state, in which the

^{2}A″ state differs from the

^{3}A′ not only by an ionized electron but also by a swapped occupied-virtual orbital pair.

## IV. CONCLUSION

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.

## Acknowledgments

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.