The information contained within ground-state one- and two-electron reduced-density matrices (RDMs) can be used to compute wave functions and energies for electronically excited states through the extended random phase approximation (ERPA). The ERPA is an appealing framework for describing excitations out of states obtained via the variational optimization of the two-electron RDM (2-RDM), as the variational 2-RDM (v2RDM) approach itself can only be used to describe the lowest-energy state of a given spin symmetry. The utility of the ERPA for predicting near-edge features relevant to x-ray absorption spectroscopy is assessed for the case that the 2-RDM is obtained from a ground-state v2RDM-driven complete active space self-consistent field (CASSCF) computation. A class of killer conditions for the CASSCF-specific ERPA excitation operator is derived, and it is demonstrated that a reliable description of core-level excitations requires an excitation operator that fulfills these conditions; the core-valence separation (CVS) scheme yields such an operator. Absolute excitation energies evaluated within the CASSCF/CVS-ERPA framework are slightly more accurate than those obtained from the usual random phase approximation (RPA), but the CVS-ERPA is not more accurate than RPA for predicting the relative positions of near-edge features. Nonetheless, CVS-ERPA is established as a reasonable starting point for the treatment of core-level excitations using variationally optimized 2-RDMs.

## I. INTRODUCTION

The development of powerful x-ray sources and sophisticated measurement technology over the years has enhanced the utility of near-edge x-ray absorption fine structure (NEXAFS)^{1} as a tool for probing local chemical information. Near-edge features contain a wealth of information about the electronic structure, providing insight into the oxidation state, coordination number, and covalency,^{2,3} while time-resolved NEXAFS can even be used to follow reaction intermediates with elemental specificity.^{4} The interpretation of experimental data from NEXAFS is often aided through the use of quantum-chemical computations. Many electronic structure methods have been applied to this problem, including time-dependent density functional theory (TDDFT),^{2,3,5–9} linear response^{10–14} and equation-of-motion^{15–18} coupled-cluster theory, second-order algebraic diagrammatic construction [ADC(2)],^{19–22} and active-space-based multiconfigurational approaches.^{23–30} While TDDFT is often the method of choice, owing to its relatively low computational cost, quantitative agreement with experimental data can require a sophisticated treatment of both static and dynamical correlation effects, as well as relativistic effects.^{31} Unfortunately, standard approaches for capturing static correlation, such as the complete active space self-consistent field (CASSCF) method,^{32–35} are associated with steep computational costs that limit their application to small systems.

The limitations of configuration-interaction (CI)-based CASSCF stem from the exponential complexity of the active-space wave function. In contrast, the number of unknowns in the active-space two-electron reduced-density matrix (2-RDM) grows only quartically with the size of the active space. For this reason, variational 2-RDM (v2RDM)-driven CASSCF^{36,37} algorithms can be applied to much larger active spaces than can be considered with CI-based CASSCF. For example, most CI-based algorithms are limited to active spaces comprised of at most 18 electrons in 18 orbitals, although active spaces as large as 22 electrons in 22 orbitals have recently been realized using massively parallel algorithms.^{38} On the other hand, an active space as large as 50 electrons in 50 orbitals can routinely be considered with a v2RDM-driven algorithm. Unfortunately, despite its nice scaling properties, v2RDM-CASSCF is not directly applicable to the NEXAFS problem, as v2RDM-based methods can only describe ground electronic states of a given spin symmetry. In light of this issue, several strategies have been proposed to extract excited-state information from ground-state reduced-density matrices (RDMs), including the Hermitian operator method^{39–43} explicitly time-dependent RDM methods,^{44,45} and the extended random phase approximation (ERPA).^{46–48} The present work focuses on the application of ERPA to the NEXAFS problem.

The ERPA is a specific application of Rowe’s equation of motion (EOM)^{49,50} to the case in which excited-state wave functions are assumed to be well described by single excitations out of the ground state. The ERPA leads to a nonsymmetric generalized eigenvalue problem that can be formulated in terms of the ground-state one-electron reduced-density matrix (1-RDM) and the 2-RDM, making it an ideal candidate for computing excitation energies from variationally obtained RDMs. While Rowe’s equation of motion is itself formally exact, the accuracy of excitation energies computed from its application within the ERPA will depend on (i) how well ground state 1- and 2-RDMs resemble those of the exact ground state and (ii) how well the actual excited-state wave functions are represented by single excitations out of the ground state.

In this work, we explore the utility of the ERPA for describing core-level excitations out of ground states described at the v2RDM-CASSCF level of theory. In Sec. II, we present the working equations for the ERPA and discuss the computational cost associated with solving the ERPA equations when using a CASSCF reference. We also derive a set of killer conditions that must be considered for a reliable description of core-level excitations. These conditions can be satisfied by removing all excitations that do not involve a core orbital [i.e., within the core-valence separation (CVS) approximation^{17,51–54}]. Section III then provides the details of our computations. In Sec. IV, we assess the quality of the core-level spectra obtained from the ERPA and examine the effects of electron correlation, *N*-representability, core-valence separation, and the choice of orbital basis (i.e., canonical versus energy-optimized orbitals) on the excitation energies.

## II. THEORY

Following Refs. 46 and 47, the ERPA is easily derived as a specific case of Rowe’s equation of motion.^{49} First, consider the ground electronic state, |Ψ_{0}⟩, and an excited state, |Ψ_{n}⟩ (*n* > 0), both of which are eigenfunctions of the Hamiltonian

with eigenvalues *E*_{0} and *E*_{n}, respectively. We introduce an excitation operator, $\xd4n\u2020$, that defines the excited state as

and a deexcitation operator, *Ô*_{n}, the action of which on the ground state yields zero,

where *ω*_{n} = *E*_{n} − *E*_{0}, and subsequently left-multiply Eq. (5) by ⟨Ψ_{0}|*Â* to obtain

Here, ⟨Ψ_{0}|*Â* represents an arbitrary state within the manifold of states defined by ⟨Ψ_{0}|*Ô*_{n}. If the killer condition is satisfied, Eq. (6) can be reexpressed in the more convenient form

which is Rowe’s formally exact equation of motion.

Within the ERPA, the excitation operator is restricted to include only single excitations out of the ground state. We choose the singlet spin-adapted excitation operator

where $\xe2q\u2020$ and *â*_{p} represent fermionic creation and annihilation operators, respectively, the indices *p* and *q* are spatial orbital labels (specifically, the natural orbitals of the ground state), and an overbar (or lack thereof) represents an electron of *β* (or *α*) spin. Note that the summation includes excitations between all orbitals; this structure contrasts with that of the more familiar random phase approximation (RPA), which involves only particle/hole and hole/particle transitions. Furthermore, the excitation operator excludes terms where *p* = *q*, which is a consequence of the killer condition.^{46} Specifically, it is the consequence of the following constraint implied by, but weaker than, Eq. (4):

In the natural orbital basis, Eq. (9) can be satisfied by restricting the labels in Eq. (8) so that *p* ≠ *q*.^{46}

To obtain the working equations of the ERPA, the arbitrary excitation operator, *Â*, is chosen to be

The double commutator on the left-hand side (LHS) of Eq. (11) is a two-particle operator, the expectation value of which is expressible as a sum over the elements of the ground state 2-RDM, and the right-hand side depends on only the natural orbital occupation numbers ($nr=\u27e8\Psi 0|\xe2r\u2020\xe2r\u2020|\Psi 0\u27e9=\u27e8\Psi 0|\xe2r\xaf\u2020\xe2r\xaf\u2020|\Psi 0\u27e9$). In this form, the construction of the LHS of Eq. (11) should scale as the sixth power of the size of the one-electron basis set, and the solution of the generalized eigenvalue equation for all states shares the same scaling.

The excitation operator given by Eq. (8) is the most general operator for the ERPA problem. In practice, the actual form of the operator depends upon the choice of the ground-state wave function. For example, in the case that the ground-state is a restricted Hartree-Fock wave function, the ERPA reduces to the usual RPA, with $\xd4n\u2020$ comprised of only particle/hole and hole/particle transitions. The present work is concerned with the case that the ground-state wave function is obtained from CASSCF, so we partition our orbitals into restricted (doubly occupied), active (partially occupied), and external (empty) orbitals, and we note that excitations among restricted or external orbitals cannot contribute to any excited-state wave function. An appropriate excitation operator for CASSCF-specific ERPA is then

where the orbital labels *i* and *a* correspond to restricted and external orbitals, respectively, and the labels *t* and *u* correspond to orbitals belonging to the active space. The application of ERPA with this excitation operator can be thought of as an uncontracted multireference configuration interaction method limited to single-electron transitions, with the exception that the ERPA also includes deexcitations. With this operator, the evaluation of the LHS of Eq. (11) and the solution of the generalized eigenvalue problem still scale as the sixth power with system size, but the prefactor is significantly reduced, relative to the general case. The dimension of the matrix representation of the LHS of Eq. (11) is 2*k*_{r}*k*_{a} + 2*k*_{r}*k*_{e} + 2*k*_{a}*k*_{e} + *k*_{a}(*k*_{a} − 1), where *k*_{r}, *k*_{a}, and *k*_{e} represent the number of restricted, active, and external orbitals, respectively; the solution of the generalized eigenvalue problem scales as the cube of this number.

As will be discussed in Sec. IV, the last term in Eq. (12), which involves active/active-type excitations, can lead to severe violations in the killer condition [Eq. (4)] that degrade the quality of the ERPA excitation energies. We therefore consider an additional class of killer conditions beyond Eq. (9) that can be derived by the projection of Eq. (4) onto an arbitrary state $\Psi 0|\xe2r\u2020\xe2s$,

These constraints, combined with Eq. (9), are stronger than Eq. (9) alone, but they are still weaker than the killer condition itself. In the case that the ground-state wave function has the form of a CASSCF wave function, we consider nine separate constraints defined by the space to which the operators $\xe2r\u2020$ and *â*_{s} belong, four of which are not trivially satisfied. The constraints

imply

where $Duwtv2$ and $Duw\xaftv\xaf2$ represent elements of the active-space 2-RDM, defined as

Additionally, the constraints

and

lead to

and

Equation (21) is identical to the condition that would be obtained from the same considerations in the case that the reference wave function was a Hartree-Fock wave function (i.e., within the RPA). The only way to satisfy this condition is to set all restricted/external deexcitation weights to zero; if these terms are eliminated from the excited-state operator expansion for RPA, that method reduces to the Tamm-Dancoff approximation (TDA) or CI with single excitations (CIS). Equations (22) and (23) are generalizations of this condition for CASSCF reference functions that can only be satisfied by eliminating the relevant deexcitations from the operator expansion. These conditions both reduce to Eq. (21) in the limit of zero electron correlation within the active space (when *n*_{t} → 0 or *n*_{t} → 1) and can thus be classified as TDA-like restrictions on deexcitations.

In this work, we ignore the TDA-like restrictions of Eqs. (21)–(23) and focus on Eq. (15). As will be shown in Sec. IV, we observe large violations of Eq. (15) for some K-edge features in small molecules. An ERPA procedure that explicitly enforces these conditions could improve the description of these states, but we find that reasonable results can be obtained simply by invoking the CVS approximation, removing from the ERPA expansion all transitions that do not involve a core orbital. Assuming that all valence electrons reside within the active space, the CVS-ERPA excitation operator can be defined as

In addition to automatically satisfying the conditions given by Eq. (15) [and Eq. (23)], the CVS approximation prevents accidental degeneracies between core-excited states and valence excitations that lie far above the ionization threshold.

Before moving on, we should highlight two important properties of Eq. (11). First, the excitation energies obtained from this equation are size intensive. We have verified this property numerically for the system of non-interacting LiH molecules described in Ref. 55. For these tests, the ground-state was treated at the full-valence v2RDM-driven CASSCF level of theory, and the excitation operator employed within the ERPA was that given in Eq. (12). Second, for a stable ground-state reference function, Eq. (11) can be recast as a half-sized nonsymmetric eigenvalue problem with real, positive eigenvalues. In the case that the orbitals are optimized for the reference state, the stability condition is satisfied. However, should one employ a different orbital basis, the stability condition could be violated. In principle, the lack of symmetry in the LHS of Eq. (11) could then lead to complex eigenvalues, as well as to matrix defects relevant to the description of conical intersections.^{56,57}

## III. COMPUTATIONAL DETAILS

In Sec. IV, we compare ERPA excitation energies to those obtained from the usual RPA, TDDFT (with the B3LYP functional), and time-dependent equation of motion (EOM) second-order approximate coupled cluster theory (TD-EOM-CC2).^{58} RPA and TDDFT excitation energies were computed using the General Atomic and Molecular Electronic Structure System (GAMESS) software package,^{59} and TD-EOM-CC2 computations were performed using a plugin to the Psi4 package.^{60} TD-EOM-CC2 computations employed Cholesky-decomposed electron repulsion integral (ERI) tensors with a tight decomposition threshold of 10^{−12} E_{h}. The molecular geometries employed within all computations were optimized for the ground state using density functional theory (with the B3LYP functional^{61}) and the cc-pVQZ basis set.^{62} All other ground and excited state computations were performed using the cc-pVTZ, cc-pCVTZ,^{63} aug-cc-pVTZ,^{64} and cc-pVQZ basis sets.

The ground-state 1- and 2-RDMs entering the ERPA equations were obtained from v2RDM-driven CASSCF computations, as implemented in a plugin^{37} to Psi4. Table I provides the active spaces employed within these computations. Optimized RDMs satisfy either two-particle (PQG) *N*-representability conditions^{65} or two-particle and partial three-particle (T2) conditions.^{66,67} All v2RDM-CASSCF computations employed the density fitting approximation to the ERI tensor. Computations within the cc-pVTZ, aug-cc-pVTZ, and cc-pVQZ basis sets used the corresponding jk-type fitting basis,^{68} while computations within the cc-pCVTZ basis used the def2-QZVpp-jkfit fitting basis.^{69} The DF approximation reduces the scaling of the orbital optimization step in v2RDM-CASSCF from $O(k5)$ to $O(k4)$, where *k* represents the dimension of the one-electron basis. All ERPA computations also employed the DF approximation, which reduces the formal scaling of the construction of some terms that contribute to the matrix representation of the double commutator on the LHS of Eq. (11) from $O(k6)$ to $O(k5)$. Even within the DF approximation, though, the overall scaling of the construction of this matrix and of the solution to the ERPA eigenvalue problem is unchanged from that described in Sec. II.

. | Point . | Restricted . | Active orbitals . | ||
---|---|---|---|---|---|

Molecule . | group . | orbitals . | f.v.^{a}
. | f.v.+3s^{b}
. | f.v.+3s3p^{c}
. |

CO | $C2v$ | [2, 0, 0, 0] | [4, 0, 2, 2] | [6, 0, 2, 2] | [8, 0, 4, 4] |

H_{2}CO | $C2v$ | [2, 0, 0, 0] | [5, 0, 2, 3] | [7, 0, 2, 3] | [9, 0, 4, 5] |

HCN | $C2v$ | [2, 0, 0, 0] | [5, 0, 2, 2] | [7, 0, 2, 2] | [9, 0, 4, 4] |

N_{2}O | $C2v$ | [3, 0, 0, 0] | [6, 0, 3, 3] | [9, 0, 3, 3] | [10, 0, 7, 7] |

CH_{4} | $C2v$ | [1, 0, 0, 0] | [4, 0, 2, 2] | [5, 0, 2, 2] | [6, 0, 3, 3] |

C_{2}H_{2} | D_{2h} | [1, 0, 0, 0, 0, 1, 0, 0] | [3, 0, 1, 1, 0, 3, 1, 1] | [4, 0, 1, 1, 0, 4, 1, 1] | [5, 0, 2, 2, 0, 5, 2, 2] |

N_{2} | D_{2h} | [1, 0, 0, 0, 0, 1, 0, 0] | [2, 0, 1, 1, 0, 2, 1, 1] | [3, 0, 1, 1, 0, 3, 1, 1] | [4, 0, 2, 2, 0, 4, 2, 2] |

. | Point . | Restricted . | Active orbitals . | ||
---|---|---|---|---|---|

Molecule . | group . | orbitals . | f.v.^{a}
. | f.v.+3s^{b}
. | f.v.+3s3p^{c}
. |

CO | $C2v$ | [2, 0, 0, 0] | [4, 0, 2, 2] | [6, 0, 2, 2] | [8, 0, 4, 4] |

H_{2}CO | $C2v$ | [2, 0, 0, 0] | [5, 0, 2, 3] | [7, 0, 2, 3] | [9, 0, 4, 5] |

HCN | $C2v$ | [2, 0, 0, 0] | [5, 0, 2, 2] | [7, 0, 2, 2] | [9, 0, 4, 4] |

N_{2}O | $C2v$ | [3, 0, 0, 0] | [6, 0, 3, 3] | [9, 0, 3, 3] | [10, 0, 7, 7] |

CH_{4} | $C2v$ | [1, 0, 0, 0] | [4, 0, 2, 2] | [5, 0, 2, 2] | [6, 0, 3, 3] |

C_{2}H_{2} | D_{2h} | [1, 0, 0, 0, 0, 1, 0, 0] | [3, 0, 1, 1, 0, 3, 1, 1] | [4, 0, 1, 1, 0, 4, 1, 1] | [5, 0, 2, 2, 0, 5, 2, 2] |

N_{2} | D_{2h} | [1, 0, 0, 0, 0, 1, 0, 0] | [2, 0, 1, 1, 0, 2, 1, 1] | [3, 0, 1, 1, 0, 3, 1, 1] | [4, 0, 2, 2, 0, 4, 2, 2] |

^{a}

The full valence space.

^{b}

f.v. plus orbitals of the symmetry of the 3s orbital on each heavy atom.

^{c}

f.v. plus orbitals of the symmetry of the 3s and 3p orbitals on each heavy atom.

## IV. RESULTS AND DISCUSSION

### A. Principal K-edge features

Experimentally obtained K-edge features corresponding to 1*s* → *π ^{*}*-type excitations are provided in Table II for a set of small molecules containing

*π*bonds (CO, H

_{2}CO, HCN, C

_{2}H

_{2}, N

_{2}, and N

_{2}O). Also provided are errors in excitation energies obtained from TD-EOM-CC2, TDDFT, RPA, ERPA, and CVS-ERPA; ERPA computations employed RDMs generated from full-valence v2RDM-CASSCF. All computed results reported here were obtained using the cc-pCVTZ basis. We first consider the results obtained with the excitation operator of Eq. (12) and discuss the role of the CVS approximation below.

. | . | . | . | . | . | ERPA . | CVS-ERPA . | ||
---|---|---|---|---|---|---|---|---|---|

Molecule . | Transition . | Expt. . | CC2 . | TDDFT . | RPA . | PQG . | PQG+T2 . | PQG . | PQG+T2 . |

CO | C 1s → π ^{*} | 287.4^{a} | 2.2 | −11.3 | 7.0 | 5.6 | 5.6 | 5.2 | 5.0 |

O 1s → π ^{*} | 533.6^{b} | 2.0 | −13.8 | 16.5 | 8.3 | 8.8 | 6.9 | 6.9 | |

H_{2}CO | C 1s → π ^{*} | 285.6^{c} | 3.1 | −10.4 | 8.8 | 7.2 | 7.2 | 7.0 | 7.2 |

O 1s → π ^{*} | 530.8^{c} | 2.0 | −14.1 | 15.1 | 7.7 | 8.1 | 6.6 | 6.7 | |

C_{2}H_{2} | C 1s → π ^{*} | 285.8^{d} | 2.7 | −10.5 | 10.2 | 6.9 | 6.9 | 6.8 | 6.9 |

HCN | C 1s → π ^{*} | 286.4^{e} | 2.7 | −10.6 | 9.6 | 6.7 | 6.6 | 6.5 | 6.6 |

N 1s → π ^{*} | 399.7^{e} | 2.4 | −12.1 | 12.3 | 7.4 | 7.5 | 7.0 | 7.0 | |

N_{2} | N 1s → π ^{*} | 401.0^{f} | 2.3 | −12.5 | 11.2 | 6.8 | 6.9 | 6.4 | 6.4 |

N_{2}O | N_{t} 1s → π ^{*} | 401.1^{g} | 2.5 | −12.2 | 12.4 | 7.1 | 7.4 | 6.6 | 6.9 |

N_{c} 1s → π ^{*} | 404.7^{g} | 2.8 | −12.4 | 11.4 | 6.6 | 6.9 | 6.6 | 6.9 | |

O 1s → π ^{*} | 534.8^{g} | 1.4 | −14.3 | 18.5 | 9.3 | 9.9 | 8.2 | 8.8 | |

Mean unsigned error | 2.0 | 12.2 | 12.1 | 7.2 | 7.4 | 6.7 | 6.9 | ||

Maximum unsigned error | 2.7 | 14.3 | 18.5 | 9.3 | 9.9 | 8.2 | 8.8 |

. | . | . | . | . | . | ERPA . | CVS-ERPA . | ||
---|---|---|---|---|---|---|---|---|---|

Molecule . | Transition . | Expt. . | CC2 . | TDDFT . | RPA . | PQG . | PQG+T2 . | PQG . | PQG+T2 . |

CO | C 1s → π ^{*} | 287.4^{a} | 2.2 | −11.3 | 7.0 | 5.6 | 5.6 | 5.2 | 5.0 |

O 1s → π ^{*} | 533.6^{b} | 2.0 | −13.8 | 16.5 | 8.3 | 8.8 | 6.9 | 6.9 | |

H_{2}CO | C 1s → π ^{*} | 285.6^{c} | 3.1 | −10.4 | 8.8 | 7.2 | 7.2 | 7.0 | 7.2 |

O 1s → π ^{*} | 530.8^{c} | 2.0 | −14.1 | 15.1 | 7.7 | 8.1 | 6.6 | 6.7 | |

C_{2}H_{2} | C 1s → π ^{*} | 285.8^{d} | 2.7 | −10.5 | 10.2 | 6.9 | 6.9 | 6.8 | 6.9 |

HCN | C 1s → π ^{*} | 286.4^{e} | 2.7 | −10.6 | 9.6 | 6.7 | 6.6 | 6.5 | 6.6 |

N 1s → π ^{*} | 399.7^{e} | 2.4 | −12.1 | 12.3 | 7.4 | 7.5 | 7.0 | 7.0 | |

N_{2} | N 1s → π ^{*} | 401.0^{f} | 2.3 | −12.5 | 11.2 | 6.8 | 6.9 | 6.4 | 6.4 |

N_{2}O | N_{t} 1s → π ^{*} | 401.1^{g} | 2.5 | −12.2 | 12.4 | 7.1 | 7.4 | 6.6 | 6.9 |

N_{c} 1s → π ^{*} | 404.7^{g} | 2.8 | −12.4 | 11.4 | 6.6 | 6.9 | 6.6 | 6.9 | |

O 1s → π ^{*} | 534.8^{g} | 1.4 | −14.3 | 18.5 | 9.3 | 9.9 | 8.2 | 8.8 | |

Mean unsigned error | 2.0 | 12.2 | 12.1 | 7.2 | 7.4 | 6.7 | 6.9 | ||

Maximum unsigned error | 2.7 | 14.3 | 18.5 | 9.3 | 9.9 | 8.2 | 8.8 |

^{a}

Experimental result is taken from Ref. 70.

^{b}

Experimental result is taken from Ref. 71.

^{c}

Experimental result is taken from Ref. 72.

^{d}

Experimental result is taken from Ref. 73.

^{e}

Experimental result is taken from Ref. 74.

^{f}

Experimental result is taken from Ref. 75.

^{g}

Experimental result is taken from Ref. 76.

We find that RPA and TDDFT consistently overestimate and underestimate, respectively, the position of the K-edge features considered; the mean unsigned error in the excitation energies computed using these approaches are 12.1 eV for RPA and 12.2 eV for TDDFT. The maximum unsigned error obtained with RPA (18.5 eV) is slightly larger than that of TDDFT (14.3 eV). The ERPA improves upon the RPA for all molecules. The mean unsigned errors are reduced to 7.2 eV and 7.4 eV when enforcing the PQG or PQG+T2 *N*-representability conditions, respectively, while the maximum errors are reduced to only 9.3 eV (PQG) and 9.9 eV (PQG+T2). Interestingly, the ERPA excitation energies appear to be insensitive to the degree to which the ground-state RDMs are *N*-representable, as evidenced by the similar performance of the approach when enforcing two- or three-particle *N*-representability conditions in the underlying v2RDM-CASSCF computations. Aside from the modest improvement in the magnitudes of the excitation energies, the ERPA also appears to provide a more balanced description of multiple K-edges within a given molecule than is obtained with RPA or TDDFT. For example, in H_{2}CO, the errors in the carbon and oxygen 1*s* → *π ^{*}* excitation energies predicted by RPA are 8.8 eV and 15.1 eV, respectively, while those predicted by ERPA with PQG constraints are 7.2 eV and 7.7 eV. The basis set dependence of the 1

*s*→

*π*excitation energies is quite small, and the conclusions drawn from the data in Table II are unchanged when considering the cc-pVTZ, cc-pVQZ, and aug-cc-pVTZ basis sets. The errors in the 1

^{*}*s*→

*π*excitation energies computed in these basis sets can be found in the supplementary material.

^{*}We note that TD-EOM-CC2 consistently outperforms all other methods considered here; the mean unsigned error in the TD-EOM-CC2 energies is only 2.0 eV in the cc-pCVTZ basis set. Even greater accuracy could be obtained from a higher level of coupled cluster (CC) theory [e.g., CC with single and double excitations (CCSD)]. EOM and linear-response CCSD typically yield K-edge features for second-row atoms that are accurate to roughly 1–2 eV,^{12,18} and the inclusion of triple excitations can reduce errors to less than 1 eV.^{12}

The performance of the ERPA appears promising, but some problematic cases reveal a shortcoming of the method that is related to the CASSCF-specific excitation operator defined in Eq. (12); these cases can be identified by violations in the killer condition. For example, consider N_{2}O described by the cc-pCVTZ basis set. ERPA computations with RDMs that satisfy PQG conditions predict two features whose wave functions contain significant 1*s* → *π ^{*}* contributions at similar energies (544.1 eV and 542.1 eV), which is problematic for two reasons. First, the spectrum is only expected to have one such feature. Second, the dominant term in both expansions corresponds to a

*σ*→

*π*-type excitation; this transition is the one between orbitals that are active in the full-valence CASSCF reference function, with similar occupation numbers (

*n*

_{σ}= 0.988 and

*n*

_{π}= 0.985). It is entirely possible that this transition contributes to the true wave function, but it is unlikely that it would be the dominant term in the expansion. The large expansion coefficients [Eq. (8)] associated with this transition indicate a breakdown in the ERPA; the expansion coefficients are approximately 16 and 7 for the states at 542.1 eV and 544.1 eV, respectively.

The breakdown of the ERPA in this case may stem from the near-degeneracy of the occupations of the orbitals in question. Regardless of its precise cause, the breakdown is easily identifiable by large violations in the killer condition. We can quantify the degree to which the killer condition is violated with the root mean squared (RMS) error in Eq. (15) for all possible operators, $\xe2t\u2020\xe2u$. For the two problematic states at the oxygen K-edge in N_{2}O, the RMS errors in Eq. (15) are 1.9 × 10^{−1} (at 542.0 eV) and 8.6 × 10^{−3} (at 544.1 eV). These errors are substantially larger than those observed in other molecules; the RMS errors in Eq. (15) at the oxygen K-edge in CO and H_{2}CO are 6.5 × 10^{−6} and 1.1 × 10^{−4}, respectively. Unlike the wave functions for K-edge states in N_{2}O, those in CO and H_{2}CO are dominated by core excitations, as expected. We therefore conclude that the difficulty in assigning the edge in N_{2}O reflects violations in the killer condition and a breakdown in the ERPA, perhaps caused by the near-degeneracy of the occupations of the active orbitals.

Table II includes errors in excitation energies computed using CVS-ERPA, which should reduce killer-condition violations by excluding all active/active-type excitations. The CVS approximation yields improved excitation energies for all molecules, reducing the error from ERPA by as much as 1.8 eV in the case of the oxygen K-edge in CO (with RDMs that satisfy the PQG+T2 conditions). The mean and maximum errors for CVS-ERPA with RDMs that satisfy two-particle *N*-representability conditions are reduced to only 6.7 eV and 8.2 eV, respectively, while those for CVS-ERPA with RDMs that satisfy partial three-particle conditions are reduced to 6.9 eV and 8.8 eV. We again find that the CVS-ERPA provides a more consistent description of multiple features within a given molecule than that afforded by other methods. Consider, for example, the 1*s* → *π ^{*}* transitions of the central and terminal nitrogen atoms in N

_{2}O (denoted N

_{t}and N

_{c}in Table II, respectively). Using the CVS-ERPA and RDMs that satisfy the PQG conditions, the error in both of the excitation energies is 6.6 eV, meaning that CVS-ERPA exactly predicts the separation between the features. No other method considered here predicts a separation between these features that agrees with experiment to this degree. Finally, and perhaps most importantly, by disregarding all non-core excitations, the K-edge can be more clearly identified in all molecules, including in the problematic case of the oxygen K-edge of N

_{2}O discussed above. The consideration of the killer condition is thus extremely important in the framework of the ERPA, and all remaining ERPA-derived results discussed in this work will use the CVS-ERPA excitation operator given by Eq. (24).

### B. Rydberg features

We now consider Rydberg-type K-edge features in the small molecules considered above, as well as in CH_{4}. Table III provides experimentally-obtained excitation energies and errors in computed excitation energies for 1*s* → 3*s* features in all molecules and the 1*s* → 3*p* feature in methane. Again, all ERPA computations employed RDMs generated from full-valence v2RDM-CASSCF. As measured by the mean and maximum unsigned errors, TDDFT performs slightly better for this set of excitations than it does for the 1*s* → *π ^{*}* transitions; the mean and maximum unsigned errors for TDDFT are 10.4 and 13.1 eV, respectively. On the other hand, RPA and CVS-ERPA perform considerably worse. The mean and maximum unsigned errors for RPA are 17.4 eV and 25.5 eV, respectively. The CVS-ERPA improves over RPA by almost 5 eV for the mean unsigned error [12.5 eV (PQG) and 12.6 eV (PQG+T2)], and CVS-ERPA reduces the maximum error observed to 20.1 (20.2) eV when the procedure employs RDMs that satisfy the PQG (PQG+T2) conditions.

. | . | . | . | . | CVS-ERPA . | |
---|---|---|---|---|---|---|

Molecule . | Transition . | Expt. . | TDDFT . | RPA . | PQG . | PQG+T2 . |

CO | C 1s → 3s | 292.4^{a} | −6.1 | 14.9 | 12.7 | 12.9 |

O 1s → 3s | 538.9^{b} | −12.8 | 25.5 | 17.2 | 17.5 | |

H_{2}CO | C 1s → 3s | 290.2^{c} | −10.9 | 12.2 | 10.2 | 10.3 |

O 1s → 3s | 535.4^{c} | −12.1 | 21.2 | 15.0 | 15.3 | |

C_{2}H_{2} | C 1s → 3s | 287.7^{d} | −10.4 | 12.6 | 10.4 | 10.0 |

HCN | C 1s → 3s | 289.1^{e} | −11.4 | 11.3 | 9.3 | 8.8 |

N 1s → 3s | 401.8^{e} | −11.6 | 19.1 | 13.9 | 15.2 | |

N_{2} | N 1s → 3s | 406.2^{f} | −4.5 | 24.4 | 20.1 | 20.2 |

N_{2}O | N_{t} 1s → 3s | 404.0^{g} | −10.8 | 17.8 | 12.7 | 13.2 |

N_{c} 1s → 3s | 407.5^{g} | −11.4 | 16.9 | 12.0 | 11.1 | |

O 1s → 3s | 536.7^{g} | −13.1 | 24.7 | 10.1 | 10.2 | |

CH_{4} | C 1s → 3s | 287.1^{d} | −10.9 | 13.0 | 10.8 | 9.8 |

C 1s → 3p | 288.0^{d} | −10.2 | 11.9 | 8.5 | 8.6 | |

Mean unsigned error | 10.5 | 17.4 | 12.5 | 12.6 | ||

Maximum unsigned error | 13.1 | 25.5 | 20.1 | 20.2 |

. | . | . | . | . | CVS-ERPA . | |
---|---|---|---|---|---|---|

Molecule . | Transition . | Expt. . | TDDFT . | RPA . | PQG . | PQG+T2 . |

CO | C 1s → 3s | 292.4^{a} | −6.1 | 14.9 | 12.7 | 12.9 |

O 1s → 3s | 538.9^{b} | −12.8 | 25.5 | 17.2 | 17.5 | |

H_{2}CO | C 1s → 3s | 290.2^{c} | −10.9 | 12.2 | 10.2 | 10.3 |

O 1s → 3s | 535.4^{c} | −12.1 | 21.2 | 15.0 | 15.3 | |

C_{2}H_{2} | C 1s → 3s | 287.7^{d} | −10.4 | 12.6 | 10.4 | 10.0 |

HCN | C 1s → 3s | 289.1^{e} | −11.4 | 11.3 | 9.3 | 8.8 |

N 1s → 3s | 401.8^{e} | −11.6 | 19.1 | 13.9 | 15.2 | |

N_{2} | N 1s → 3s | 406.2^{f} | −4.5 | 24.4 | 20.1 | 20.2 |

N_{2}O | N_{t} 1s → 3s | 404.0^{g} | −10.8 | 17.8 | 12.7 | 13.2 |

N_{c} 1s → 3s | 407.5^{g} | −11.4 | 16.9 | 12.0 | 11.1 | |

O 1s → 3s | 536.7^{g} | −13.1 | 24.7 | 10.1 | 10.2 | |

CH_{4} | C 1s → 3s | 287.1^{d} | −10.9 | 13.0 | 10.8 | 9.8 |

C 1s → 3p | 288.0^{d} | −10.2 | 11.9 | 8.5 | 8.6 | |

Mean unsigned error | 10.5 | 17.4 | 12.5 | 12.6 | ||

Maximum unsigned error | 13.1 | 25.5 | 20.1 | 20.2 |

^{a}

Experimental result is taken from Ref. 70.

^{b}

Experimental result is taken from Ref. 71.

^{c}

Experimental result is taken from Ref. 72.

^{d}

Experimental result is taken from Ref. 73.

^{e}

Experimental result is taken from Ref. 74.

^{f}

Experimental result is taken from Ref. 75.

^{g}

Experimental result is taken from Ref. 76.

It is unclear exactly why the performance of the CVS-ERPA is so different for Rydberg-type features than for 1*s* → *π ^{*}* transitions. What is clear, however, is that this difference directly impacts the utility of ERPA for evaluating the relative positions of these features at the K-edge. Indeed, while CVS-ERPA yields slightly better absolute excitation energies than those provided by RPA, the overall spectral shape predicted by CVS-ERPA and RPA are similar. As a measure of the shape of the absorption curve, we present in Table IV the errors in computed gaps between the 1

*s*→

*π*and Rydberg features tabulated above, as compared to those from experiment. According to this metric, neither RPA nor CVS-ERPA perform particularly well, with slightly worse results obtained from CVS-ERPA. Mean errors for both methods are on the order of 5–6 eV, depending on the choice of the basis set. On the other hand, TDDFT provides much more reasonable results; the B3LYP functional yields the gaps between these features with an accuracy of roughly 1–2 eV, on average.

^{*}. | Basis . | . | . | CVS-ERPA^{a}
. | ||
---|---|---|---|---|---|---|

. | set . | TDDFT . | RPA . | f.v.^{b}
. | f.v.+3s^{c}
. | f.v.+3s3p^{d}
. |

Mean | cc-pVTZ | 1.8 (12.1) | 5.5 (12.1) | 6.4 (6.7) | 6.7 (7.3) | 6.1 (8.9) |

cc-pCVTZ | 2.0 (12.2) | 6.2 (12.1) | 6.3 (6.7) | 6.7 (7.3) | 6.0 (8.9) | |

aug-cc-pVTZ | 1.4 (12.5) | 4.2 (12.1) | 5.9 (6.8) | 5.7 (7.5) | 4.9 (8.9) | |

cc-pVQZ | 0.9 (12.3) | 5.5 (12.1) | 6.4 (6.7) | 6.4 (7.3) | 5.7 (9.0) | |

Maximum | cc-pVTZ | 9.6 (14.3) | 10.1 (18.5) | 14.1 (8.2) | 11.6 (8.3) | 11.2 (12.4) |

cc-pCVTZ | 8.0 (14.3) | 13.1 (18.5) | 13.7 (8.2) | 11.5 (8.4) | 10.7 (12.4) | |

aug-cc-pVTZ | 3.4 (15.3) | 6.1 (18.5) | 10.9 (8.3) | 9.7 (9.3) | 8.4 (12.3) | |

cc-pVQZ | 2.7 (14.4) | 7.9 (18.4) | 12.4 (8.2) | 10.8 (8.4) | 9.9 (12.4) |

. | Basis . | . | . | CVS-ERPA^{a}
. | ||
---|---|---|---|---|---|---|

. | set . | TDDFT . | RPA . | f.v.^{b}
. | f.v.+3s^{c}
. | f.v.+3s3p^{d}
. |

Mean | cc-pVTZ | 1.8 (12.1) | 5.5 (12.1) | 6.4 (6.7) | 6.7 (7.3) | 6.1 (8.9) |

cc-pCVTZ | 2.0 (12.2) | 6.2 (12.1) | 6.3 (6.7) | 6.7 (7.3) | 6.0 (8.9) | |

aug-cc-pVTZ | 1.4 (12.5) | 4.2 (12.1) | 5.9 (6.8) | 5.7 (7.5) | 4.9 (8.9) | |

cc-pVQZ | 0.9 (12.3) | 5.5 (12.1) | 6.4 (6.7) | 6.4 (7.3) | 5.7 (9.0) | |

Maximum | cc-pVTZ | 9.6 (14.3) | 10.1 (18.5) | 14.1 (8.2) | 11.6 (8.3) | 11.2 (12.4) |

cc-pCVTZ | 8.0 (14.3) | 13.1 (18.5) | 13.7 (8.2) | 11.5 (8.4) | 10.7 (12.4) | |

aug-cc-pVTZ | 3.4 (15.3) | 6.1 (18.5) | 10.9 (8.3) | 9.7 (9.3) | 8.4 (12.3) | |

cc-pVQZ | 2.7 (14.4) | 7.9 (18.4) | 12.4 (8.2) | 10.8 (8.4) | 9.9 (12.4) |

^{a}

CVS-ERPA computations employed RDMs that satisfied the PQG *N*-representability conditions.

^{b}

A full valence active space.

^{c}

f.v. plus orbitals of the symmetry of the 3s orbital on each heavy atom.

^{d}

f.v. plus orbitals of the symmetry of the 3s and 3p orbitals on each heavy atom.

One factor that potentially contributes to the poor performance of CVS-ERPA observed here is that the full-valence active space does not include the 3*s* and 3*p* orbitals. Table IV also provides errors in the spacing between the 1*s* → *π ^{*}* and Rydberg features obtained from CVS-ERPA with two slightly larger active spaces consisting of the full valence space plus orbitals of the symmetry of the 3s or 3s and 3p orbitals on each heavy atom (see Table I for details). The errors in the excitation energies for 1

*s*→

*π*transitions are also provided in parentheses. We find that the CVS-ERPA predictions for the relative positions of the K-edge features improve slightly with the larger active spaces, and the accuracy of CVS-ERPA approaches that of the usual RPA with the largest active space. However, the quality of the description of the 1

^{*}*s*→

*π*transitions decreases with increasing size of the active space size.

^{*}### C. The role of the orbital basis

The data discussed to this point demonstrate that the CVS-ERPA and RPA approaches offer comparable descriptions of K-edge features in small molecules. CVS-ERPA provides a modest improvement in absolute excitation energies for both 1*s* → *π ^{*}* and Rydberg-type excitations, but the methods predict similar relative positions of the excitations. Two features of the CVS-ERPA framework distinguish it from RPA, which lead to the small differences in the numerical performance of the approaches. Both distinctions are rooted in the underlying v2RDM-CASSCF-based description of the ground state. First, the RDMs that enter CVS-ERPA contain nondynamical correlation information lacking in RPA based on a Hartree-Fock reference. Second, the orbitals employed within CVS-ERPA are optimized for the correlated ground state. We can differentiate between the effects of correlation and orbital optimization within the CVS-ERPA by employing RDMs generated from ground-state active-space v2RDM computations that do not involve any orbital optimization [i.e., v2RDM-driven complete active space CI (CASCI)]. Figure 1 illustrates the mean unsigned errors for RPA and CVS-ERPA in several basis sets, and CVS-ERPA computations either employed v2RDM-CASSCF or v2RDM-CASCI references [labeled CVS-ERPA(CASSCF) and CVS-ERPA(CASCI), respectively]. In all basis sets, orbital optimization leads to a significant improvement in the quality of the CVS-ERPA excitation energies. Interestingly, in the case of the aug-cc-pVTZ basis set, it appears that an ERPA that accounts for static correlation alone offers essentially no advantage over the usual RPA.

The fact that the performance of CVS-ERPA is coupled to the choice of orbital basis is not surprising; it is well known that orbital relaxation effects are considerable for core-excitated states. In principle, the performance of CVS-ERPA could be improved by choosing a more appropriate set of orbitals than those obtained from Hartree-Fock or v2RDM-CASSCF. From the perspective of CASSCF, the most obvious strategy would be to employ a state-averaged (SA)^{77,78} v2RDM-CASSCF procedure that optimizes the orbitals for the ground and target excited states simultaneously. Excited-state RDMs required for such a procedure are accessible within the equation of motion framework. Alternatively, one could pre-relax the orbitals as is done within the static exchange (STEX) approximation^{79–81} or employ CIS natural orbitals as approximation to SA-CASSCF.^{82} While each of these strategies may yield improved orbitals for the description of core excitations, within the context of the ERPA, they all share the same potential shortcoming. As mentioned above, the ERPA problem is nonsymmetric, but it can be recast as a symmetric problem half-sized problem with real eigenvalues, provided that the reference is stable. Unfortunately, a guarantee of stability requires that the orbitals be optimized for the ground state, and the use of STEX or (approximate) SA-CASSCF orbitals could lead to instabilities in the ERPA equations.

### D. Oscillator strengths

Finally, we consider the quality of CVS-ERPA-derived oscillator strengths. Within Rowe’s equation of motion formalism and the ERPA, the transition dipole moment that defines the oscillator strength is given by

where $\mu ^\xi $ is the *ξ*-component of the dipole operator (*ξ* ∈ *x*, *y*, *z*). Figure 2 illustrates the oscillator strengths computed at the TDDFT, RPA, and CVS-ERPA levels of theory within the cc-pCVTZ basis set; CVS-ERPA computations employed RDMs from full-valence v2RDM-CASSCF that satisfied either PQG or PQG+T2 *N*-representability conditions. For 1*s* → *π ^{*}* excitations [Fig. 2(a)], CVS-ERPA oscillator strengths resemble those from TDDFT more than those from RPA. In particular, when using RDMs that satisfy PQG

*N*-representability conditions, CVS-ERPA oscillator strengths tend to agree with those from TDDFT to a greater degree than when using RDMs that satisfy PQG and T2 conditions. On the other hand, oscillator strengths for Rydberg-type excitations [Fig. 2(b)] computed using CVS-ERPA appear to resemble those from RPA more than those from TDDFT.

## V. CONCLUSIONS

We have developed a procedure for computing core-level excitation energies using a combination of the variational optimization of the ground-state 2-RDM and the extended random phase approximation. Because we lack *N*-representability conditions that differentiate between ground- and excited-state RDMs, the direct optimization of an excited-state 2-RDM is still an open problem. As a result, several non-variational approaches have been proposed that can extract excited-state information from variationally obtained ground-state RDMs. The ERPA is one such method.

We have benchmarked the quality of ERPA- and CVS-ERPA-derived 1*s* → *π ^{*}*, 1

*s*→ 3

*s*, and 1

*s*→ 3

*p*excitations in a set of small molecules for the case that the ground-state 2-RDMs are obtained from a v2RDM-driven CASSCF procedure. In general, we find that CVS-ERPA offers a modest improvement over its uncorrelated limit, the RPA, when considering absolute excitation energies. In particular, the description of 1

*s*→

*π*-type transitions that the CVS-ERPA offers is somewhat better than that provided by both RPA and TDDFT (with the B3LYP functional). However, the CVS-ERPA appears to be a less promising approach for the description of Rydberg-type excitations. In this case, CVS-ERPA absolute excitation energies are still slightly better than those from RPA, but mean and maximum unsigned errors from CVS-ERPA are larger than those from TDDFT. Furthermore, CVS-ERPA is slightly less accurate than RPA for predicting the relative positions of the principal and Rydberg features. Similar conclusions can be drawn from the inspection of oscillator strengths. For 1

^{*}*s*→

*π*-type excitations, CVS-ERPA-derived oscillator strengths are in better agreement with those from TDDFT than with those from RPA (presumably, the TDDFT oscillator strengths are the more reliable ones). However, the CVS-ERPA does not consistently offer any improvement over RPA for oscillator strengths corresponding to 1

^{*}*s*→ 3

*s*and 1

*s*→ 3

*p*excitations.

Interestingly, we have shown that the (CVS-)ERPA appears to be insensitive to the degree to which the ground-state 2-RDMs are *N*-representable. In terms of both excitation energies and oscillator strengths, results from the (CVS-)ERPA are quite similar when utilizing RDMs that satisfy either the two-particle (PQG) constraints or two-particle plus partial three-particle (PQG+T2) constraints. This result suggests that the imposition of computationally demanding higher-order *N*-representability conditions on RDMs that enter the (CVS-)ERPA may not be as crucial as it is for the ground-state problem itself.

We derived new conditions implied by the killer condition that should be fulfilled by ERPA-derived excited-state wave functions. One subset of these conditions can be satisfied by removing active/active-type excitations from the ERPA excitation operator. We showed that the CVS-ERPA excitation operator (which removes all active/active-type excitations) consistently yields excitation energies that are superior to those obtained with the original CASSCF-specific ERPA excitation operator. The truncated operator also simplifies the assignment of K-edge features in some molecules.

Finally, we demonstrated that the EOM-CC2 framework provides a description of 1*s* → *π ^{*}* excitations that is clearly superior to that from TDDFT, RPA, or CVS-ERPA, and we note that even more accurate results could be obtained by considering higher levels of CC theory. Along these lines, the success of the CC hierarchy suggests a clear route to improving CVS-ERPA. The inclusion of double excitations within the ERPA will undoubtedly improve the description of the core-excited states. Without any additional approximations, an ERPA procedure that includes doubles could be devised with knowledge of the ground-state four-electron reduced-density matrix. Furthermore, the present formulation of the ERPA could be improved by accounting for orbital relaxation effects, either through the STEX approximation or a combined SA-CASSCF/ERPA procedure.

## SUPPLEMENTARY MATERIAL

See supplementary material for errors in 1*s* → *π ^{*}* excitation energies computed at the TD-EOM-CC2, TDDFT, RPA, and ERPA levels of theory within the cc-pVTZ, cc-pCVTZ, aug-cc-pVTZ, and cc-pVQZ basis sets.

## ACKNOWLEDGMENTS

This material is based upon work supported by the National Science Foundation under Grant No. CHE-1554354.