The accurate and efficient description of strongly correlated systems remains an important challenge for computational methods. Doubly occupied configuration interaction (DOCI), in which all electrons are paired and no correlations which break these pairs are permitted, can in many cases provide an accurate account of strong correlations, albeit at combinatorial computational cost. Recently, there has been significant interest in a method we refer to as pair coupled cluster doubles (pCCD), a variant of coupled cluster doubles in which the electrons are paired. This is simply because pCCD provides energies nearly identical to those of DOCI, but at mean-field computational cost (disregarding the cost of the two-electron integral transformation). Here, we introduce the more complete pair extended coupled cluster doubles (pECCD) approach which, like pCCD, has mean-field cost and reproduces DOCI energetically. We show that unlike pCCD, pECCD also reproduces the DOCI wave function with high accuracy. Moreover, pECCD yields sensible albeit inexact results even for attractive interactions where pCCD breaks down.

## I. INTRODUCTION

While the description of the ground state of weakly correlated systems is by now fairly routine, the same cannot be said for strongly correlated problems. Because coupled cluster theory^{1–3} offers exceptional accuracy for the description of weak dynamic correlations, we would like to use some variant of coupled cluster theory as well for the strong static correlations; in this way, one could seamlessly merge the two ideas to provide a more powerful technique which should be accurate for both regimes. Unfortunately, the construction of coupled cluster techniques for strongly correlated systems is a work in progress.

A particularly interesting recent development is the notion of pair coupled cluster doubles (pCCD),^{4–9} which takes the simple coupled cluster wave function and makes the dramatic simplification that the only allowed excitations are of a paired form in which two electrons are removed from the same spatial orbital and placed in some other spatial orbital. In other words, pCCD is coupled cluster doubles restricted to include only seniority zero determinants, where the seniority of a determinant is the number of singly occupied spatial orbitals. This restriction greatly decreases the cost of the coupled cluster calculation to mean field or $O ( M 3 ) $ if one ignores the two-electron integral transformation, but rather paradoxically, despite this simplification, the pCCD wave function yields results very close to the doubly occupied configuration interaction (DOCI),^{10–16} which includes such pair excitations to all excitation levels. This DOCI method is not new and includes many powerful geminal wave functions, including the antisymmetrized geminal power (AGP), the antisymmetric product of strongly orthogonal geminals (APSG), and many others. For many problems of interest, DOCI is able to describe the basics of the strong correlations and to the extent that pCCD reproduces DOCI, so too does pCCD.

We note, however, that the coincidence between pCCD and DOCI is not entirely universal. For the attractive pairing Hamiltonian

where *σ* indexes spins and *G* > 0, we observe that as the pairing strength *G* becomes larger and the mean-field solution develops an instability toward a number symmetry broken Hartree-Fock-Bogoliubov state, pCCD breaks down dramatically, overcorrelating wildly before eventually returning complex energies.^{17} This suggests that pCCD may not be able to describe the kinds of strong correlations needed to model phonon-mediated superconductivity, for example.

Moreover, even when the pCCD energy accurately reproduces DOCI, the wave functions may not. In coupled cluster theory, the Hamiltonian is similarity transformed, yielding a non-Hermitian effective Hamiltonian $ H \u0304 =exp ( \u2212 T ) \u2009H\u2009exp ( T ) $. Because $ H \u0304 $ is non-Hermitian, it has different left- and right-hand eigenvectors,

where *Z* creates excitations to the left and |0〉 is the mean-field reference. This structure in turn translates to different left- and right-hand wave functions for the original Hamiltonian, which in pCCD are taken to be

where *T* and *Z*, respectively, create pair excitations and pair de-excitations when acting to the right. In Ref. 9, we and others showed that the overlap

is close to unity, but this does not separately test the pCCD left-hand and right-hand wave functions and, as we will see below, while the pCCD right-hand wave function is close to the DOCI state, the same is not necessarily true of the pCCD left-hand wave function. This can cause the pCCD density matrices to be less accurate approximations of the actual DOCI density matrices than we would like.

This manuscript introduces the pair extended coupled cluster method (pECCD) which seeks to remedy these deficiencies. Where pCCD is the seniority zero version of coupled cluster doubles, pECCD is the seniority zero version of the extended coupled cluster doubles method of Arponen and Bishop.^{18–24} While extended coupled cluster has not seen a great deal of use due to its large computational cost ($O ( M 10 ) $ for extended coupled cluster doubles^{22,25}), pECCD has the same mild mean-field $O ( M 3 ) $ scaling with system size displayed by pCCD, though with a rather larger prefactor. We will sketch the pCCD method in Sec. II and provide a few results in Sec. III to show the basic performance of the method before providing our conclusions and prospects for further development in Sec. IV. We should note that the equations needed for efficient computation of the pECCD energy and wave function amplitudes are exceedingly lengthy despite their low computational cost; here, we provide a simpler but computationally more demanding expression for the energy (from which the amplitude equations follow) and have not included the amplitude equations at all. The Appendix contains the more efficient $O ( M 3 ) $ energy expression, and code to solve the amplitude equations has been made available in the supplementary material.^{26}

We should emphasize that DOCI treats only a portion of the Hilbert space of the problem, and all correlations outside of the seniority zero sector must be included in some manner if one expects to reach quantitative accuracy. That is, DOCI is valuable conceptually because it includes many strong correlation effects, but it should not be understood as a stand-alone solution to the Schrödinger equation. Accordingly, while DOCI (and therefore also pCCD and pECCD) is a useful piece of the correlation puzzle, other pieces are missing and must be in some way incorporated. Figure 1 shows total energies from restricted Hartree-Fock (RHF), DOCI, and full configuration interaction (FCI) in the 6-site Hubbard Hamiltonian and in an active space of N_{2}. Indeed, both DOCI and its non-variational pCCD approximation describe strong correlations quite accurately, though clearly residual correlations beyond seniority zero are needed for quantitative agreement with exact results.

## II. THEORY

The basics of the pair extended coupled cluster method are simple. We define a pair excitation operator,

and a pair de-excitation operator,

where the pair creation operator $ P p \u2020 $ is given by

and the pair annihilation operator *P _{p}* is its adjoint. Note that these operators are nilpotent ($ P p 2 =0$) and that more general pairing schemes than this simple singlet pairing within the same spatial orbital are possible.

^{27}

Regardless, having defined an excitation operator *T* and a de-excitation operator *Z*, the pECCD energy is given by

where in the second line, we have used the facts that the de-excitation operator *Z* annihilates the vacuum to the right and that the similarity transformed Hamiltonian $ H \u0304 = e \u2212 T \u2009H\u2009 e T $ is a six-body operator that creates up to hextuple excitations to the right, so that no more than hextuple de-excitations (created by *Z*^{3}) are needed for the case under consideration. The amplitudes *t _{ia}* and

*z*are obtained by solving

_{ai}We will take the physical Hamiltonian to be

where 〈*p*|*h*|*q*〉 and 〈*pq*|*v*|*rs*〉 are one- and two-electron integrals; note that the two-electron integrals here are not antisymmetrized. Greek letters index spin while Latin letters index spatial orbitals, with *i*, *j*, *k*, … denoting occupied orbitals, *a*, *b*, *c*, … denoting virtual orbitals, and *p*, *q*, *r*, … denoting arbitrary orbitals.

We can greatly simplify the derivation of the pECCD energy (and therefore amplitude equations) by replacing the physical Hamiltonian *H* with the portion *H*^{δΩ=0} which preserves seniority. For a two-body Hamiltonian, we would have

where *H*^{δΩ=2} and *H*^{δΩ=4} couple determinants whose seniorities differ by two and by four, respectively. Determinants differing by an odd seniority have different electron numbers, while determinants differing by an even seniority greater than four differ by a triple excitation or higher (and thus cannot be coupled by a two-body operator). Because every spatial orbital in pCCD or pECCD has seniority zero (i.e., all orbitals are either doubly occupied or empty, as singly occupied orbitals have seniority one), it will suffice for our purposes to determine what we will call $ H 0 \delta \Omega = 0 $, which is the part of *H*^{δΩ=0} that preserves the seniority of each orbital. We caution, however, that $ H 0 \delta \Omega = 0 $ omits terms which change the seniority of individual orbitals while preserving the total seniority; these terms do not contribute to closed-shell pCCD or pECCD but could contribute to a kind of ROHF-based generalization of pCCD or pECCD.

We derive our expression for $ H 0 \delta \Omega = 0 $ in the Appendix and merely quote it here,

where the necessary integrals are

and where the number and spin operators are given by

where $ \sigma \u2192 $ is the vector of Pauli matrices. The Heisenberg-like term $ \u2211 p \u2260 q K p q \u2009 S \u2192 p \u22c5 S \u2192 q $ will not contribute in our closed-shell case because closed-shell orbitals have spin zero. The number operator and pair creation and annihilation operators satisfy SU(2) commutation relationships

One can verify that

where |0〉 is a closed shell or restricted open-shell determinant (i.e., a determinant which uses the same spatial orbitals for different spins).

Given the seniority-preserving Hamiltonian, the pECCD energy can be obtained simply by computing the pECCD density matrices. As discussed in Ref. 9, the density matrices of seniority zero methods are sparse, and in fact the only non-zero elements are precisely those we need to evaluate the expectation value of $ H 0 \delta \Omega = 0 $. We may of course also use density matrices to take other expectation values; in particular, they can be used in orbital optimization.^{9} In this work, we will use orbitals optimized for pCCD, because in our experience, the pCCD- and pECCD-optimized orbitals are virtually identical and the pCCD orbital optimization is somewhat less expensive because the pCCD density matrices are simpler to compute.

To compute density matrix elements, we define a pECCD expectation value via

Given this expectation value, we find that the density matrices are given by

where $ \delta \u0304 p q =1\u2212 \delta p q $. Then, the energy is just given by

where indices *p* and *q* label arbitrary orbitals and where we have disregarded the Heisenberg Hamiltonian term as its expectation value vanishes. One can verify that the pCCD energy is properly reproduced by taking the pECCD energy and omitting terms of $O ( Z 2 ) $ or $O ( Z 3 ) $.

While constructing the density matrices given in Eq. (18) would appear to require $O ( M 6 ) $ computational time, the summation restrictions can be lifted by adding and removing terms which, with sufficient exertion, allows one to evaluate the density matrices and thus the energy in $O ( M 3 ) $ time after introducing intermediates. For example, we may write

in terms of intermediates,

Similarly, the amplitude equations can be solved in $O ( M 3 ) $ operations with the appropriate definition of intermediates. We have checked the correctness of our $O ( M 3 ) $ implementation by comparison to the explicit $O ( M 6 ) $ result for random input *T* and *Z* amplitudes and integrals *h*, *v*, and *w* and have verified our $O ( M 3 ) $ amplitude equations by comparing analytic and numerical derivatives of *E*_{pECCD} for random inputs.

We should emphasize that while the pCCD and pECCD density matrices both adopt a quasi-diagonal form, the pECCD two-particle density matrices are much more complicated. This is simply because the pECCD left-hand wave function (see below) is more sophisticated; it contains excitations to all even orders and thereby has more flexibility in fitting DOCI than does pCCD, even though pCCD and pECCD have the same number of parameters to optimize. Note that the *z* amplitudes of pCCD do not affect the energy directly, whereas those of pECCD do.

## III. RESULTS

Following Ref. 9, we will compare the pCCD, pECCD, and DOCI energies for a variety of systems, defining, for example,

and will also assess the quality of the pCCD and pECCD wave functions by evaluating

where for both pCCD and pECCD, the right-hand wave function is

while the left-hand wave functions are

for pECCD and pCCD, respectively. Recall that

where |DOCI_{k}〉 is the *k*th excited DOCI state; thus, we will have *S* ≈ 1 provided that the DOCI excited states have minimal overlap with either the left-hand or right-hand state of pECCD or pCCD. It will also prove fruitful to look in more detail; however, at the individual left- and right-hand overlaps,

where the normalization constants $ N L $ and $ N R $ are such that the left- and right-hand states are individually normalized to unity so that, for example,

These allow us to assess separately the quality of the left- and right-hand wave functions.

As Ref. 9 makes clear, with orbital optimization, pCCD is exact for two-electron singlets. The same is of course true for pECCD. Thus, both are exact for H_{2} and nearly exact for LiH, with little to distinguish the two approaches in the latter case. It will therefore be more fruitful to focus on systems with more strongly correlated electrons.

We begin, then, with chains of equally spaced hydrogen atoms. The strong correlations in these systems seem to be described reasonably well by DOCI. Figure 2 shows that while pCCD reproduces DOCI quite well energetically, pECCD does so even better. More relevantly, while both pCCD and pECCD accurately describe the right-hand wave function (i.e., we see that |DOCI〉 ≈ exp(*T*)|0〉), the left-hand wave function of pCCD is a fairly poor approximation to the DOCI state, while the left-hand wave function of pECCD is again very accurate. This should not be too surprising, as the pCCD left-hand wave function 〈0|(1 + *Z*) exp(−*T*) consists only of the reference and doubly excited determinants, while the pECCD left-hand wave function includes all the same higher excitations that DOCI adds. Indeed, that pECCD accurately reproduces the left-hand wave function of DOCI while pCCD sometimes does not seems to be a fairly general feature. Note that while an accurate left-hand wave function is not needed for accurate energies, errors in the left-hand wave function may translate into errors for properties other than the energy. In other words, while both pCCD and pECCD match the DOCI energy, we would expect only pECCD to match DOCI for arbitrary properties. We should also perhaps emphasize the smallness of the various discrepancies; the pECCD energy differs from DOCI by less than 0.001 kcal/mol per electron, and the left- and right-hand states have overlaps with DOCI differing from 1 by about 10^{−8}, implying that the coefficients of DOCI excited states in the pECCD ground state are less than 10^{−4}. For all practical purposes, pECCD reproduces DOCI exactly.

Similar results are seen in the symmetric double dissociation of H_{2}O (Fig. 3) and in the dissociation of N_{2} (Fig. 4). While the pCCD energy and right-hand wave function are close to those of DOCI, pECCD is closer yet; meanwhile, the pECCD left-hand wave function may be a much better approximation to DOCI than the pCCD left-hand wave function. This is true not just in the basis of energetically optimized orbitals, but is also true with canonical Hartree-Fock orbitals (see, e.g., Fig. 5).

We have mentioned that the close coincidence between pCCD and DOCI breaks down for the attractive pairing Hamiltonian of Eq. (1) which in the language of the SU(2) generators discussed earlier is just

Note that for *G* > 0, pairs attract one another. This is in contrast to the electronic Hamiltonian of Eq. (12), where the pairing interaction specified by the integrals *v _{pq}* > 0 is purely repulsive. Because the pairing Hamiltonian contains only the SU(2) pairing (or pseudospin) generators, it is solved exactly by DOCI. Conveniently, however, a more compact, exact solution was found by Richardson,

^{28–30}which permits the generation of exact energies for systems far too large to be practicably solved by DOCI. In Fig. 6, we consider a 40-site pairing Hamiltonian with equally spaced levels (

*ϵ*=

_{p}*p*) at half-filling. One can see that near

*G*≈ 0.2, pCCD begins to deviate significantly from DOCI, and for

*G*≳ 0.3, we can find no real solution to the pCCD equations (solutions with complex

*T*amplitudes and complex correlation energies exist, but are of limited physical interest). While pECCD also begins to break down somewhat, it provides a much more accurate description of the correlations in the pairing Hamiltonian for large

*G*. Note that for the half-filled forty-site Hamiltonian under consideration, a broken number symmetry mean field appears at

*G*≈ 0.22, not coincidentally close to the value at which pCCD begins to break down. While pECCD is not a panacea, it provides results much superior to pCCD for the pairing Hamiltonian and competitive with pCCD based on the broken-symmetry mean-field

_{c}^{17}for

*G*not too large. To the extent that non-pair correlations renormalize the interactions in the physical Hamiltonian via exp(−

*T*

_{Ω≠0})

*H*exp(

*T*

_{Ω≠0}), where

*T*

_{Ω≠0}is the non-pair portion of

*T*, the fact that pECCD much more accurately describes attractive interactions than does pCCD may be important in the modeling of realistic materials.

As discussed earlier, because the pECCD left-hand wave function has much better overlap with the DOCI wave function than does the left-hand wave function of pCCD, one would expect pECCD to yield density matrices closer to the DOCI density matrices than does pCCD. This is indeed the case. Figure 7 shows fractional errors in the entries of the two-particle density matrix contributions $\u3008 P p \u2020 \u2009 P q \u3009$ for N_{2} in the recoupling regime (*R*_{N−N} ≈ 2.2 Å). Explicitly, we are plotting the fractional error 1 − Γ/Γ_{DOCI} in the matrix elements of the two-particle density matrix for all elements Γ_{DOCI} larger than 10^{−4}, since very small elements of Γ_{DOCI} are unlikely to have significant contributions to expectation values. We see that typically, the pCCD values differ from those of DOCI by ∼1%, while the errors in the pECCD density matrix elements are an order of magnitude or more smaller. Note that we have Hermitized the pCCD and pECCD density matrices to simplify the comparison; this has a much larger effect on pCCD than on pECCD because the pECCD density matrix is more nearly Hermitian. For example, the Frobenius norms $ Tr ( \Gamma \u2009 \Gamma T ) $ of the antisymmetric parts of two-particle density matrices discussed in Fig. 7 are 5.15 × 10^{−2} for pCCD and 7.68 × 10^{−5} for pECCD.

## IV. CONCLUSIONS

In many cases, a configuration interaction restricted to paired excitations but not restricted by excitation level can provide an accurate accounting for static correlation effects, particularly once the orbitals used to define the pairing and the reference determinant are optimized. The idea that this doubly occupied configuration interaction could describe many forms of static correlation is not a new one; indeed, DOCI itself was introduced over forty years ago.^{10–13} But while pair-excited configuration interaction attracted a great deal of early interest, the method was basically abandoned due to its exponential scaling with system size, which renders DOCI unsuitable for practical calculations.

With pair coupled cluster doubles and pair extended coupled cluster doubles, we now possess two models in the coupled cluster family which generally provide results almost indistinguishable from those of DOCI but with mean-field computational cost. While the pCCD energy and right-hand wave function typically reproduce DOCI very well, the left-hand wave function, being a linear combination of the reference and doubly excited determinants, is not always of particularly good quality. Moreover, for attractive interactions (or at least for the attractive pairing interaction), pCCD does not accurately reproduce DOCI once the mean-field is unstable toward number symmetry breaking (though we emphasize again that for the physical electronic Hamiltonian, the pairing-type interaction is repulsive and can only become attractive upon renormalization with some Hamiltonian transformation). Both these difficulties can be ameliorated by using pair extended coupled cluster doubles, which seems to offer much more accurate results for attractive interactions where it is within a few percent of DOCI for weakly broken number symmetry, as well as a much superior left-hand wave function (and therefore superior density matrices and expectation values). The computational scaling of pECCD is the same as pCCD, namely, $O ( M 3 ) $, though the pECCD energy and amplitude equations are significantly more complicated and the method is accordingly roughly an order of magnitude more expensive than is pCCD itself. Nonetheless, the mean-field computational scaling permits routine pECCD calculations on systems with hundreds of basis functions, for which we can reliably anticipate getting results of essentially DOCI quality for both the energy and for other observables. We hope, then, that pECCD will be a valuable tool for the description of strongly correlated systems, particularly when pairing interactions become attractive and number symmetry is not strongly broken, and that it will form an excellent starting point from which to add correlations which break pairs.

There are, of course, important drawbacks of pECCD as well. The most significant is that pECCD does not include dynamic correlation. One can attempt to fix this with the addition of a simple density functional correlation energy,^{31} or by relaxing the restriction to seniority zero, whether by freezing amplitudes^{8,9} or by including all possible singlet pairing terms,^{27} or potentially by including a perturbative account for higher seniority sectors. In other words, the restriction to the seniority zero sector of Hilbert space allows us to accurately reproduce DOCI, but DOCI is only a portion of the full configuration interaction, and the correlations from the remaining sectors of Hilbert space must still be incorporated to yield accurate results in general systems. Moreover, the restriction to seniority zero is sometimes too severe even for the description of strong correlations; some systems simply require higher seniority sectors for the accurate treatment of static correlation effects, as is seen, for example, in the dissociation of N_{2}.^{16} We may hope to treat these problems by lifting the restriction that the pairing scheme pairs the two spinorbitals corresponding to the same spatial orbital. This might even open the door to the use of broken spin symmetry in seniority zero methods for strongly correlated systems.

## Acknowledgments

This work was supported by the National Science Foundation (No. CHE-1462434). G.E.S. is a Welch Foundation chair (C-0036). We thank Jorge Dukelsky for helpful discussions.

### APPENDIX A: EFFICIENT ENERGY EXPRESSION

The density matrices given in Eq. (18) are evaluated in $O ( M 6 ) $ computational effort due to the presence of restrictions on summation indices. As we have noted, one can carefully add zero to the equations in such a way as to eliminate the summation restrictions, which then permits an efficient evaluation of the density matrices (and therefore the energy). Here, we quote our final result for the energy, which can be obtained by a straightforward but tedious series of manipulations,

The energy, then, can be evaluated with $O ( M 2 ) $ cost in terms of the various intermediate quantities which appear in the foregoing equation.

The intermediates we have defined are

These intermediates can be evaluated efficiently in $O ( M 3 ) $ time through a combination of matrix multiplication $ A \u22c5 B p q = \u2211 r A p r \u2009 B r q $ and element-by-element multiplication [(** A** ×

**)**

*B*_{pq}=

*A*

_{pq}*B*]. It is of course possible to further simplify the energy expression by defining more complicated intermediates.

_{pq}### APPENDIX B: USEFUL DERIVATIONS

Here, we sketch a few derivations which the reader may find helpful but which are not necessary for understanding the thrust of the manuscript.

#### 1. Derivation of $ H 0 \delta \Omega = 0 $

Recall that $ H 0 \delta \Omega = 0 $ is the portion of the Hamiltonian which preserves seniorities of individual spatial orbitals. Here, we wish to derive the expression quoted in Eq. (12).

Let us start, then, with the one-electron part of the Hamiltonian. In order to preserve orbital seniority, when we remove an electron from orbital *p*, we must return it to orbital *p*. Thus, we have

where the spinorbital number operators are

and their sum is the spatial orbital number operator *N _{p}* given in Eq. (14a).

The two-electron part of the Hamiltonian is slightly more complicated. We could remove two electrons from orbital *p* and place them in orbital *q*, or we could remove one electron from orbital *p* and another from orbital *q*≠*p*, in which case, we must place electrons back in orbitals *p* and *q*; in the latter case, we must include the possibility of an exchange where the first electron is removed from *p* but placed in *q*. All told, we have

where the spin index $ \sigma \u0304 $ is the opposite of the index *σ* and where we have made use of the pair creation operator $ P p \u2020 $ given in Eq. (7) and the adjoint operator *P _{q}* and have introduced the spin raising and lowering operators

To simplify further, we note that

where

and *σ*^{z} is a Pauli matrix. This means that

whence

in terms of the spin vector operators defined in Eq. (14b) so that overall, we may make the replacement

as desired.

#### 2. Density matrices

In order to compute the pECCD energy, we need to evaluate density matrix elements. In other words, we need to take the expectation value

where we recall that |0〉 is the single-determinant reference. As noted in Ref. 9, the density matrices of seniority zero methods are sparse, and in fact the only non-zero entries of the full one- and two-particle density matrices can be constructed from 〈*N _{p}*〉

_{pECCD}, 〈

*N*

_{p}*N*〉

_{q}_{pECCD}, and $\u3008 P p \u2020 \u2009 P q \u3009 pECCD $.

Computing the density matrix elements is simplified by constructing similarity transformations of the number operators *N _{p}* and the pair creation and annihilation operators $ P p \u2020 $ and

*P*. We find

_{q}where we have taken advantage of the commutator expansion

and have used the SU(2) commutation relationships given in Eq. (15) and transcribed here for convenience,

We have added a few summation index restrictions which are unnecessary (but permitted) in view of the nilpotency of the pair creation and annihilation operators, simply because they may help clarify the origins of the summation restrictions and factors such as $ \delta \u0304 a b $ and $ \delta \u0304 i j $ in the final density matrices of Eq. (18).