Communication: Random-phase approximation excitation energies from approximate equation-of-motion coupled-cluster doubles Free

The random-phase approximation (RPA) plays a foundational role in quantum chemistry, condensed-matter physics, materials science, and nuclear physics.^1–3 As a theory of the ground-state correlation energy, the particle-hole RPA (ph-RPA) is an infinite-order resummation of all time-independent ring diagrams, which critically controls the leading-order divergence in the energy of metals at high density.^2,4,5 Especially when combined with density functional theory via the adiabatic connection fluctuation-dissipation theorem,^6,7 the ph-RPA also provides a good description of long-range dispersion interactions.^3,8–11 Similarly, the particle-particle RPA (pp-RPA) is an infinite resummation of particle-particle and hole-hole ladder diagrams, which are important at low density,^12,13 and provides a promising treatment of the exchange-correlation energy in density functional theory.¹⁴ Both flavors of the RPA correspond to approximate bosonic treatments of harmonic electronic fluctuations.^1,15

The terms appearing in the ph-RPA and pp-RPA correlation energy are a subset of those included in coupled-cluster theory with double excitations (CCD). Therefore, approximate solutions of the CCD equations, known as ring-CCD and ladder-CCD, can be used to calculate the ph-RPA and pp-RPA correlation energy, respectively. For the ph-RPA, this was shown numerically by Freeman for the electron gas¹⁶ and proven analytically by Scuseria, Henderson, and Sorensen;¹⁷ see also Ref. 18. The equivalence between pp-RPA and ladder-CCD was proven simultaneously in Refs. 15 and 19.

Alternatively, the RPA may be viewed as a theory of dynamical response functions. The ph-RPA approximates the dynamical polarizability, a context in which it is known to be identical to time-dependent Hartree or Hartree-Fock.²⁰ For finite systems, such as molecules, the ph-RPA leads to reasonably accurate electronic excitations²¹ and underlies the successful time-dependent density functional theory.^22–26 For solids, the ph-RPA polarizability correctly predicts the properties of the collective plasmon excitation² and forms the basis for screening in the popular GW approximation.²⁷ Analogous to the correlation energy, the ph-RPA polarizability is a resummation of all time-dependent ring diagrams. Less studied, the pp-RPA approximates the dynamical response function associated with a pairing field, which does not conserve particle number and probes excitation energies associated with a system of N ± 2 electrons. Excitation energies calculated with pp-RPA have been investigated in recent years by Yang and co-workers^28,29 and shown to provide promising accuracy for states with double-excitation character.

These observations suggest a relation between excitation energies calculated with the RPA and those calculated with approximate versions of excited-state coupled-cluster theory. In this manuscript, I provide the precise recipe for this analogy, showing that the ph-RPA excitation energies (with or without exchange) can be obtained from an approximation to electronic-excitation equation-of-motion CCD (EE-EOM-CCD), and the pp-RPA excitation energies can be obtained from an approximation to double-electron-affinity (DEA-) or double-ionization-potential (DIP-) EOM-CCD. The former relation was partially addressed in early work by Emrich,³⁰ and the latter one is completely unexplored to the best of my knowledge.

Both the ph-RPA and pp-RPA can be used to approximate the two-particle Green’s function, for different choices of the time ordering. The two-particle Green’s function is defined by¹³ 

where ψ^† and ψ are creation and annihilation field operators, respectively; Ψ₀ is the N-electron ground state, T is the time-ordering operator, and (1) = (r₁, t₁) is a combined space and time variable; throughout, I only consider fixed particle number and zero temperature. The two-particle Green’s function describes a number of physical processes associated with addition or removal of two particles, the specific details of which depend on the choice of relative time ordering. The dynamical polarizability is one particular time ordering corresponding to the density-density response function,^2,13

where δn(1) = ψ^†(1)ψ(1) − ⟨ψ^†(1)ψ(1)⟩ is the density fluctuation away from the ground-state density and G is the one-particle Green’s function. In the frequency domain, the poles of the polarizability occur at all N-electron excitation energies Ω_ν, with residues given by the square of the transition densities |⟨Ψ₀|n(r)|Ψ_ν⟩|.² 

In the usual diagrammatic route,^12,13 the RPA polarizability is expressed as $Π=Π0+Π0v¯Π$⁠, where $v¯$ is the antisymmetrized Coulomb kernel and Π⁰ is the irreducible polarizability of a noninteracting particle-hole pair Π⁰(1, 2) = −iG⁰(1, 2)G⁰(2, 1). This generates the conventional form of the RPA polarizability as a sum of all time-dependent ring diagrams and particle-hole ladder diagrams (the exchange interaction leading to particle-hole ladder diagrams can also be considered a simple vertex correction in Π). The noninteracting one-particle Green’s function is defined in terms of the solution of some mean-field problem,

where here and throughout the indices, i, j, k, l are used to denote occupied orbitals and a, b, c, d are used to denote unoccupied orbitals. The locations of the poles of the RPA polarizability, i.e., the excitation energies, are the eigenvalues of the well-known ph-RPA matrix, given in Subsection II B.

An alternative time ordering of G₂ describes the simultaneous addition or removal of two electrons,

which is clearly related to the system’s response to a non-number-conserving pairing field, V = ψψ + ψ^†ψ^†. The poles of K occur at (N ± 2)-electron excited-state energies with respect to the N-electron ground-state energy, $EnN±2−E0N$⁠, with residues given by $|⟨Ψ0|ψ†ψ†|ΨnN−2⟩|2$ or $|⟨Ψ0|ψψ|ΨnN+2⟩|2$⁠. Analogous to the above, approximating K⁰ as an independent particle (or hole) pair, K⁰(1, 2) = −iG⁰(1, 2)G⁰(1, 2) and using the RPA form $K=K0+K0v¯K$ leads to the RPA pairing response function as a sum of all ladder diagrams.

In order to precisely relate the RPA excitation energies to those of an approximate EOM-CCD calculation, in Sec. II B I perform a downfolding of the ph-RPA matrix into the single particle-hole excitation subspace; in Sec. II C I show that this matrix is identical to the one obtained from EE-EOM-CCD in the single particle-hole excitation subspace when the ground-state double excitation amplitudes satisfy the ring-CCD equations and the similarity transformation of the Fock operator is neglected. Using the same approach, I show that an analogous result relates the pp-RPA excitation energies to those of an approximate DEA/DIP-EOM-CCD. Having established the algebraic equivalence of the RPA excitation energies and those from approximate EOM-CCD, in Sec. II D I analyze the time-dependent Goldstone diagrams in the RPA polarizability and pairing response function, comparing to their construction in the coupled-cluster framework, with special attention paid to the non-Tamm-Dancoff diagrams. In Sec. II E, I discuss the effect of similarity-transforming the Fock operator and relate it to a self-consistent dressing of the single-particle propagator.

The ph-RPA eigenvalue problem is given by the system of equations^1,17

where

and Ω is a diagonal matrix of positive ph-RPA excitation energies; here and throughout, I assume the typical case where the RPA stability matrix is positive definite but see Ref. 15 for a more general discussion. The antisymmetrized two-electron integrals are defined by ⟨pq||rs⟩ = ⟨pq|rs⟩−⟨pq|sr⟩, with

Formally solving the second equation, −B^*X −A^*Y = YΩ, gives

Using this expression to replace Y in the first of the ph-RPA equations leads to two equivalent forms that are eigenvalue problems for X only,

Therefore, the matrix on the left-hand side, which only has support in a single particle-hole excitation subspace (and not the subspace twice as large), has all of the positive ph-RPA excitation energies as its eigenvalues. As written, Eqs. (9) are not practical because the construction of the downfolded matrix requires knowledge of the Y component of all eigenvectors and either all eigenvalues Ω or the X component of all eigenvectors; however, the matrix in Eqs. (9) will be shown to be identical to an approximate matrix derived from EE-EOM-CCD.

Similar to the above, the pp-RPA eigenvalue problem is given by the system of equations^1,15,19 (cf. Ref. 5)

where

Ω₊ is a diagonal matrix of pp excitation energies, and Ω₋ is a diagonal matrix of hh excitation energies. Following identical algebra as above leads to the two eigenvalue equations,

These two eigenvalue problems separately give the excitation energies associated with (N ± 2)-electron systems and will be shown to be identical to an approximate solution of the DEA- and DIP-EOM-CCD equations, respectively.

In typical equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD),^30–32 the T₁ and T₂ amplitudes are obtained from the ground-state CCSD equations and the EOM eigenvalue problem is obtained by projecting the similarity-transformed normal-ordered Hamiltonian, $H¯N≡e−THeT−ECC$⁠, into a basis of singly and doubly excited determinants. By contrast, to construct the relation with the ph- or pp-RPA requires only the similarity transformation due to an approximate T₂ projected only in the minimal space of singly excited determinants, corresponding to neutral ph excitations (for ph-RPA) or doubly charged pp/hh excitations (for pp-RPA). Carrying out this procedure for the case of neutral excitations leads to

where³³ 

However, the similarity-transformed Fock operators lead to a shift of the single-particle energies in a manner which is absent in the RPA (see Sec. II E); neglecting this effect of T₂ gives

Using the definitions of CC intermediate W_ibaj and the A and B matrices leads to

where $[T]ia,jb=tijab$⁠. As shown in Refs. 1 and 17, the approximate ring-CCD equations

can be solved in closed form in terms of the eigenvectors of the ph-RPA equation (5), T = YX⁻¹. Therefore the similarity-transformed Hamiltonian, using T₂ amplitudes that satisfy the ring-CCD equations, has exactly the ph-RPA eigenvalues when truncated to the single-excitation subspace and transformation of the Fock operator is neglected. Likewise, the EOM single-excitation operator $R1(ν)=∑airia(ν)aa†ai$⁠, which gives the EOM-CC eigenstate, has amplitudes that are exactly equal to the columns of X, i.e., $ria(ν)=Xia,ν$⁠.

Although the poles of the ph-RPA polarizability are identical in the RPA and approximate EOM-CCD, the same is not true for the residues, i.e., the transition amplitudes; they are in agreement with the lowest order in the Coulomb interaction.³⁰ The RPA X and Y matrices are defined as transition density matrices,

In biorthogonal CC theory,³¹ these transition density matrices can be evaluated using the Λ-equations, and to the lowest order, the de-excitation operator amplitudes are the same as those of the T₂ amplitudes, $Λabij≈tijab$⁠. Using this approximation, it is simple to show

which agrees with the RPA interpretation after recalling the identities R = X and TX = Y.

In the same spirit, the DEA- and DIP-EOM-CCSD formalisms can be analyzed with an approximate solution of the ground-state CCSD equations and a minimal subspace for the EOM eigenvalue problem. Neglecting again the transformation of the Fock operator leads to the matrices in the 2p or 2h spaces

The CC intermediates can be straightforwardly shown to be

leading to the approximate EOM matrices

where $[T+]ij,ab=tijab$⁠. As shown in Refs. 15 and 19, the solution of the ladder-CCD equations,

can be given in terms of the solutions of the pp-RPA problem, Eq. (10), $tijab=−[Y+X+−1]ij,ab=[Y−X−−1]ab,ij*$⁠. Therefore the similarity-transformed Hamiltonian, using T₂ amplitudes that satisfy the ladder-CCD equations, has exactly the pp-RPA eigenvalues when the DEA- or DIP-EOM problems are truncated to the 2p or 2h subspaces and transformation of the Fock operator is neglected. The eigenvectors are analogously related, $rij=[X−]ij$⁠, $rab=[X+]ab$⁠.

In the language presented, the realization of the excited-state RPA problem follows from two distinct approximations (temporarily ignoring the transformation of the Fock operator): (1) an approximate solution of the ground-state CCD equations and (2) an approximate solution of the EOM eigenvalue problem. Making only one approximation, but not the other, generates two intermediate post-RPA theories for excitation energies with desirable formal features, as shown schematically in Fig. 1 for the case of EE-EOM-CCD leading to the ph-RPA. Using the exact solution of the CCD (or CCSD) amplitude equations and an approximate EOM treatment in the single excitation subspace produces the minimal fermionic theory containing RPA physics; I call this approach EOM(S)-CCD. Alternatively, using the approximate ring- or ladder-CCD ground-state solution combined with an EOM eigenvalue problem that includes doubly excited determinants (2p2h for EE-EOM or 3p1h/3h1p for DEA-/DIP-EOM) produces a flavor of “second-RPA” (or extended-RPA) theories, which include energy shifts and lifetimes due to the interaction with doubly excited configurations; I call this approach EOM-rCCD or EOM-lCCD. Making both approximations leads to an approach I call EOM(S)-rCCD or EOM(S)-lCCD, which is identical to the usual ph- or pp-RPA, as shown above (when the similarity transformation of the Fock operator is neglected). Viewed in this language, EOM-CCD is a properly fermionic theory containing both RPA and second-RPA physics.

Preliminary numerical results (not shown) indicate that the two intermediate theories (EOM(S)-CCD and EOM-rCCD/EOM-lCCD), despite their desirable formal features, do not produce quantitatively accurate excitation energies due to an unbalanced treatment of the ground state and the excited state. For example, even though EOM(S)-CCD is a properly fermionic theory that includes more diagrams, it can be viewed as a more accurate treatment of the ground state than the excited state, leading to an overestimation of the excitation energy—which is indeed observed numerically. In this sense, the RPA benefits from a favorable cancellation of errors in the ground and excited states, which is most clear in the language of approximate EOM-CCSD but relatively opaque in the traditional diagrammatic picture. Further work is needed in these directions.

The time-dependent Goldstone diagrams of the RPA polarizability are straightforward to enumerate as all ring diagrams with all possible time-orderings. In order to compare with coupled-cluster theory, a diagrammatic analysis of the coupled-cluster polarization propagator is required,³⁴ along the lines of Refs. 35 and 36 for the one-particle Green’s function. While future work will present a more detailed analysis and numerical results, the diagrams of the coupled-cluster polarization propagator can be analyzed by cutting the diagram after each vertex; each connected diagram at previous time can be classified as generated by the ground-state cluster operators, the Λ operators, or the EOM excitation operators.

Figure 2 presents some example RPA ring and ladder diagrams included through third order in perturbation theory. Vertical cuts, indicated by dashed lines, indicate that the three diagrams in Fig. 2(a) are described solely by the single-excitation EOM operator R₁. These are all examples of forward-time-ordered ring diagrams, i.e., those resulting from the Tamm-Dancoff approximation (TDA). When antisymmetrized vertices are assumed (exchange is included), then particle-hole ladder diagrams are included and the poles of the polarizability are at the excitation energies produced by configuration interaction with single excitations.

The diagrams shown in Figs. 2(b) and 2(c) are examples of non-TDA ring and ladder diagrams, due to the permuted time ordering. Graphical analysis shows that this diagram is generated through a combination of the T₂ and R₁ operators. It is straightforward to show that all non-TDA diagrams included in the RPA response functions can be deconstructed in the same manner, using disconnected products of T₂ and R₁, but never the EOM R₂ double excitation operator; this is why it was sufficient in Sec. II C to analyze the EOM eigenvalue equation in the single-excitation subspace only. Therefore, the RPA polarizability diagrams are exactly those produced by the EOM(S)-rCCD approach and the RPA ladder diagrams are those produced by the EOM(S)-lCCD approach. Incidentally, retaining the R₂ double excitation operator leads to (among many other effects) dressing of the single-particle propagator via forward time-ordered ring diagrams and ring-diagram screening of the particle-hole interaction.

All equations, as presented above, include exchange. Exchange can be trivially removed by neglecting the antisymmetrization of the two-electron integrals in the ring-CCD equations (leading to “direct” ring-CCD) and in the EOM eigenvalue problem (with a factor of 2 arising from the product of two antisymmetrized objects). This leads to a time-dependent Hartree theory of excitation energies, which is the more common variant of the RPA polarizability in the condensed-matter physics literature. Retaining exchange leads to particle-hole ladder diagrams in addition to the ring diagrams shown in Fig. 2. The particle-hole ladder diagrams are required for a description of excitonic effects in molecules or solids and are responsible for a reduction in the excitation energies compared to the time-dependent Hartree theory that only includes direct ring diagrams.

As discussed above, the similarity transformation of the Fock operator due to T₂ has been neglected, i.e., in Eqs. (14), (21), and (22). Retaining the effects of this transformation in the ground-state and excited-state calculations amounts to the use of Brueckner orbitals and orbital energies.^37,38 In diagrammatic language, the similarity transformation of the Fock operator is responsible for a dressing of the single-particle propagator G used to construct the irreducible propagator, which is still of the RPA form Π^⋆ = −iGG. In particular, the resulting G in the ph-RPA theory is very similar to that obtained from a self-energy calculated via the self-consistent GW approximation.^39,40 However, there are two major differences with respect to the GW approximation: first, only one of the two time-orderings is accounted for in the self-energy (as discussed above, the other time-ordering is non-self-consistently captured by R₂) and second, only an asymmetric subset of non-TDA ring diagrams are included in the screened Coulomb interaction. This asymmetry is identical to the behavior our group has recently described⁴¹ arising in the CCSD Green’s function approach based on ionization potential and electron affinity EOM-CCSD.^35,36 Figure 3 shows an example diagram generated by retaining this transformation.

To summarize, I have shown that the relation between the RPA and CCSD ground states can be extended to all excited states, with a particular set of additional approximations in the EOM-CCSD equations; specifically, the ph-RPA excitation energies are obtained from an approximate EE-EOM-CCSD calculation and the pp-RPA excitation energies are obtained from an approximate DEA/DIP-EOM-CCSD calculation. The exact equivalences presented here have been verified numerically, using modified implementations of the RPA and EOM-CCSD methodologies in the PySCF software package.⁴² 

In the same way that previous work¹⁷ established ground-state CCD as the natural generalization of the RPA with correct fermionic behavior, the present work proposes EOM(S)-CCD as the simplest fermionic theory of excited states that contains RPA physics. This generalization comes with a cost: for a single low-lying excited state, an RPA calculation scales as N⁴, whereas an EOM(S)-CCD calculation (as well as an EOM(S)-rCCD or EOM(S)-lCCD calculation) scales as N⁶ due to the ground-state CCD step; the ensuing EOM(S) step only has N⁴ scaling. This overall N⁶ scaling is no worse than that of EOM-CCSD, which is clearly preferred for a few low-lying excited states. However, the cost to obtain all excited states is N⁶ for both RPA and EOM(S)-CCD, to be compared to N⁸ for EOM-CCSD (for all excited states with dominant single-excitation character), which may be important for certain spectral quantities.

In addition to providing a properly fermionic theory, the present manuscript establishes the RPA polarizability diagrams as a strict subset of those from EOM-CCSD. In this sense, the CC hierarchy is a natural post-RPA route, distinct from time-dependent density functional theory and, importantly, systematically improvable. It is hoped that this connection will lead to fruitful developments in the simulation of excited states, especially in the condensed phase where RPA physics is essential. For example, various CC-derived polarizabilities can be used for a more accurate treatment of screening in the GW approximation, leading to a well-defined class of vertex corrections. Similarly, a comparison of EOM-CCSD excited states to those predicted by the GW + Bethe-Salpeter equation approach will provide further insight and sow deeper connections between the condensed-matter and quantum chemistry communities. Work along both of these lines is currently in progress.

I thank Alan Lewis and Bryan Lau for useful conversations and comments on this manuscript, and an anonymous reviewer for bringing results of Ref. 30 to my attention. This work was supported in part by startup funds from the University of Chicago and by the Air Force Office of Scientific Research under Award No. FA9550-18-1-0058.