Excited states via coupled cluster theory without equation-of-motion methods: Seeking higher roots with application to doubly excited states and double core hole states

A conceptually simple approach to solving the Schrödinger equation is to diagonalize the Hamiltonian represented by the many-particle basis set spanning the entire Hilbert space. While this full configuration interaction (FCI) approach (or exact diagonalization) is formally exact, it becomes quickly unfeasible due to the exponentially growing dimension of the Hilbert space.¹ 

Coupled-cluster (CC) theory, which is usually limited to singles and doubles (i.e., CCSD), has been a popular approximate solver to the Schrödinger equation. Unlike truncated CI methods, truncated CC methods are size-consistent and therefore can be reliably applied to large systems and reach the thermodynamic limit. Most of the CC applications have been focused on approximating the ground state (GS) of systems, and therefore CC methods are usually considered to be ground state methods.² 

There is a way to compute excitation energies of CC wavefunctions based on the equation-of-motion (EOM-CC)³ formalism or the linear response (LR-CC)⁴ formalism. EOM-CCSD is one of the widely used CC excited state (ES) methods, which provides very accurate excitation energies for states dominated by single-excitations. The accuracy of EOM-CCSD for valence single excitations is about 0.1–0.2 eV. However, EOM-CCSD commonly fails to predict excitation energies for states with a significant amount of double-excitation character, and the typical error is about 1 eV or even greater than this. These failures could be avoided if the desired excited state is in a different irreducible representation from that of the ground state, since one could just employ a ground state CCSD calculation. However, if there is no point group symmetry in the system or the desired state is in the same irreducible representation, this workaround is no longer an option. The failure of EOM-CCSD for doubly excited states is largely due to the lack of relaxation of doubles amplitudes, which can be usually achieved by having triple excitations (i.e., EOM-CCSDT).^5,6 Since EOM-CCSDT has a cost, which scales $O(N8)$⁠, much research has been dedicated to improving the double-excitation gaps of EOM-CCSD by approximating the effect of connected triples either via an $O(N8)$ scaling method with a smaller prefactor or via an $O(N7)$ scaling method. These methods include EOM-CCSDT-n,^7,8 EOM-CCSD(T),⁷ EOM-CCSD(⁠$T̃$⁠),⁸ EOM-CCSD(T′),⁸ CC3,^9,10 and CCSDR(3).¹¹ 

Another challenging class of excited states for EOM-CC are core-ionized states. In most cases, these core ionization energies can be well described by the EOM with ionization potential (EOM-IP) approaches.^3,12 However, one has to obtain a large number of eigenvectors to cover the energy range for core ionizations, which can be very time-consuming for an $O(N5)$ method. In principle, one may use a modified Davidson procedure that can target high-lying interior roots.¹³ However, due to the strong coupling between core states and continuum states, numerical instability is often encountered.¹⁴ There are tricks to remedy this problem to an extent via core valence separation (CVS),^14–17 but it does not solve the inherent drawbacks of EOM-IP-CCSD. In other words, CVS-EOM-IP-CCSD fails when EOM-IP-CCSD fails. In particular, for some core-ionized states, the EOM with excitations up to doubles is not sufficient.

Recently, there has been increasing interest in so-called ΔSCF methods^18–23 as an alternative to the linear-response mean-field approaches such as CI singles (CIS) and time-dependent density functional theory (TDDFT).²⁴ In this category, the most popular approach is based on the maximum overlap method (MOM) developed by Gilbert, Besley, and Gill.²² The resulting approximate excited states from ΔSCF are not orthogonal to the approximate ground state. This seems to be suboptimal since the exact excited state should be orthogonal to the exact ground state. However, extensive benchmarks have so far suggested that the nonorthogonality of approximate wavefunction methods is not problematic to get good energetics. Furthermore, it is possible to diagonalize the Hamiltonian with those nonorthogonal determinants to obtain orthogonal states in the end. This approach is called nonorthogonal CI (NOCI).^25,26

Similar in spirit to ΔSCF, it is possible to obtain approximate solutions to exact excited states using the CC wavefunction parametrization. We call this approach ΔCC, and this is the focus of our work. In ΔCC, one computes the ground state CCSD (GS-CCSD) and an excited state CCSD (ES-CCSD) energies and takes a difference between them to compute the corresponding excitation gap. Performing an ordinary ground state CCSD calculation on an excited reference determinant leads to a desired ES-CCSD energy. Just like ΔSCF targets an excited SCF solution, ΔCC targets an excited CC solution that starts from an excited reference state. We emphasize that ΔCC is not a new approach and has been known in the literature for a while.^27–39 In particular, there are seminal works by Kowalski and co-workers that attempt to find higher roots in CC methods using the homotopy method.^31,35 They also established connections between these roots and excited states in FCI for model systems such as H₄. ΔCC has been underappreciated because of the obscure nature of CC amplitude solutions. In particular, the higher roots of the CC amplitude equation are difficult to assign to a specific state. It is also often very difficult to converge the CC amplitude equation, and multiple CC roots sometimes correspond to the same FCI state.³⁸ 

While these drawbacks make ΔCC not so appealing in general, we will show that ΔCC can be an accurate tool for excited states that are dominated by one Hartree-Fock (HF) state. CCSD with perturbative triples [CCSD(T)] is a de facto standard method for the ground state of systems with one dominant determinant. One may expect CCSD(T) to work well as long as the underlying electronic structure has only one dominant determinant, which does not need to be the ground state. In such cases, we expect the excited state CCSD(T) energies to be quite accurate and even similar in quality to that of the ground state calculation. We found excited states dominated by one double-excitation to be a perfect candidate for this approach. This is largely because the state assignment becomes much easier since it is dominated by one determinant. As mentioned earlier, EOM-CCSD fails to describe such states with dominant double-excitations so ΔCC can be an excellent alternative with the same $O(N6)$ cost.

It is also worthwhile to note that ΔCC has been used in the literature to compute core ionization energies.^40–43 Similarly to the double excitations, this is due to the ease of assigning proper states as well as relatively more stable amplitude iterations. The amplitude convergence can often become problematic, but this issue can be completely removed by a CVS-like treatment, which freezes core hole orbitals as proposed in Ref. 43. The resulting CVS-ΔCC is a good computational tool for targeting core-ionized states at the cost of ground state CCSD calculations while retaining the full flexibility of the CC wavefunction. In this work, we will focus on the computation of double ionization potentials (DIPs), which currently not many methods are able to compute. In particular, the CVS implementation of EOM-DIP-CCSD^44,45 is unavailable at the time of writing this manuscript. Furthermore, we will illustrate that EOM-DIP-CCSD does not retain the full flexibility of CCSD, and it is not an exact approach for computing electronic energies for 2-electron systems when starting from a 4-electron reference.

The goal of this paper is to (1) revive the idea of ΔCC with the emphasis on targeting doubly excited states and double core hole (DCH) states and (2) present numerical data on small molecules to support this idea.

Coupled-cluster (CC) wavefunctions use an exponential parametrization,

where $T^$ is the CC cluster operator defined as

with $τ^μ$ being the excitation operator, which creates |Φ_μ⟩ from the reference determinant |Φ₀⟩, and t_μ is the cluster amplitude. The CC ansatz then follows

assuming that |Ψ⟩ is an eigenstate of $H^$ and E_CC is the corresponding eigenvalue (i.e., energy),

The amplitudes {t_μ} are obtained by solving

where |Φ_μ⟩ denotes the μth excited state determinant,

Up to this point, we have not assumed whether we are trying to approximate the ground state or one of the excited states. In fact, the only assumption that has been made is that the state |Ψ⟩ is an eigenstate of a given Hamiltonian.

The bias toward the reference state |Φ₀⟩ built in the exponential parametrization controls which state we are targeting. The exponential parametrization is expanded to

Since nonstrongly correlated systems typically have amplitudes smaller than 1, the largest component in a usual CC wavefunction is the reference state |Φ₀⟩. This is known for model problems due to the work by Kowalski and co-workers.³¹ However, with the advances in ΔSCF methods,²² it is meaningful to revisit this idea for more complex chemical systems. We will denote such excited CC states as ES-CC states, whereas the ground state CC state will be referred to as the GS-CC state. The energy difference between the GS-CC and ES-CC states defines the ΔCC approach for electronic excitation energies. This viewpoint can also be easily extended to number-changing excitations such as ionization potential (IP) and electron attachment (EA).

It is important to note two main limitations of these ΔCC methods for electronic excited (EE) states, IP, and EA. First, traditional CC (TCC) methods are not capable of describing strong electron correlation, so ΔTCC methods are limited to states of single-reference character. Those with multireference character need more sophisticated CC approaches that can handle strong correlation. Examples of such approaches include the CC valence bond with singles and doubles (CCVB-SD),^46,47 parametrized CCSD (pCCSD),⁴⁸ and distinguishable cluster SD (DCSD),⁴⁹. For the purpose of this paper, we will focus on the application of TCC approaches to states described well by a single determinant. Applying more advanced CC approaches to multireference problems will be an interesting topic for future study. We will refer to TCC simply as CC for the rest of this paper.

Second, the computation of transition properties such as oscillator strengths and transition dipole moments is not straightforward and seems to scale exponentially with system size. Any transition properties between GS-CCSD and ES-CCSD states should technically involve a CC state for both bra and ket (first-order derivatives for each). Moreover, orbitals of GS-CCSD are not orthogonal to any orbitals of ES-CCSD in general. The evaluation of transition properties therefore formally scales exponentially with system size if done exactly. This contrasts with EOM approaches where the bra state is not a CC state, instead it is only a linear wavefunction with the same set of orbitals as the GS-CCSD state. One may consider linearizing both of the CC states to evaluate transition properties to get an approximate answer, but the exact evaluation of such properties is still highly desirable. For the purpose of this work, we will compute only energies and leave the computation of transition properties to the future study.

Furthermore, ΔCC methods often require auxiliary calculations to guide the selection of the reference determinant for a desired state. This can often be hinted by other economical approaches such as MP2 or CC2, but there is currently no generally applicable approach to this problem. We also note that ΔCC methods can often lead to artificial spatial symmetry breaking when targeting single core-hole (SCH) states in symmetric molecules. This is another manifestation of their inability to describe excited states with multireference character.

For a given ground state CC wavefunction, one can solve a Hamiltonian eigenvalue problem in the linear response space. We first define the CC Lagrangian,

where the subscript C implies that it involves only “connected” diagrams⁵⁰ and the bra is defined as

with the deexcitation operator $Λ^$ being

Evidently, we have $L=ECC$ for t such that the CC amplitude equation [Eq. (5)] is satisfied. Then, the equation-of-motion (EOM) Hamiltonian (or the CC Jacobian) can be derived from the linear response of this Lagrangian,⁴ 

where t₀ is a set of amplitudes that satisfies the “ground state” CC amplitude equation. Since $L$ is a linear function of λ, J is independent from λ. The EOM is linear response because it is a derivative of an energy expression with respect to the wavefunction parameter for both bra and ket.

In EOM-CCSD, Eq. (11) is formed in the space of singles and doubles. Evidently, EOM-CCSD cannot describe any excited states that mainly contain triples and higher excitations. What may not be immediately obvious is that EOM-CCSD, in practice, cannot describe excited states with strong double excitation character. In our view, there are two aspects of Eq. (11) that should be highlighted: (1) orbitals are determined for the ground-state SCF calculation and are fixed and (2) the CC amplitudes, t, are also determined for the ground state and are also fixed. This naturally imposes constraints on EOM calculations, and reducing the effects of those constraints requires higher-order excitations (in this case triples). This has, of course, been well-known in the community, and a method such as EOM-CC(2,3) is motivated by this observation.⁵¹ EOM-CC(2,3) takes the ground state CCSD wavefunction and forms the CC Jacobian in the space of singles, doubles, and triples. This is not really a linear-response method since it goes beyond the ground state parameter space, but it has shown to improve the accuracy of EOM-CCSD greatly, especially for states with strong double excitation character. As we will see later, without any nonperturbative connected triples, ΔCCSD and ΔCCSD(T) can perform significantly better than EOM-CCSD.

This formalism can be extended to Fock space to treat the number of electrons different from that of the ground state. The method that is relevant to the present work is the EOM ionization potential (EOM-IP) methods. In EOM-IP-CCSD, the “singles” operator (1p) removes an electron and the “doubles” operator (1h2p) removes two electrons from occupied orbitals and adds an electron to one of the unoccupied orbitals. By performing EOM-IP-CCSD on an N-electron system, one can obtain the energies of the corresponding (N − 1)-electron system and therefore ionization potentials. Removing an electron from a molecule must be accompanied by sufficient orbital relaxation. This is implicitly done by the 1h2p operator which resembles the singles operator for (N − 1) systems. Interestingly, EOM-IP-CCSD effectively has only “singles”-type excitations from a (N − 1) reference state via the 1h2p operator and no higher excitations. Therefore, the flexibility of EOM-IP-CCSD is smaller than that of CCSD or EOM-CCSD in terms of describing correlation between electrons.

A similar conclusion can be drawn for EOM double IP (EOM-DIP) methods. EOM-DIP-CCSD employs the “singles” operator (2p) which removes two electrons and the “doubles” operator (1h3p) which removes three electrons from occupieds and adds an electron back to virtuals. From an (N − 2) reference state, this EOM-CC state effectively has only “singles” excitations. A majority of those singles would account for orbital relaxation, and only little correlation effect would be gained from using EOM-DIP-CCSD.

The limited flexibility of EOM-DIP-CCSD can be most clearly understood by considering a model problem that contains four electrons and four orbitals. If we apply EOM-DIP-CCSD to this system, one would generate some determinants within the two-electron Hilbert space but not all. In Fig. 1(b), we explicitly show four determinants that are unreachable via EOM-DIP-CCSD if one uses the ground state determinant with four electrons shown in Fig. 1(a). This is somewhat disappointing because CCSD is exact for 2-electron systems. EOM-DIP-CCSD is not exact for 2-electron systems when starting from a 4-electron reference. On the other hand, if one were to compute DIPs of a 4-electron system via ΔCCSD, at least the 2-electron system energy is exactly treated via CCSD. The remaining error is then solely from the CCSD error in the ground state.

A spin-flip EOM method, EOM-SF-CCSD, can access a different spin-manifold by spin-flipping from a higher spin-manifold to a lower spin-manifold. This is most commonly used to describe diradical systems,³ but one may use it to describe doubly excited states. Not every doubly excited state can be accurately described by EOM-SF-CCSD. Generally speaking, doubly excited states which promote at least one electron to lowest unoccupied molecular orbital (LUMO) is within the scope of EOM-SF-CCSD. This is most obvious from an example for a singlet doubly excited state. One starts from a triplet ground state which typically singly occupies highest occupied molecular orbital (HOMO) and LUMO. Therefore, at the reference level, there is already an electron promoted to the LUMO, which may help to describe doubly excited states accurately via spin-flips. The same logic applies to the EOM double EA (EOM-DEA) method. In EOM-DEA-CCSD, doubly excited states that involve the promotion of two electrons to the LUMO can be within the scope of the method. Similarly to EOM-SF-CCSD, this is because one starts from a (N + 2)-electron determinant which typically doubly occupies an orbital that corresponds to the LUMO of the N-electron determinant.

For the ground state of H₂, CCSD is exact since it includes all possible excitations of two electrons in the system. Similarly, EOM-CCSD is exact for every state of H₂ and therefore with conventional ground-state CCSD (i.e., GS-CCSD) and EOM-CCSD one can get all of the electronic states of H₂ exactly for a given basis set. We will show that it is possible to reproduce those exact energies with the excited-state CCSD (i.e., ES-CCSD) method without severe numerical issues.

First, in Fig. 2, we present the results for H₂ with the STO-3G basis set. With this basis set, there are only four states in the M_S = 0 sector. All of these states can be obtained by running CCSD calculations with a carefully chosen reference determinant along with initial guess amplitudes. For the ground state and the doubly excited state $(1σu)2$⁠, the MP2 amplitude guess was used. For the singly excited states (singlet and triplet), we use a guess of $Thhll=±1$⁠, respectively, where $Thhll$ denotes the doubles amplitude for 1σ_g → 1σ_u. These reference determinants are also shown in Fig. 3. This strategy was enough to obtain the numerical data presented in Fig. 2.

The same principles can be applied to a larger basis set calculation (cc-pVTZ)⁵² as shown in Fig. 4. A distinct feature of the targeted ES-CCSD method is that it follows a root of the same character throughout the potential energy surface (PES). This is most obvious from the ES-CCSD state obtained from the $(1σg)1(2σg)1$ reference shown in Fig. 4. At R = 0.5 Å, it starts as the 6th excited state of EOM-CCSD, stays on the same state, and eventually becomes the third excited state of EOM-CCSD as the bond gets elongated. Around R = 1.25 Å, an avoided crossing appears between the third and fourth excited states. The GS-CCSD energies for these two states switch near this avoided crossing. This is natural for a targeted excited state since it follows a state of desired character.

We consider three molecules (CH₂, ethylene, and formaldehyde) and investigate their singlet dark excited states. Those excited states are characterized by promoting two electrons from HOMO to LUMO. A procedure to perform targeted ES-CCSD calculations is as follows:

1. Perform a ground state HF calculation.

2. Promote two electrons from HOMO to LUMO and converge this non-Aufbau determinant via the MOM algorithm.

3. Perform the CCSD or CCSD(T) calculations on top of the optimized non-Aufbau determinant.

In the third step, numerical instability may occur when starting from MP1 amplitudes as a guess. As a simple fix to this problem, we rescale the MP1 guess by a factor of 0.5 or less. In the future, one may use regularized MP1 amplitudes as a guess.^53–56 Furthermore, we start to update the amplitudes via the direct inversion of the iterative subspace (DIIS) algorithm^57–59 from the beginning as opposed to starting after a few iterations. With these tricks, we were able to converge the ES-CC calculations in this work without numerical difficulties.

The ground state of methylene (or carbene) is triplet. The singlet ground state (1¹A₁) for CH₂ is therefore an excited state. We optimized the geometry of CH₂ on this electronic surface with ωB97X-D and the def2-QZVPPD basis set. Interestingly, the next excited state with the same term symbol (i.e., 2¹A₁) has strong double excitation character. This doubly excited state is dominated by a closed-shell single determinant, and therefore it is a perfect candidate for the ΔCC methods. Furthermore, it is possible to perform brute-force methods such as the semistochastic heat-bath CI (SHCI) method and a second-order perturbation correction (SHCI+PT2) on this system.⁶⁰ As such, we compare ΔCCSD and EOM-CCSD against near-exact SHCI results. We employed the frozen-core approximation for the results presented in this section.

In Fig. 5, the excitation energies for the (1¹A₁ → 2¹A₁) transition are presented for various methods computed with the aug-cc-pVQZ basis set. For Δ methods, we used a reference that doubly occupies the 1¹A₁ LUMO. In other words, we used a reference with a transition of $(3a1)2→(1b2)2$⁠. As shown in Fig. 5(a), the use of this non-Aufbau state with ground state orbitals yields ill-behaved ΔSCF and ΔMP2 energies. This erratic behavior does not appear in the case of ΔCCSD and ΔCCSD(T) due to the singles operator. EOM-CCSD shows an error of 1.89 eV, which is much larger than its typical error for valence single excitations. In contrast, ΔCCSD and ΔCCSD(T) show remarkably accurate excitation energies whose errors are less than 0.1 eV. This is because the 2¹A₁ state is mainly dominated by one closed-shell determinant, which can be accurately described by CCSD.

In Fig. 5(b), we examine the effect of orbital relaxation for the Δ methods. The non-Aufbau determinant was subsequently optimized to lower the HF energy using the MOM algorithm.²² ΔSCF and ΔMP2 improve significantly when using an orbital-optimized excited state determinant. However, in the case of ΔCCSD and ΔCCSD(T), the results are more or less the same as before. This is expected because the important effect of orbital optimization can be incorporated through single excitations. In passing, we note that other variants of EOM-CCSD (SF, DEA, and DIP) may work well for this excited state. In particular, from the triplet ground state, one single spin-flip would be sufficient for EOM-SF-CCSD to access this state. States mainly dominated by one spin-flip are usually accurately described by EOM-SF-CCSD. It will be interesting to investigate these variants and assess their accuracy in the future.

With advances in brute-force approaches, it is now possible to produce high-quality benchmarks for small molecules. An example of such benchmarks is Loos and co-workers’ recent study where they used a brute-force selected CI (sCI) approach⁶¹ to produce reference energies for doubly excited states of a total of 14 small molecules: acrolein, benzene, beryllium, carbon dimer, carbon trimer, ethylene, formaldehyde, glyoxal, hexatriene, nitrosomethane, nitroxyl, pyrazine, and tetrazine.

Interestingly, most of these molecules exhibit multireference character in the doubly excited state, and they are challenging for single-reference CCSD to describe properly. For instance, Loos and co-workers considered the 2s² → 2p² excitation in Be. This state has three equally important determinants, $2px2$⁠, $2py2$⁠, and $2pz2$⁠. Therefore, this state is of multireference by nature and not easy to describe with CCSD. More complications arise for butadiene and its isoelectronic species, acrolein and glyoxal. These molecules have a low-lying dark state that has three significant configurations, one of which is a doubly excited configuration. This dark state is likewise beyond the scope of conventional single-reference CCSD.

On the other hand, the doubly excited states of ethylene and formaldehyde were found to be well described by a single determinant. Therefore, these are perfect candidates for the ΔCC approach. Given the CH₂ results discussed above, we obtained the excitation gaps for Δ methods with optimized ground-state and excited-state orbitals.

In the following benchmarks, we shall compare our results against benchmark numbers reported by Loos and co-workers.⁶¹ For smaller basis sets, they produced near-exact excitation gaps based on sCI with second-order correction (sCI+PT2) and extrapolated full CI (exFCI) methods.⁶² This should be adequate in assessing the quality of Δ methods for small basis sets. For larger basis sets, Loos and co-workers produced EOM-CC3 excitation energies. As we will see, EOM-CC3 is less accurate than ΔCCSD(T) when compared to sCI+PT2 and exFCI in smaller basis sets. Nevertheless, EOM-CC3 is a widely used iterative $O(N7)$ correlated excited state method that can yield qualitatively correct excitation gaps for doubly excited states. As such, we will also compare ΔCC methods to EOM-CC3.

In the case of ethylene, as shown in Fig. 6, EOM-CCSD significantly overestimates the gap by 2–3 eV. This highlights the failure of EOM-CCSD for doubly excited states. The doubles amplitudes, R₂, in EOM-CCSD are not enough to describe this state, and it is necessary to incorporate triple excitations to reach reasonable accuracy as for instance in EOM-CCSDT with aug-cc-pVDZ. The role of quadruples is relatively unimportant in this case. The use of a reference determinant, $(1b1u)2→(1b2g)2$⁠, yields remarkably accurate excitation energies with ΔCCSD and ΔCCSD(T). The largest T₁ and T₂ amplitudes from the ES-CCSD calculation are 0.0472 and 0.1640, respectively. These small amplitudes mean that the underlying excited state is largely dominated by the reference non-Aufbau determinant. ΔCCSD(T) excitation energies are within the error bar of exFCI. Compared to sCI+PT2, the errors are 0.01 eV and 0.03 eV for aug-cc-pVDZ and aug-cc-pVTZ basis sets, respectively. This is better in accuracy than EOM-CC3 whose error is 0.49 eV for both basis sets.

Since the doubly excited state of ethene is well described by a single determinant, relatively accurate gaps from ΔSCF are not unexpected. What is surprising is the striking underestimation of the gap in ΔMP2. One would think that for problems for which a single determinant is qualitatively correct MP2 should perform well. While this is commonly true, in the case of ΔMP2, the orbital optimization of a non-Aufbau determinant often leads to a very small gap between HOMO and LUMO.⁶³ Consequently, the MP2 correlation energy for such determinants would be heavily overestimated (i.e., more negative than it should be). As an attempt to remedy this problem, we applied a recently developed regularized MP2 method (κ-MP2).^53–56 Since small energy gaps will be damped away, the resulting correlation energy is stable even for those non-Aufbau determinants. As shown in Fig. 6, Δκ-MP2 excitation energies are similar to those of ΔCCSD, which highlights the utility of κ-MP2 for excited state simulations.

In Fig. 7, we present another successful application of ΔCC methods. The doubly excited state of formaldehyde is largely dominated by one determinant. Similarly to previous examples, EOM-CCSD significantly overestimates the excitation gap. The error of EOM-CCSD is about 4–5 eV in this case. EOM-CCSDT greatly improves but incorporating quadruples (i.e., EOM-CCSDTQ) is necessary to reach near-exact results. Directly targeting the excited state with CCSD and CCSD(T) using a non-Aufbau reference determinant [$(2b1)2→(2b2)2$] handles this state nearly exactly. The largest T₁ and T₂ amplitudes from the ES-CCSD calculation are −0.1157 and 0.1620, respectively, which implies that the underlying excited state is largely dominated by the reference non-Aufbau determinant. With aug-cc-pVTZ, ΔCCSD and ΔCCSD(T) yield an error of 0.07 eV compared to sCI+PT2. ΔCCSD overestimates, whereas ΔCCSD(T) underestimates the gap. This is better than EOM-CC3, which overestimates the gap by 0.81 eV. Given the accuracy of ΔCCSD(T), we conclude that the role of connected quadruples in describing this state can be made negligible with a properly chosen reference deterimnant. We also note that ΔSCF produces a qualitatively correct gap and ΔMP2 does not exhibit the overcorrelation problem previously shown in the case of ethylene. Therefore, Δκ-MP2 does not offer any improvement. In fact, Δκ-MP2 performs about 0.2 eV worse than ΔMP2.

In summary, not every doubly excited state requires an explicit treatment for triples unlike what was stated in Loos and co-workers’ work.⁶¹ It is only those states that are dominated by more than one determinant that require a more sophisticated treatment than single-reference CCSD and CCSD(T). For doubly excited states with one dominant determinant, we showed that CCSD and CCSD(T) can directly target such states by simply employing a non-Aufbau determinant as a reference state. The errors of ΔCCSD and ΔCCSD(T) were found to be less than 0.1 eV for the systems considered in this work.

Core-ionized states are another class of excited states that can be effectively handled by ΔCC methods. In fact, this was noted in the literature several times^40–43 and was recently revived by Zheng and Cheng.⁴³ In particular, Zheng and Cheng benchmarked single core hole (SCH) states for various small molecules and found about 0.13 eV standard deviation for ΔCCSD(T) in the ionization energies with respect to experimental values. Interested readers are referred to Ref. 43 for further information about their work.

What we will focus in this work is the use of ΔCCSD(T) for double core hole (DCH) states. Following the prescription by Zheng and Cheng for SCH states, we first obtain an (N − 2) electron reference state and freeze two unoccupied core orbitals for numerical stability. The removal of unoccupied core orbitals is similar in spirit to the CVS^{14–16,64,65} treatment in EOM-CC, and it explicitly prevents the CC wavefunction from collapsing to the ground state of the same number of particles. In our case, a double excitation from the HOMO to the unoccupied core orbitals would yield much lower energy than the desired core-ionized state. This is the source of numerical instability. The approach which freezes core hole orbitals will be referred to as CVS-ΔCC.

Investigating DCH states to probe chemical environment was first proposed by Cederbaum and co-workers.^66–68 Compared to SCH states, DCH states are much more sensitive to the chemical environment. A classic example that illustrates this point is the series of hydrocarbons, C₂H₂, C₂H₄, and C₂H₆.⁶⁶ Creating a SCH state by removing an electron from a carbon atom in these molecules results in IPs that differ only by tens of electronvolts from each other. On the other hand, DIPs exhibit a difference over 4 eV or so per C–C bond. This highlights the utility of DCH states in probing the chemical environment. Since Cederbaum’s proposal, DCH states have also been experimentally realized.^69–77 In particular, two-site DCH (TSDCH) states are sensitive to the chemical structure, so obtaining TSDCH states in experiments has become a focus.^73,75 A single-site DCH (SSDCH) state can be readily obtained from a closed-shell (N − 2) reference determinant, whereas TSDCH states are inherently of open-shell singlet character. In this section, we will study both kinds of DCH states and apply the ΔCC approach to obtain their electronic energies.

The algorithm to obtain DCH states in CVS-ΔCC is as follows:

1. Perform an SCF calculation on an N-electron system.

2. Localize core orbitals (and optionally valence orbitals separately) if there is more than one atom for the chemical element of interest. We employed Boys localization⁷⁸ for this step, and other localization schemes are also possible.

3. Identify core orbitals that will be made to be unoccupied.

4. Remove two electrons from hand-selected core orbitals, and perform an (N − 2)-electron SCF calculation using the MOM algorithm or Newton’s method.

5. Perform a CCSD calculation on the converged (N − 2)-electron reference. Note that it is necessary to freeze the two core hole orbitals in this step for numerical stability.

We note that the CVS approach naturally requires frozen-core calculations for the N-electron ground state calculation.

We will investigate the SSDCH states of five small molecules, CO, CH₄, NH₃, N₂, and CO₂, and compare the DIP values computed from CVS-ΔCCSD and CVS-ΔCCSD(T) with those of experiments. All geometries were obtained from geometry optimization with ωB97X-D⁷⁹ and aug-cc-pCVTZ.^52,80

In Table I, we present DIPs for SSDCH states using ΔCC methods with increasing the size of basis set. In the case of ionizing two electrons from a carbon atom, we observe roughly 2 eV of correlation effects in CO. On the other hand, the correction effect plays a smaller role for CH₄. In both molecules, increasing the size of the basis set reduces the DIPs. CO is well within the error bar of the experimental value, partly because the experimental error bar is quite large. For CH₄, CVS-ΔCCSD and CVS-ΔCCSD(T) exhibit an error on the order of 1 eV. Due to the lack of relativistic effect treatment in our calculations, this error is not so surprising⁸¹ and we will leave more thorough benchmarks for future study. Nonetheless, ΔSCF, CVS-ΔCCSD, and CVS-ΔCCSD(T) all yield the correct trend that CO’s DIP is several eVs higher than CH₄’s DIP. This qualitative conclusion holds even at the SCF level.

For SSDCH states involving nitrogen core vacancies, we investigated NH₃ and N₂. Their experimental estimates are about 10 eV apart. Similar to the previous cases, we observe smaller DIPs with larger basis sets. With the aug-cc-pCV5Z basis set, we observe about 1.5 eV error for all the methods. The result improves as we go from SCF to CVS-ΔCCSD(T). A major source of error is again the lack of relativistic effects. Nevertheless, all these methods successfully capture qualitative differences between these two chemical species. Specifically, the DIPs of NH₃ and N₂ differ by about 10 eV.

In the case of oxygen, we investigated CO and CO₂. The DIPs of these molecules were experimentally shown to differ by about 5 eV. This qualitative difference is well described by all of the methods examined here. However, the quantitative agreement between CVS-ΔCC methods and experimental values was only within several electronvolts as before in other systems.

In passing, we note that we neglected strong correlation present in double core hole states. Specifically, in N₂, there are two possible ways to obtain a DCH state on a nitrogen atom. Likewise, there are two equally important choices for CO₂ for generating a DCH state on an oxygen atom. One may think that this would require mixing of two such references. In fact, this effect was studied in the context of NOCI with singles for simulating core-valence excitations in our group.⁸² Given the quantitative agreement between CVS-CCSD(T) and experimental values for IPs reported in Ref. 43, we suspect that such a strong correlation effect may not be important in simulating core-ionized states. This can also be found in experimental results where single core hole states are usually localized on one atom even when there are multiple atoms of the same chemical element.⁸³ 

Unlike SSDCH, TSDCH states are all open-shell states. Two core holes are created at different atomic sites, and therefore two unpaired open-shell electrons will remain. There are two possible ways to spin-couple these two electrons: singlet and triplet. We will obtain rough energetics of these states by employing a broken-symmetry HF reference state whose ⟨S²⟩ is close to 1.0.

Such a reference state is well-suited for Yamaguchi’s approximate spin projection (AP).⁸⁴ With AP, one can obtain a spin-pure energy for the singlet state. Specifically,

where E_BS is the broken-symmetry state energy, E_S=0 is the singlet energy, E_S=1 is the triplet energy, and the coefficient α is given by

which uses the ⟨S²⟩ values of BS and S = 1 states (assuming $⟨S2⟩S=0=0$⁠). Clearly, Eq. (12) requires not only a broken-symmetry calculation but also a M_S = 1 calculation to obtain E_S=1 and $⟨S2⟩S=1$⁠. Within our CVS approach, an M_S = 1 TSDCH calculation would require a different number of frozen core and virtual orbitals for α and β spin sectors. This is uncommon to run for most quantum chemistry packages available at the moment, and therefore we will leave the use of AP for TSDCH states for future study.

The TSDCH states of four small molecules, CO, CO₂, N₂, and N₂O, are investigated. We report the DIP values computed from CVS-ΔCCSD and CVS-ΔCCSD(T) and compare them with those of experiments in Table II. All geometries were obtained from geometry optimization with ωB97X-D and aug-cc-pCVTZ.

First, we study the TSDCH states where one core hole is localized on carbon and the other one is localized on oxygen in CO and CO₂. Experimentally, these two molecules have DIPs that are 6 eV apart from each other. The difference is quite small at the SCF level as two values are only 2 eV apart. With CVS-ΔCCSD, the difference becomes 3.2 eV and CVS-ΔCCSD(T) yields a difference of 3.5 eV. While these are not quantitatively accurate, they all still correctly reproduce the qualitative behavior observed experimentally.

Next, we investigate the TSDCH states in N₂ and N₂O by creating one core hole on each nitrogen. The DIPs of these two molecules are only 2 eV apart and almost within the experimental error bar from each other. Nonetheless, our goal is to reproduce the fact that the DIP of N₂ is slightly larger than the DIP of N₂O. At the SCF level, the trend is reversed. With ΔSCF, N₂ has a DIP that is 2 eV lower than that of N₂O. With CVS-ΔCCSD and CVS-ΔCCSD(T), a correct trend is reproduced. At the CCSD level, the DIP of N₂ is only 0.3 eV higher than that of N₂O, whereas the difference becomes 0.9 eV at the CCSD(T) level.

As we can see, even without the relativistic treatment and spin-projection, we observe good qualitative agreement between CVS-ΔCC methods and experiments. It will be valuable to revisit these systems with proper relativistic corrections and Yamaguchi’s AP and try to observe a quantitative agreement between theory and experiments. We note that ⟨S²⟩ at both SCF and CCSD levels lies between 0 and 2, which asserts that these states are suitable for Yamaguchi’s AP.

In this section, we applied the ΔCC method to both single-site and two-site double core hole states. Similarly to the doubly excited states studied in Sec. IV, there was no difficulty encountered as long as the underlying CC state we are targeting is single-reference by nature. Furthermore, the satisfactory numerical stability was ensured by freezing core holes for CC calculations. The resulting CVS-ΔCC methods were tested on a variety of small molecular systems. While it was difficult to make a quantitative comparison between CVS-ΔCC and experiments due to the lack of relativistic treatment and large error bars in experimental values, both CVS-ΔCCSD and CVS-ΔCCSD(T) captured qualitative trends observed in experiments even when ΔSCF failed to do so. Furthermore, given small differences between CVS-ΔCCSD and CVS-ΔCCSD(T), it appears that the role of electron correlation can be fully captured at the CVS-CCSD level. We were not able to make comparisons to other available approaches because EOM-DIP-CCSD cannot obtain those highly excited core hole states at a reasonable cost. Furthermore, a production-level CVS-EOM-DIP-CCSD implementation is currently unavailable.

In this work, we revisited the long-standing idea of using the coupled-cluster (CC) wavefunction to directly target excited states that may be beyond the scope of equation-of-motion (EOM) approaches. In particular, we focused on using CC with singles and doubles (CCSD) to describe (1) doubly excited states and (2) double core hole states.

For doubly excited states, we show that it is possible to directly target an excited state through the ground state formalism of CCSD without numerical difficulties as long as the targeted state is dominated by one determinant. We achieve this simply by employing a non-Aufbau reference determinant that is orbital-optimized at the mean-field level via the maximum overlap method. A directly targeted CCSD and CCSD with a perturbative triples [CCSD(T)] excited state was shown to yield excellent excitation gaps for CH₂, ethylene, and formaldehyde. In particular, ΔCCSD(T) was shown to yield near-exact excitation gaps when compared to brute-force methods in a small basis set. This is quite promising since EOM-CCSD typically exhibits an error greater than 1 eV for these states. Furthermore, ΔCCSD(T) was found to be more accurate than EOM third-order approximate CC (EOM-CC3).

Likewise, double core hole states (DCHs) can be directly obtained from the ground state CCSD formalism. This is also done by using a non-Aufbau reference determinant that has double core holes. To ensure numerical stability, those core holes were frozen in correlation calculations. The resulting ΔCC ansatz is referred to as core-valence separation (CVS)-ΔCC and was benchmarked over double ionization potentials (DIPs) of small molecular systems (CO, CO₂, N₂, N₂O, and NH₃). Without relativistic corrections, CVS-ΔCCSD and CVS-ΔCCSD(T) were not able to reach quantitative accuracy when compared to experimental values. Nonetheless, they were able to estimate correct trends even when the mean-field method (ΔSCF) could not.

With the success of ΔCC described here, some interesting new directions become apparent. A more thorough investigation of open-shell singlet states in conjunction with Yamaguchi’s spin-projection will be interesting. Currently, for valence excitations, we investigated states dominated by a closed-shell determinant. Two-site DCHs were investigated without spin-projection. With spin-projection, a broader class of states will be accessible and spin-pure energies for two-site DCHs can be obtained. Second, the use of more sophisticated CC methods such as the CC valence bond with singles and doubles (CCVB-SD) for targeting excited states with a multireference character will be interesting. Finally, in addition to core-ionized states, targeting core-valence excited states will be a promising candidate to apply the techniques described in this work. Some of these are currently underway in our group.

This work was supported by the Director, Office of Science, Office of Basic Energy Sciences of the U.S. Department of Energy under Contract No. DE-ACO2-05CH11231. We thank Anna Krylov for stimulating discussions on useful applications of ΔCC for core-ionized states. J.L. thanks Soojin Lee for constant encouragement.