Particle-particle and quasiparticle random phase approximations: Connections to coupled cluster theory

Gustavo E. Scuseria Thomas M. Henderson Department of Chemistry and Department of Physics and Astronomy, Rice University, Houston, TX 77005-1892 Ireneusz W. Bulik Department of Chemistry, Rice University, Houston, TX 77005-1892

(Updated August 22, 2025)

Abstract

We establish a formal connection between the particle-particle (pp) random phase approximation (RPA) and the ladder channel of the coupled cluster doubles (CCD) equations. The relationship between RPA and CCD is best understood within a Bogoliubov quasiparticle (qp) RPA formalism. This work is a follow-up to our previous formal proof on the connection between particle-hole (ph) RPA and ring-CCD. Whereas RPA is a quasibosonic approximation, CC theory is a correct bosonization in the sense that the wavefunction and Hilbert space are exactly fermionic. Coupled cluster theory achieves this goal by interacting the ph (ring) and pp (ladder) diagrams via a third channel that we here call “crossed-ring” whose presence allows for full fermionic antisymmetry. Additionally, coupled cluster incorporates what we call “mosaic” terms which can be absorbed into defining a new effective one-body Hamiltonian. The inclusion of these mosaic terms seems to be quite important. The pp-RPA and qp-RPA equations are textbook material in nuclear structure physics but are largely unknown in quantum chemistry, where particle number fluctuations and Bogoliubov determinants are rarely used. We believe that the ideas and connections discussed in this paper may help design improved ways of incorporating RPA correlation into density functionals based on a CC perspective.

I Introduction

In a previous paper,Scuseria et al. (2008) we established a connection between the particle-hole random phase approximation (RPA) and the ring channel of the coupled cluster doubles (CCD) equations. Here, we extend this analysis to the ladder channel of CCD and demonstrate a rigorous connection with the particle-particle (pp) RPA equations. RPA is a quasibosonic approximation in the sense that fermion products are treated as bosons when in reality these operators satisfy Lie algebra commutation rules that are neither fermionic nor bosonic. The RPA quasibosonic approximation contaminates the fermionic Hilbert space with bosonic states, thus leading to systematic overestimation of ground state correlation energies. We here also show that from a Bogoliubov quasiparticle (qp) RPA perspective, the particle-hole and particle-particle RPA channels get added together but do not interact. On the other hand, CCD can be interpreted as a “correct bosonization” because a third channel, here referred to as “crossed ring,” allows the other two channels to interact and closes the equations in a manner that exactly preserves the fermionic nature of the wavefunction and the Hilbert space of the problem. The particle-particle RPA and quasiparticle RPA equations are textbook material in nuclear structure physics but are largely unknown in quantum chemistry, where particle number fluctuations and HFB determinants are rarely used. Coupled cluster theory has undoubtedly been very successful in quantum chemistry and given the current interest of using RPA to improve DFT,Ren et al. (2012); Eshuis et al. (2012) we believe that the ideas and connections discussed in this paper may help design improved ways of incorporating RPA effects into functionals based on a coupled cluster perspective.

II Different Flavors of RPA

II.1 Particle-Hole RPA

The traditional RPA formulation can be interpreted as an attempt to treat particle-hole fermionic excitations as bosons. RPA is a theory that can be used both for excited states and ground state correlation. The standard derivation of the RPA equations usually follows a so-called equation of motion approach.Ring and Schuck (1980) Here, we pursue a somewhat different perspective. Products of particle-hole excitation fermionic operators are assumed to be bosons

a_{a}^{\dagger}\,a_{i}\rightarrow b_{\beta}

(1)

that satisfy commutation rules

	$\displaystyle\left[b_{\beta},b_{\beta^{\prime}}^{\dagger}\right]$	$\displaystyle=\delta_{\beta\beta^{\prime}}$		(2)
	$\displaystyle\left[b_{\beta},b_{\beta^{\prime}}\right]$	$\displaystyle=\left[b_{\beta}^{\dagger},b_{\beta^{\prime}}^{\dagger}\right]=0.$		(3)

Note that these commutation rules are satisfied on average when using a single-determinant reference.

This quasiboson approximation is the central approximation of RPA. In reality, the $b_{\beta}$ fermionic product operators satisfy unitary group $U(M)$ commutation rules that are those of a Lie algebra

\left[a_{a}^{\dagger}\,a_{i},a_{j}^{\dagger}\,a_{b}\right]=\delta_{ij}\,a_{a}^{\dagger}\,a_{b}-\delta_{ab}\,a_{j}^{\dagger}\,a_{i}.

(4)

We follow the traditional notation where spin orbitals $i,j,k,l$ are occupied (holes) and $a,b,c,d$ are unoccupied (particles) in a reference determinant. Indices $p,q,r,s,\ldots$ refer to unspecified spin orbitals. All repeated indices are summed.

In RPA, the quartic fermionic Hamiltonian

H=h_{pq}\,a_{p}^{\dagger}\,a_{q}+\frac{1}{4}\langle pq\|rs\rangle\,a_{p}^{\dagger}\,a_{q}^{\dagger}\,a_{s}a_{r}

(5)

is, after particle-hole transformation, interpreted as a quadratic bosonic Hamiltonian with number-violating terms:

H=E_{\mathrm{HF}}+A_{\beta\beta^{\prime}}\,b_{\beta}^{\dagger}\,b_{\beta^{\prime}}+\frac{1}{2}\,\left(B_{\beta\beta^{\prime}}b_{\beta}^{\dagger}b_{\beta^{\prime}}^{\dagger}+B_{\beta\beta^{\prime}}^{\star}\,b_{\beta}\,b_{\beta^{\prime}}\right).

(6)

Here, the Hermitian $\bm{A}=\bm{A}^{\dagger}$ and symmetric $\bm{B}=\bm{B}^{\mathsf{T}}$ matrices are those of standard RPA:


$\displaystyle A_{\beta\beta^{\prime}}$	$\displaystyle\rightarrow A_{ai,bj}=f_{ab}\,\delta_{ij}-f_{ji}\,\delta_{ab}+\langle aj\\|ib\rangle,$	(7a)
$\displaystyle B_{\beta\beta^{\prime}}$	$\displaystyle\rightarrow B_{ai,bj}=\langle ij\\|ab\rangle,$	(7b)

where $\bm{f}$ is the Fock matrix. Using the bosonic commutations relations, one can recognize the Hamiltonian as

H=E_{\mathrm{HF}}+\frac{1}{2}\begin{pmatrix}\bm{b}^{\dagger}&\bm{b}\end{pmatrix}\begin{pmatrix}\bm{A}&\bm{B}\\ \bm{B}^{\star}&\bm{A}^{\star}\end{pmatrix}\begin{pmatrix}\bm{b}\\ \bm{b}^{\dagger}\end{pmatrix}-\frac{1}{2}\,\mathrm{Tr}(\bm{A}).

(8)

This is a quadratic bosonic Hamiltonian that can be solved exactly (diagonalized) via a Bogoliubov transformation for bosons.Blaizot and Ripka (1985) The solution is obtained via a non-Hermitian diagonalization problem

\begin{pmatrix}\bm{A}&\bm{B}\\ -\bm{B}^{\star}&-\bm{A}^{\star}\end{pmatrix}\begin{pmatrix}\bm{X}&\bm{Y}^{\star}\\ \bm{Y}&\bm{X}^{\star}\end{pmatrix}=\begin{pmatrix}\bm{X}&\bm{Y}^{\star}\\ \bm{Y}&\bm{X}^{\star}\end{pmatrix}\begin{pmatrix}\bm{\omega}&\bm{0}\\ \bm{0}&-\bm{\omega}^{\star}\end{pmatrix}.

(9)

For the solution to be physically meaningful, we must have

\bm{M}=\begin{pmatrix}\bm{A}&\bm{B}\\ \bm{B}^{\star}&\bm{A}^{\star}\end{pmatrix}>=0

(10)

so that all of the RPA eigenvalues $\bm{\omega}$ and the correlation energy

E_{c}=\frac{1}{2}\sum_{\omega>0}\omega-\frac{1}{2}\,\mathrm{Tr}(\bm{A})

(11)

are real.

General quadratic fermionic Hamiltonians (i.e. those containing number-conserving and number-violating terms as the bosonic Hamiltonian above) are diagonalized via Bogoliubov canonical transformations defining a Hartree-Fock-Bogoliubov (HFB) determinant, which is a product state of quasiparticles. The simplest possible bosonic wavefunction is defined by the diagonalization of $\bm{A}$ , which is the so-called Hartree-Bose (HB) problem. Hartree-Bose yields a condensate wavefunction where every boson occupies the same orbital so that the density matrix is simply $\rho_{pq}=z_{p}z_{q}^{\star}=(\bm{z}\,\bm{z}^{\dagger})_{pq}$ . Including $\bm{B}$ but retaining the condition that every boson occupies the same orbital specified by $\bm{z}$ leads to a coherent state in Fock space, and is the result of a simple shift canonical transformation. These two wavefunctions are essentially equivalent in the context of RPA. The full quadratic bosonic Hamiltonian diagonalization involves a Bogoliubov bosonic transformation (shift plus rotation) that introduces fluctuations over the condensate; these fluctuations can be thought of as correlations. When we approximate fermionic product operators as bosons, we take advantage of this bosonic structure to obtain ground state correlations via a quadratic Hamiltonian diagonalization which yields a correlated wavefunction essentially because the bosons are composite fermions. One disadvantage of this approach is that the dimension of the bosonic Hamiltonian is proportional to the square of the number of fermions (strictly, the number of particle-hole excitations).

From this perspective, the RPA for fermionic excitations and ground-state correlations is a bosonic treatment that includes two fundamental choices:

•

The effective bosonic excitation operators in particle-hole RPA are written as

$b_{\nu}^{\dagger}\rightarrow Q_{\nu}^{\dagger}=X_{ai}^{\nu}\,a_{a}^{\dagger}\,a_{i}-Y_{ai}^{\nu}\,a_{i}^{\dagger}\,a_{a}.$ (12)
•

The fermionic reference state in traditional particle-hole RPA is simply the Hartree-Fock determinant, $|\mathrm{HF}\rangle$ .

Other choices for excitation operators and the reference state are possible.

The fact that products of fermion operators are not bosons manifests in RPA in myriad ways. To name a few and without going into details, the lack of a killer condition, the difficulty in obtaining a self-consistent ground-state RPA approximation, the indefinition of entire blocks of the two-particle reduced density matrix (RDM), the appearance of non-representable RDMs, violations to Pauli’s principle, and the presence of unphysical (bosonic) states in the fermionic Hilbert space of the problem are all symptoms of the same condition: the quasiboson approximation.

II.2 Particle-Particle RPA

In this approach, one simply considers non-number–conserving excitations and de-excitationsRing and Schuck (1980); Blaizot and Ripka (1985)

Q_{\nu}^{\dagger}=\frac{1}{2}X^{\nu}_{ab}\,a^{\dagger}_{a}\,a^{\dagger}_{b}-\frac{1}{2}Y^{\nu}_{ij}a^{\dagger}_{j}\,a^{\dagger}_{i}.

(13)

As shown below, this leads to a different particle-particle RPA problem, one that is seldom discussed in quantum chemistry. From a coupled cluster perspective, the contractions in particle-particle RPA are ladders rather than the rings in particle-hole RPA. A formal analytic proof is also presented showing the equivalence between particle-particle RPA and ladder-CCD, which is a natural follow-up to our previous proof of the equivalence between particle-hole RPA and ring-CCD.Scuseria et al. (2008) In particle-particle RPA too, the commutation rules between the excitation operators are approximated as being bosonic whereas in reality they are those of an $SO(2M)$ Lie algebra, one that includes the number conserving ph subalgebra of $U(M)$ .

While one can cast particle-particle RPA into a symplectic eigenvalue problem, and from it extract a ground-state correlation energy, we do not show the detailed expressions in this subsection. Rather, we prefer to wait until we have discussed the more general quasiparticle RPA, which subsumes both particle-hole and particle-particle RPA into a single diagonalization problem.

II.3 Quasiparticle RPA

The particle-hole and particle-particle forms of the excitation operators are special cases of a more general quasiparticle excitation operator that one can write asRing and Schuck (1980)

Q_{\nu}^{\dagger}=\frac{1}{4}\mathcal{X}_{pq}^{\nu}\,\alpha_{p}^{\dagger}\,\alpha_{q}^{\dagger}-\frac{1}{4}\mathcal{Y}_{pq}^{\nu}\,\alpha_{q}\,\alpha_{p}

(14)

where the $\alpha$ are Bogoliubov canonically transformed fermionic quasiparticle operators

	$\displaystyle\begin{pmatrix}\bm{\alpha}\\ \bm{\alpha}^{\dagger}\end{pmatrix}$	$\displaystyle=\begin{pmatrix}\bm{U}^{\dagger}&\bm{V}^{\dagger}\\ \bm{V}^{\mathsf{T}}&\bm{U}^{\mathsf{T}}\end{pmatrix}\begin{pmatrix}\bm{a}\\ \bm{a}^{\dagger}\end{pmatrix}=\bm{W}^{\dagger}\begin{pmatrix}\bm{a}\\ \bm{a}^{\dagger}\end{pmatrix}$		(15)
	$\displaystyle\bm{W}^{\dagger}\,\bm{W}$	$\displaystyle=\bm{W}\,\bm{W}^{\dagger}=\bm{1}$		(16)

preserving anticommutation rules. In this quasiparticle basis, the Hamiltonian takes the form

\begin{split}H&=H^{0}+H^{11}_{pq}\,\alpha_{p}^{\dagger}\,\alpha_{q}+\frac{1}{2}\left(H^{20}_{pq}\,\alpha^{\dagger}_{p}\,\alpha^{\dagger}_{q}+H_{qp}^{20\star}\,\alpha_{q}\,\alpha_{p}\right)+\frac{1}{4}\,H^{22}_{pqrs}\,\alpha_{p}^{\dagger}\,\alpha_{q}^{\dagger}\,\alpha_{r}\,\alpha_{s}\\ &+H^{40}_{pqrs}\,\alpha_{p}^{\dagger}\,\alpha_{q}^{\dagger}\,\alpha_{r}^{\dagger}\,\alpha_{s}^{\dagger}+H^{40\star}_{sqrp}\,\alpha_{s}\,\alpha_{r}\,\alpha_{q}\,\alpha_{p}+H^{31}_{pqrs}\,\alpha_{p}^{\dagger}\,\alpha_{q}^{\dagger}\,\alpha_{r}^{\dagger}\,\alpha_{s}+H^{31\star}_{sqrp}\,\alpha_{s}^{\dagger}\,\alpha_{r}\,\alpha_{q}\,\alpha_{p}.\end{split}

(17)

Detailed expressions for the matrix elements can be found in Appendix E of Ref. Ring and Schuck, 1980. From the quasiparticle mean-field approximation we obtain quasiparticle energies $E_{k}$ .

Assuming a quasiboson approximation for the fermionic quasiparticle products of Eqn. 14 and taking $|\mathrm{HFB}\rangle$ as a reference, one finds an RPA problem of the form

\begin{pmatrix}\bm{\mathcal{A}}&\bm{\mathcal{B}}\\ -\bm{\mathcal{B}}^{\star}&-\bm{\mathcal{A}}^{\star}\end{pmatrix}\begin{pmatrix}\bm{\mathcal{X}}&\bm{\mathcal{Y}}^{\star}\\ \bm{\mathcal{Y}}&\bm{\mathcal{X}}^{\star}\end{pmatrix}=\begin{pmatrix}\bm{\mathcal{X}}&\bm{\mathcal{Y}}^{\star}\\ \bm{\mathcal{Y}}&\bm{\mathcal{X}}^{\star}\end{pmatrix}\begin{pmatrix}\bm{\omega}&\bm{0}\\ \bm{0}&-\bm{\omega}\end{pmatrix}

(18)

where the indices now run over all pairs. The RPA matrices are


$\displaystyle\mathcal{A}_{pq,rs}$	$\displaystyle=\langle\Big{[}\alpha_{q}\alpha_{p}\Big{[}H,\alpha{}^{\dagger}_{r}\alpha{}^{\dagger}_{s}\Big{]}\Big{]}\rangle=\left(E_{p}+E_{q}\right)\delta_{pr}\delta_{qs}+H_{pqrs}^{22}$	(19a)
$\displaystyle\mathcal{B}_{pq,rs}$	$\displaystyle=\langle\Big{[}\alpha_{q}\alpha_{p}\Big{[}H,\alpha{}_{s}\alpha{}_{r}\Big{]}\Big{]}\rangle=4!H_{pqrs}^{40}.$	(19b)

Again, $\bm{\mathcal{A}}$ is Hermitian and $\bm{\mathcal{B}}$ is symmetric. Note that the quasiparticle RPA matrix is $\bigl{(}\begin{smallmatrix}\bm{\mathcal{A}}&\bm{\mathcal{B}}\\ -\bm{\mathcal{B}}^{\star}&-\bm{\mathcal{A}}^{\star}\end{smallmatrix}\bigr{)}=\bm{\eta}\,\bm{\mathcal{M}}$ , where $\bm{\eta}=\bigl{(}\begin{smallmatrix}\bm{1}&\bm{0}\\ \bm{0}&\bm{-1}\end{smallmatrix}\bigr{)}$ is a symplectic metric and $\bm{\mathcal{M}}$ is the HFB orbital Hessian.

The quasiparticle RPA matrix has the standard symplectic form, and therefore has properties similar to particle-hole RPA. In particular, qp-RPA leads to a Ricatti equationScuseria et al. (2008)

\bm{\mathcal{B}}^{\star}+\bm{\mathcal{A}}^{\star}\,\bm{\mathcal{T}}+\bm{\mathcal{T}}\,\bm{\mathcal{A}}+\bm{\mathcal{T}}\,\bm{\mathcal{B}}\,\bm{\mathcal{T}}=\bm{0}

(20)

with

\bm{\mathcal{T}}=\bm{\mathcal{Y}}\,\bm{\mathcal{X}}^{-1}

(21)

and a corresponding correlation energy

E_{c}=\frac{1}{2}\sum_{\omega>0}\omega-\frac{1}{2}\,\mathrm{Tr}(\bm{\mathcal{A}}).

(22)

From the form of $Q_{\nu}^{\dagger}$ used in qp-RPA, it is evident that qp-RPA does not treat the $31$ and $13$ blocks of $H$ ; in other words $H^{31},$ which is responsible for connecting the particle-particle and particle-hole channels, does not appear in qp-RPA. Below we discuss how single-reference CCD theory makes these channels interact, essentially by demanding that these channels (together with a crossed-ring channel we will introduce later) lead to the same fermionic wave function amplitudes. We christen this process as a “correct bosonization.”

Let us now consider a special case of quasiparticle RPA, namely, the limit when the HFB reference determinant reduces to Hartree-Fock. For systems with purely repulsive electron-electron interactions, this is the variationally optimal result.Bach et al. (1994); Scuseria and Tsuchimochi (2009) In this case, the RPA matrix greatly simplifies, and takes the blocked form

\bm{\eta}\,\bm{\mathcal{M}}=\begin{pmatrix}\bm{\mathcal{A}}_{oo,oo}&\bm{0}&\bm{0}&\bm{0}&\bm{0}&\bm{\mathcal{B}}_{oo,vv}\\ \bm{0}&\bm{\mathcal{A}}_{ov,ov}&\bm{0}&\bm{0}&\bm{\mathcal{B}}_{ov,ov}&\bm{0}\\ \bm{0}&\bm{0}&\bm{\mathcal{A}}_{vv,vv}&\bm{\mathcal{B}}_{vv,oo}&\bm{0}&\bm{0}\\ \bm{0}&\bm{0}&-\bm{\mathcal{B}}_{oo,vv}^{\star}&-\bm{\mathcal{A}}_{oo,oo}^{\star}&\bm{0}&\bm{0}\\ \bm{0}&-\bm{\mathcal{B}}_{ov,ov}^{\star}&\bm{0}&\bm{0}&-\bm{\mathcal{A}}_{ov,ov}^{\star}&\bm{0}\\ -\bm{\mathcal{B}}_{vv,oo}^{\star}&\bm{0}&\bm{0}&\bm{0}&\bm{0}&-\bm{\mathcal{A}}_{vv,vv}^{\star}\end{pmatrix}

(23)

where subscripts “oo”, “ov”, and “vv” refer to occupied-occupied, occupied-virtual, and virtual-virtual, respectively. The indices here run over unique orbital pairs, so “oo” has indices $i<j$ , “ov” has indices $ia$ , and “vv” has indices $a<b$ . We should note that when we have a genuine HFB reference which does not conserve particle number, “occupied” and “virtual” lose their meaning and the indices of the quasiparticle RPA matrix correspond to quasiparticle creation and annihilation operators, as seen from Eqn. 14.

To make contact with our later discussion, we note that the matrices $\bm{\mathcal{A}}_{ov,ov}$ and $\bm{\mathcal{B}}_{ov,ov}$ are the matrices $\bm{A}$ and $\bm{B}$ of particle-hole RPA, and we will define


$\displaystyle\bm{\mathcal{A}}_{oo,oo}$	$\displaystyle=\bm{D},$	(24a)
$\displaystyle\bm{\mathcal{A}}_{vv,vv}$	$\displaystyle=\bm{C},$	(24b)
$\displaystyle\bm{\mathcal{B}}_{vv,oo}$	$\displaystyle=-\bar{\bm{B}}.$	(24c)

Using the fact that $\bm{\mathcal{B}}_{oo,vv}=(\bm{\mathcal{B}}_{vv,oo})^{\mathsf{T}}$ , the quasiparticle RPA matrix $\bm{\mathcal{M}}$ expressed in this notation is

\bm{\eta}\,\bm{\mathcal{M}}=\begin{pmatrix}\bm{D}&\bm{0}&\bm{0}&\bm{0}&\bm{0}&-\bar{\bm{B}}^{\mathsf{T}}\\ \bm{0}&\bm{A}&\bm{0}&\bm{0}&\bm{B}&\bm{0}\\ \bm{0}&\bm{0}&\bm{C}&-\bar{\bm{B}}&\bm{0}&\bm{0}\\ \bm{0}&\bm{0}&\bar{\bm{B}}^{\dagger}&-\bm{D}^{\star}&\bm{0}&\bm{0}\\ \bm{0}&-\bm{B}^{\star}&\bm{0}&\bm{0}&-\bm{A}^{\star}&\bm{0}\\ \bar{\bm{B}}^{\star}&\bm{0}&\bm{0}&\bm{0}&\bm{0}&-\bm{C}^{\star}\end{pmatrix}

(25)

The matrix elements of $\bm{C}$ , $\bm{D}$ , and $\bar{\bm{B}}$ are


$\displaystyle D_{ij,kl}$	$\displaystyle=-(\epsilon_{i}+\epsilon_{j})\,\delta_{ik}\,\delta_{jl}+\langle kl\\|ij\rangle,$	(26a)
$\displaystyle C_{ab,cd}$	$\displaystyle=(\epsilon_{a}+\epsilon_{b})\,\delta_{ac}\,\delta_{bd}+\langle ab\\|cd\rangle,$	(26b)
$\displaystyle\bar{B}_{ab,ij}$	$\displaystyle=\langle ab\\|ij\rangle.$	(26c)

It is clear that one can decompose the quasiparticle RPA into subproblems. One subproblem gives us the usual particle-hole RPA. Taking the central block of $\bm{\mathcal{M}}$ gives us particle-particle RPA, wherein one solvesRing and Schuck (1980)

\begin{pmatrix}\bm{C}&-\bar{\bm{B}}\\ \bar{\bm{B}}^{\dagger}&-\bm{D}^{\star}\end{pmatrix}\begin{pmatrix}\bm{X}_{1}&\bm{Y}_{2}\\ \bm{Y}_{1}&\bm{X}_{2}\end{pmatrix}=\begin{pmatrix}\bm{X}_{1}&\bm{Y}_{2}\\ \bm{Y}_{1}&\bm{X}_{2}\end{pmatrix}\begin{pmatrix}\bm{\Omega_{1}}&\bm{0}\\ \bm{0}&\bm{\Omega}_{2}\end{pmatrix}.

(27)

The frequencies $\bm{\Omega}_{1}$ are positive and $\bm{\Omega}_{2}$ are negative. The remaining portion of $\bm{\mathcal{M}}$ gives hole-hole RPA, the symplectic counterpart to particle-particle RPA. Note that particle-particle and hole-hole RPA are not individually symplectic eigenvalue problems, so one cannot straightforwardly extract a plasmonic correlation energy from just one or the other. Rather, they should be grouped together as

\begin{pmatrix}\bm{D}&\bm{0}&\bm{0}&-\bar{\bm{B}}^{\mathsf{T}}\\ \bm{0}&\bm{C}&-\bar{\bm{B}}&\bm{0}\\ \bm{0}&\bar{\bm{B}}^{\dagger}&-\bm{D}^{\star}&\bm{0}\\ \bar{\bm{B}}^{\star}&\bm{0}&\bm{0}&-\bm{C}^{\star}\end{pmatrix}\begin{pmatrix}\bm{X}_{2}^{\star}&\bm{0}&\bm{0}&\bm{Y}_{1}^{\star}\\ \bm{0}&\bm{X}_{1}&\bm{Y}_{2}&\bm{0}\\ \bm{0}&\bm{Y}_{1}&\bm{X}_{2}&\bm{0}\\ \bm{Y}_{2}^{\star}&\bm{0}&\bm{0}&\bm{X}_{1}^{\star}\end{pmatrix}=\begin{pmatrix}\bm{X}_{2}^{\star}&\bm{0}&\bm{0}&\bm{Y}_{1}^{\star}\\ \bm{0}&\bm{X}_{1}&\bm{Y}_{2}&\bm{0}\\ \bm{0}&\bm{Y}_{1}&\bm{X}_{2}&\bm{0}\\ \bm{Y}_{2}^{\star}&\bm{0}&\bm{0}&\bm{X}_{1}^{\star}\end{pmatrix}\begin{pmatrix}-\bm{\Omega}_{2}&\bm{0}&\bm{0}&\bm{0}\\ \bm{0}&\bm{\Omega}_{1}&\bm{0}&\bm{0}\\ \bm{0}&\bm{0}&\bm{\Omega}_{2}&\bm{0}\\ \bm{0}&\bm{0}&\bm{0}&-\bm{\Omega}_{1}\end{pmatrix}.

(28)

From here, the plasmonic correlation energy is just

E_{c}=\frac{1}{2}\mathrm{Tr}\left(\bm{\Omega}_{1}-\bm{\Omega}_{2}-\bm{C}-\bm{D}\right).

(29)

Using the particle-particle RPA equations, one sees that

\mathrm{Tr}(\bm{C})-\mathrm{Tr}(\bm{D}^{\star})=\mathrm{Tr}(\bm{\Omega}_{1})+\mathrm{Tr}(\bm{\Omega}_{2})

(30)

so that

\mathrm{Tr}(\bm{\Omega}_{1}-\bm{C})=-\mathrm{Tr}(\bm{\Omega}_{2}+\bm{D}^{\star}).

(31)

Since $\bm{D}$ is Hermitian so that $\mathrm{Tr}(\bm{D})=\mathrm{Tr}(\bm{D}^{\star})$ , it follows that the plasmonic correlation energy associated with particle-particle/hole-hole RPA can be equivalently expressed as

E_{c}=\mathrm{Tr}(\bm{\Omega}_{1}-\bm{C})=-\mathrm{Tr}(\bm{\Omega}_{2}+\bm{D}^{\star}).

(32)

Note finally that because quasiparticle RPA in this limit can be factored into two symplectic subproblems, the correlation energy associated with quasiparticle RPA is simply additive:

E_{c,qp-RPA}=E_{c,ph-RPA}+E_{c,pp-RPA}.

(33)

II.4 Stability of RPA Problems

The three types of RPA we have discussed above all have the same formal symplectic structure when the particle-particle RPA is understood as particle-particle/hole-hole RPA. All diagonalize a matrix

\bm{\eta}\,\mathbf{M}=\begin{pmatrix}\mathbf{A}&\mathbf{B}\\ -\mathbf{B}^{\star}&-\mathbf{A}^{\star}\end{pmatrix}

(34)

where we recall that

\bm{\eta}=\begin{pmatrix}\bm{1}&\bm{0}\\ \bm{0}&-\bm{1}\end{pmatrix}.

(35)

The matrix $\mathbf{M}$ is an orbital Hessian, or, to put it another way, it is a stability matrix. For particle-hole RPA, diagonalizing $\mathbf{M}$ tests for instabilities of the reference $|\mathrm{HF}\rangle$ solution toward another Hartree-Fock state. For particle-particle/hole-hole RPA, diagonalizing $\mathbf{M}$ tests for instabilities of the reference $|\mathrm{HF}\rangle$ towards an HFB state, while for quasiparticle RPA, diagonalizing $\mathbf{M}$ tests for instabilities of the reference $|\mathrm{HFB}\rangle$ toward another HFB state. Stability in this context means that $\mathbf{M}>=0$ . Earlier, we noted that particle-hole RPA gives physically meaningful results when this condition is satisfied; the same holds for particle-particle and for quasiparticle RPA. Quite generally, when $\mathbf{M}$ is positive definite, the corresponding RPA problem has real eigenvalues and $2n$ linearly independent eigenvectors where $\mathbf{M}$ is $2n\times 2n$ .Colpa (1986)

From a minor modification of the proof in Appendix 5 of Ref. Scuseria et al., 2008, we see that a Riccati equation of the form $\mathbf{B}^{\star}+\mathbf{A}^{\star}\,\mathbf{T}+\mathbf{T}\,\mathbf{A}+\mathbf{T}\,\mathbf{B}\,\mathbf{T}=\mathbf{0}$ implies an eigenvalue-like problem

\begin{pmatrix}\mathbf{A}&\mathbf{B}\\ -\mathbf{B}^{\star}&-\mathbf{A}^{\star}\end{pmatrix}\begin{pmatrix}\mathbf{x}\\ \mathbf{y}\end{pmatrix}=\begin{pmatrix}\mathbf{x}\\ \mathbf{y}\end{pmatrix}\bm{\Delta}.

(36)

Here, we have written


$\displaystyle\mathbf{A}+\mathbf{B}\,\mathbf{T}$	$\displaystyle=\mathbf{x}\,\bm{\Delta}\,\mathbf{x}^{\dagger}$	(37a)
$\displaystyle\mathbf{y}$	$\displaystyle=\mathbf{T}\,\mathbf{x}$	(37b)

where $\mathbf{x}$ is unitary and $\bm{\Delta}$ is upper triangular. When the matrix $\bm{\eta}\mathbf{M}$ is diagonalizable, the Riccati equation implies the eigenvalue problem itself. Whether $\bm{\eta}\mathbf{M}$ is diagonalizable or not, we have

\mathrm{Tr}(\mathbf{B}\,\mathbf{T})=\mathrm{Tr}(\bm{\Delta}-\mathbf{A})

(38)

where we recall that the diagonal elements of $\bm{\Delta}$ are its eigenvalues.

In other words, when the RPA matrix $\bm{\eta}\mathbf{M}$ is diagonalizable, the eigenvectors of the RPA problem are intimately connected to a corresponding Riccati equation. When the reference is stable, the RPA matrix is diagonalizable. We cannot say much about the stability of a general Hartree-Fock or HFB determinant, but we note that for a repulsive two-body interaction, HF is never unstable toward HFB.Bach et al. (1994); Scuseria and Tsuchimochi (2009) Thus, particle-particle RPA applies to the standard Coulombic Hamiltonian should not suffer from complex correlation energies. A more detailed proof is forthcoming. In the following sections, we will use a simpler proof which assumes invertibility of $\mathbf{X}$ to derive the Riccati equation from the RPA eigenvalue problem, but the foregoing shows that we do not need $\mathbf{X}^{-1}$ to exist for a relation between RPA and a Riccati equation to be found.

III Different Flavors of CCD

Having discussed RPA at some length, we now turn to coupled cluster doubles theory.Bartlett and Shavitt (2009) As we shall see, CCD contains various pieces of the RPA problem and unifies them all in such a way as to preserve the fermionic character of the wavefunction. It will be convenient for our purposes to work with the Brueckner version of CCD theory (referred to as BD) that eliminates single excitations. This section closely follows the notation of Ref. Scuseria, 1995. The basic ingredients of the BD model are one ( $h$ ) and two-electron ( $v$ ) integrals, and cluster amplitudes ( $t$ ) in the spin-orbital basis


$\displaystyle f_{p}^{q}$	$\displaystyle=h_{p}^{q}+v_{pq}^{qk},$	(39a)
$\displaystyle v^{ij}_{ab}$	$\displaystyle=\langle ij\\|ab\rangle=\langle ij\|ab\rangle-\langle ij\|ba\rangle,$	(39b)
$\displaystyle v^{ab}_{ij}$	$\displaystyle=\langle ab\\|ij\rangle=\langle ij\\|ab\rangle^{\star}=\left(v^{ij}_{ab}\right)^{\star},$	(39c)
$\displaystyle t_{ij}^{ab}$	$\displaystyle=\langle ab\|t_{2}\|ij\rangle=-t_{ji}^{ab}=-t_{ij}^{ba}=t_{ji}^{ba},$	(39d)
$\displaystyle E$	$\displaystyle=E_{0}+E_{c},$	(39e)
$\displaystyle E_{0}$	$\displaystyle=\frac{1}{2}\,\left(h^{i}_{i}+f^{i}_{i}\right)=h^{i}_{i}+\frac{1}{2}v^{ik}_{ik},$	(39f)
$\displaystyle E_{c}$	$\displaystyle=\frac{1}{4}\,v^{ij}_{ab}\,t_{ij}^{ab}.$	(39g)

Repeated indices are always summed (even in $h_{i}^{i}$ ). Upper and lower indices can be identified as bra and ket, respectively.

We now define a Brueckner effective one-body Hamiltonian through the energy expression

E=\frac{1}{2}\left(h_{i}^{i}+F_{i}^{i}\right)

(40)

which defines the occupied-occupied block of $\bm{F}$

F_{i}^{k}=f_{i}^{k}+\frac{1}{2}\,v_{ab}^{kj}\,t_{ij}^{ab}

(41)

Using particle-hole symmetry, we impose

F_{c}^{a}=f_{c}^{a}-\frac{1}{2}\,v_{cb}^{ij}\,t_{ij}^{ab}

(42)

which is needed to complete the CCD equations in the desired form (vide infra). The occupied-virtual block of $F$ , which we force to be zero, is chosen as the $T_{1}$ equation within the BD approximation (i.e., enforcing zero $T_{1}$ amplitudes)

F_{i}^{a}=f_{i}^{a}+f_{b}^{j}\,t_{ij}^{ab}+\frac{1}{2}\,v_{bc}^{aj}\,t_{ij}^{bc}-\frac{1}{2}\,v_{ib}^{jk}\,t_{jk}^{ab}=0.

(43)

The three-index contraction terms (e.g., $v_{ab}^{kj}\,t_{ij}^{ab}$ ) appearing in the $F$ equations above are here referred to as mosaic (they have ring and ladder contractions). The BD equations become simply

\begin{split}0&=v_{ij}^{ab}+\frac{1}{2}\,t_{kl}^{ab}\,v_{ij}^{kl}+\frac{1}{2}\,v_{cd}^{ab}\,t_{ij}^{cd}+\frac{1}{4}\,t_{kl}^{ab}\,v_{cd}^{kl}\,t_{ij}^{cd}\\ &-P_{ij}\left(t_{kj}^{ab}\,F_{i}^{k}\right)+P_{ab}\left(F_{c}^{a}\,y_{ij}^{cb}\right)+P_{ij}P_{ab}\left[\left(v_{ic}^{ak}+\frac{1}{2}\,t_{il}^{ad}\,v_{cd}^{kl}\right)t_{kj}^{cb}\right]\end{split}

(44)

where $P$ is an index permutation operator (e.g., $P_{ij}=1-i\leftrightarrow j$ ).

Let us now analyze the BD amplitude equation:

•

The first term $v_{ij}^{ab}$ is called the driver.
•

The next three terms contain pp or hh (ladder) contractions only.
•

The third term is quadratic in the amplitudes; this is the highest degree of the equations.

•

The next two terms are most readily understood in the Brueckner canonical basis where $F$ is diagonal with eigenvalues $\zeta$ :

P_{ab}\left(F^{a}_{c}\,t_{ij}^{cb}\right)=F^{a}_{c}\,t_{ij}^{cb}-F^{b}_{c}\,t_{ij}^{ca}=F^{a}_{c}\,t_{ij}^{cb}+F^{b}_{c}\,t_{ij}^{ac}=\left(\zeta_{a}+\zeta_{b}\right)t_{ij}^{ab}

(45)

and

-P_{ij}\left(F_{i}^{k}\,t_{kj}^{ab}\right)=-\left(\zeta_{i}+\zeta_{j}\right)t_{ij}^{ab}.

(46)

One should note that in these terms, the indices on the eigenvalues $\zeta$ are not summed. These terms provide the denominators in perturbation theory so they are normally grouped with the driver terms. We use $\zeta$ for the eigenvalues of $F$ to emphasize that these eigenvalues are not the eigenvalues of the Fock operator, and contain a dependence on the $t$ amplitudes.

•

In the double permutation, we find all the ring (ph) contractions. There are eight terms in total:

\begin{split}P_{ij}P_{ab}\left[\left(v_{ic}^{ak}+\frac{1}{2}\,t_{il}^{ad}\,v_{cd}^{kl}\right)t_{kj}^{cb}\right]&=v_{ic}^{ak}\,t_{kj}^{cb}+\frac{1}{2}\,t_{il}^{ad}\,v_{cd}^{kl}\,t_{kj}^{cb}+v_{jc}^{bk}\,t_{ki}^{ca}+\frac{1}{2}\,t_{jl}^{bd}\,v_{cd}^{kl}\,t_{ki}^{ca}\\ &-v_{ic}^{bk}\,t_{kj}^{ca}-\frac{1}{2}\,t_{il}^{bd}\,v_{cd}^{kl}\,t_{kj}^{ca}-v_{jc}^{ak}\,t_{ki}^{cb}-\frac{1}{2}\,t_{jl}^{ad}\,v_{cd}^{kl}\,t_{ki}^{cb}.\end{split}

(47)

The two quadratic terms in the top row are identical as are the two quadratic terms in the bottom row, and using antisymmetry of $v$ and $t$ , we can simplify this slightly to

\begin{split}P_{ij}P_{ab}\left[\left(v_{ic}^{ak}+\frac{1}{2}\,t_{il}^{ad}\,v_{cd}^{kl}\right)t_{kj}^{cb}\right]&=v_{ic}^{ak}\,t_{kj}^{cb}+v_{jc}^{bk}\,t_{ki}^{ca}+t_{ki}^{ca}\,v_{cd}^{kl}\,t_{lj}^{db}\\ &-v_{ic}^{kb}\,t_{kj}^{ac}-v_{jc}^{ka}\,t_{ki}^{bc}+t_{ki}^{bc}\,v_{dc}^{kl}\,t_{lj}^{ad}.\end{split}

(48)

The terms in the first row are ring terms and are included in particle-hole RPA; those in the second row also involve particle-hole contractions, but with summation (and external) indices crossed, so we refer to these as crossed ring terms. Note that including the ring terms but excluding the crossed rings, thereby including only a portion of the antisymmetry term $P_{ij}P_{ab}[\ldots]$ , breaks the antisymmetry of the amplitude equations and therefore of the $T_{2}$ amplitudes.

III.1 Ring-CCD

Here, we merely summarize the results of Ref. Scuseria et al., 2008. We collect the driving term and the ring terms (but not the crossed ring terms) in the ring-CCD equation:


$\displaystyle 0$	$\displaystyle=v_{ij}^{ab}-P_{ij}\left(F_{i}^{k}\,t_{kj}^{ab}\right)+P_{ab}\left(F^{a}_{c}\,t_{ij}^{cb}\right)+v_{ic}^{ak}\,t_{kj}^{cb}+v_{cj}^{kb}\,t_{ik}^{ac}+t_{ik}^{ac}\,v^{kl}_{cd}\,t_{lj}^{db}$	(49a)
	$\displaystyle=v_{ij}^{ab}+\left(F^{a}_{c}\,\delta_{i}^{k}-F_{i}^{k}\,\delta^{a}_{c}\right)t_{kj}^{cb}+\left(F^{b}_{c}\,\delta_{j}^{k}-F_{j}^{k}\,\delta^{b}_{c}\right)t_{ik}^{ac}+v_{ic}^{ak}\,t_{kj}^{cb}+v_{kj}^{cb}\,t_{ik}^{ac}+t_{ik}^{ac}\,v^{kl}_{cd}\,t_{lj}^{db}.$	(49b)

The resulting $T_{2}$ amplitudes are not antisymmetric but do retain the bosonic symmetry $t_{ij}^{ab}=t_{ji}^{ba}$ . We can simplify the amplitudes equations by using the $\bm{A}$ and $\bm{B}$ matrices of particle-hole RPA, which in this notation are


$\displaystyle A_{ia,jb}$	$\displaystyle=F^{b}_{a}\,\delta^{i}_{j}-F^{i}_{j}\,\delta^{b}_{a}+v^{ib}_{aj},$	(50a)
$\displaystyle B_{ia,jb}$	$\displaystyle=v^{ij}_{ab}.$	(50b)

Note that we have used the Brueckner Hamiltonian $\bm{F}$ in defining these matrices, though one can instead use the Hartree-Fock Hamiltonian $\bm{f}$ . It is apparent that the ring-CCD equation can be written as

\bm{B}^{\star}+\bm{A}^{\star}\,\bm{T}+\bm{T}\,\bm{A}+\bm{T}\,\bm{B}\,\bm{T}=\bm{0};

(51)

replacing the Brueckner Hamiltonian with the Hartree-Fock Hamiltonian corresponds to discarding the mosaic terms in the ring-CCD equations.

On the other hand, the particle-hole RPA equations for the non-negative excitation energies are

\begin{pmatrix}\bm{A}&\bm{B}\\ -\bm{B}^{\star}&-\bm{A}^{\star}\end{pmatrix}\begin{pmatrix}\bm{X}\\ \bm{Y}\end{pmatrix}=\begin{pmatrix}\bm{X}\\ \bm{Y}\end{pmatrix}\bm{\omega}.

(52)

We have assumed that we are in a physically meaningful case where $\bm{\omega}$ is real, so that we can always choose non-negative $\bm{\omega}$ . We emphasize again that this RPA is a bosonic mean-field problem where the Hermitian $\bm{A}$ and symmetric $\bm{B}$ play the role of Fock and pairing fields, respectively.

The equivalence between ring-CCD and particle-hole RPA is most simply established when $\bm{X}$ is invertible, so that we can define $\bm{T}=\bm{Y}\,\bm{X}^{-1}$ . In that case, the RPA eigenvalue problem can be used to derive


$\displaystyle\bm{A}+\bm{B}\,\bm{T}$	$\displaystyle=\bm{X}\,\bm{\omega}\,\bm{X}^{-1}=\bm{R},$	(53a)
$\displaystyle\bm{B}^{\star}+\bm{A}^{\star}\,\bm{T}$	$\displaystyle=-\bm{T}\,\bm{R}.$	(53b)

Inserting the first equation into the second and rearranging yields

\bm{B}^{\star}+\bm{A}^{\star}\,\bm{T}+\bm{T}\,\bm{A}+\bm{T}\,\bm{B}\,\bm{T}=\bm{0}

(54)

so that from particle-hole RPA we can extract the amplitudes $\bm{T}$ which solve the ring-CCD equation. Moreover, the RPA correlation energy comes from the plasmon formula,

E_{c}^{\mathrm{RPA}}=\frac{1}{2}\,\mathrm{Tr}(\bm{\omega}-\bm{A})

(55)

while the coupled cluster correlation energy is just

E_{c}^{\mathrm{CCD}}=\frac{1}{4}\,t_{ij}^{ab}\,v^{ij}_{ab}=\frac{1}{4}\,\mathrm{Tr}(\bm{B}\,\bm{T}).

(56)

From Eqn. 53a, we see that

\mathrm{Tr}(\bm{\omega}-\bm{A})=\mathrm{Tr}(\bm{B}\,\bm{T}).

(57)

Thus, the ring-CCD and particle-hole RPA correlation energies differ by a factor of two. This reflects the fact that ring CCD is not a correct bosonization of the fermionic problem. The discrepancy in the correlation energy disappears for direct RPA, where we keep only Hartree terms in the interaction.

III.2 Ladder-CCD

We have seen that the ring-CCD problem is intimately connected to particle-hole RPA; here, we demonstrate that the ladder-CCD problem is analogously connected to particle-particle RPA. To the best of our knowledge, this connection has never been discussed in the literature.

The ladder-CCD equations are


$\displaystyle 0$	$\displaystyle=v_{ij}^{ab}-P_{ij}\left(F^{k}_{i}\,t_{kj}^{ab}\right)+P_{ab}\left(F^{a}_{c}\,t_{ij}^{cb}\right)+\frac{1}{2}\,t_{kl}^{ab}\,v^{kl}_{ij}+\frac{1}{2}\,v^{ab}_{cd}\,t_{ij}^{cd}+\frac{1}{4}\,t_{ij}^{cd}\,v_{cd}^{kl}\,t_{kl}^{ab}$	(58a)
	$\displaystyle=v_{ij}^{ab}+\left(F^{a}_{c}\,\delta^{b}_{d}+F^{b}_{d}\,\delta^{a}_{c}\right)t_{ij}^{cd}-\left(F^{k}_{i}\,\delta^{l}_{j}+F^{l}_{j}\,\delta^{k}_{i}\right)\,t_{kl}^{ab}+\frac{1}{2}\,t_{kl}^{ab}\,v^{kl}_{ij}+\frac{1}{2}\,v^{ab}_{cd}\,t_{ij}^{cd}+\frac{1}{4}\,t_{ij}^{cd}\,v_{cd}^{kl}\,t_{kl}^{ab}.$	(58b)

We can express this in terms of the matrices $\bar{\bm{B}}$ , $\bm{C}$ , and $\bm{D}$ of particle-particle RPA, which in this notation are


$\displaystyle\bar{B}_{ab,ij}$	$\displaystyle=v^{ab}_{ij}$	(59a)
$\displaystyle C_{ab,cd}$	$\displaystyle=\left(F^{a}_{c}\,\delta^{b}_{d}+F^{b}_{d}\,\delta^{a}_{c}\right)+v^{ab}_{cd},$	(59b)
$\displaystyle D_{ij,kl}$	$\displaystyle=-\left(F^{k}_{i}\,\delta^{l}_{j}+F^{l}_{j}\,\delta^{k}_{i}\right)+v^{kl}_{ij},$	(59c)

where we recall that $\bar{\bm{B}}$ is rectangular and $\bm{C}$ and $\bm{D}$ are Hermitian. Note that while $\bar{\bm{B}}$ and $\bm{B}$ have the same matrix elements, they are organized into bosonic composite indices differently.

With these definitions in hand, the ladder-CCD equations become

\bar{\bm{B}}+\bm{C}\,\bm{T}+\bm{T}\,\bm{D}^{\star}+\bm{T}\,\bar{\bm{B}}^{\dagger}\,\bm{T}=\bm{0},

(60)

where the composite indices restrict $a<b$ and $i<j$ , yielding the needed factors of $\tfrac{1}{2}$ and $\tfrac{1}{4}$ . Here, we have clearly defined $\bm{T}$ as being a $vv\times oo$ matrix. Equivalently, we could have written

\bar{\bm{B}}^{\dagger}+\tilde{\bm{T}}\,\bm{C}+\bm{D}^{\star}\,\tilde{\bm{T}}+\tilde{\bm{T}}\,\bar{\bm{B}}\,\tilde{\bm{T}}=\bm{0},

(61)

where $\tilde{\bm{T}}=\bm{T}^{\dagger}$ is $oo\times vv$ .

Recall that the particle-particle RPA problem is

\begin{pmatrix}\bm{C}&-\bar{\bm{B}}\\ \bar{\bm{B}}^{\dagger}&-\bm{D}^{\star}\end{pmatrix}\begin{pmatrix}\bm{X}_{1}&\bm{Y}_{2}\\ \bm{Y}_{1}&\bm{X}_{2}\end{pmatrix}=\begin{pmatrix}\bm{X}_{1}&\bm{Y}_{2}\\ \bm{Y}_{1}&\bm{X}_{2}\end{pmatrix}\begin{pmatrix}\bm{\Omega_{1}}&\bm{0}\\ \bm{0}&\bm{\Omega}_{2}\end{pmatrix},

(62)

and that $\bm{\Omega}_{1}$ is positive while $\bm{\Omega}_{2}$ is negative. From particle-particle RPA, we can write two Riccati equations. One writes $\bm{T}_{1}=-\bm{Y}_{1}\,\bm{X}_{1}^{-1}$ and the other $\bm{T}_{2}=-\bm{Y}_{2}\,\bm{X}_{2}^{-1}$ ; these matrices $\bm{T}_{1}$ and $\bm{T}_{2}$ are of dimension $vv\times oo$ and $oo\times vv$ , repectively.

The Riccati equation for $\bm{T}_{1}$ follows from the particle-particle RPA problem for $\bm{X}_{1}$ and $\bm{Y}_{1}$ , which yields


$\displaystyle\bm{C}+\bar{\bm{B}}\,\bm{T}_{1}$	$\displaystyle=\bm{X}_{1}\,\bm{\Omega}_{1}\,\bm{X}_{1}^{-1}=\bm{R}_{1},$	(63a)
$\displaystyle\bar{\bm{B}}^{\dagger}+\bm{D}^{\star}\,\bm{T}_{1}$	$\displaystyle=-\bm{T}_{1}\,\bm{R}_{1},$	(63b)

from which one extracts


$\displaystyle\bm{0}$	$\displaystyle=\bar{\bm{B}}^{\dagger}+\bm{D}^{\star}\,\bm{T}_{1}+\bm{T}_{1}\,\bm{C}+\bm{T}_{1}\,\bar{\bm{B}}\,\bm{T}_{1},$	(64a)
$\displaystyle\mathrm{Tr}(\bar{\bm{B}}\,\bm{T}_{1})$	$\displaystyle=\mathrm{Tr}(\bm{\Omega}_{1}-\bm{C}).$	(64b)

If we instead use the particle-particle RPA problem for $\bm{X}_{2}$ and $\bm{Y}_{2}$ , we get


$\displaystyle-\bar{\bm{B}}^{\dagger}\,\bm{T}_{2}-\bm{D}^{\star}$	$\displaystyle=\bm{X}_{2}\,\bm{\Omega}_{2}\,\bm{X}_{2}^{-1}=\bm{R}_{2},$	(65a)
$\displaystyle-\bm{C}\,\bm{T}_{2}-\bar{\bm{B}}$	$\displaystyle=-\bm{T}_{2}\,\bm{R}_{2}$	(65b)

which imply that


$\displaystyle\bm{0}$	$\displaystyle=\bar{\bm{B}}+\bm{C}\,\bm{T}_{2}+\bm{T}_{2}\,\bm{D}^{\star}+\bm{T}_{2}\,\bar{\bm{B}}^{\dagger}\,\bm{T}_{2},$	(66a)
$\displaystyle\mathrm{Tr}(\bar{\bm{B}}^{\dagger}\,\bm{T}_{2})$	$\displaystyle=-\mathrm{Tr}(\bm{\Omega}_{2}+\bm{D}^{\star}).$	(66b)

Clearly, $\bm{T}_{2}=\bm{T}_{1}^{\dagger}$ .

Finally, the ladder-CCD correlation energy is

\begin{split}E_{c}&=\mathrm{Tr}(\bar{\bm{B}}\,\bm{T}^{\dagger})=\mathrm{Tr}(\bar{\bm{B}}^{\dagger}\,\bm{T})=\mathrm{Tr}(\bm{\Omega}_{1}-\bm{C})=-\mathrm{Tr}(\bm{\Omega}_{2}+\bm{D}^{\star})\\ &=\frac{1}{2}\mathrm{Tr}\left(\bm{\Omega}_{1}-\bm{\Omega}_{2}-\bm{C}-\bm{D^{\star}}\right)\end{split}

(67)

which is exactly the result from particle-particle RPA. The formal equivalence proven here is validated numerically by the calculations we discuss later.

III.3 A Third Channel: Crossed-Ring–CCD

Comparing the ring- and ladder-CCD equations to the full CCD equations reveals that we have used the driving terms twice, and have excluded the crossed ring terms entirely. We can take these remaining pieces and form a third sort of channel. With this channel we can associate a Riccati equation and form an RPA-like problem; this RPA-like problem, however, has no clear physical significance (unlike particle-hole and particle-particle RPA).

Because we have counted the driving terms twice, we would like the crossed-ring–CCD to take the driving term with a minus sign so that the CCD equations for the three channels add to give the regular CCD equations; again, all of this is subject to caveats with regards to the mosaic terms, which should be included in decomposing the CCD equations into these three channels but which are not present in the typical RPA approach. This choice turns out not to be associated with a symplectic RPA-like matrix. Rather, we must include the interaction $v_{ij}^{ab}$ with a positive sign and the remainder of the driver term (the terms giving rise to orbital energy denominators, in other words) we can safely leave with a negative sign. This gives us

v_{ij}^{ab}-F^{a}_{c}\,t_{ij}^{cb}-F^{b}_{c}\,t_{ij}^{ac}+F^{k}_{i}\,t_{kj}^{ab}+F_{j}^{k}\,t_{ik}^{ab}-v_{ic}^{kb}\,t_{kj}^{ac}-v_{jc}^{ka}\,t_{ik}^{cb}-t_{ik}^{cb}\,v^{kl}_{cd}\,t_{lj}^{ad}=0.

(68)

This can be fruitfully rewritten as

-v_{ij}^{ab}+\left(F^{a}_{c}\,\delta^{k}_{j}-F^{k}_{j}\,\delta^{a}_{c}+v_{jc}^{ka}\right)t_{ik}^{cb}+\left(F^{b}_{c}\,\delta^{k}_{i}-F^{k}_{i}\,\delta^{b}_{c}+v_{ic}^{kb}\right)t_{kj}^{ac}-t_{ik}^{cb}\,v^{kl}_{dc}\,t_{lj}^{ad}=0.

(69)

This time, we want to define composite indices as $ib$ and $ja$ ; this prevents us from simply adding the ring and crossed-ring equations. We can define


$\displaystyle\tilde{A}_{kc,ja}$	$\displaystyle=F^{a}_{c}\,\delta^{k}_{j}-F^{k}_{j}\,\delta^{a}_{c}+v_{jc}^{ka}$	(70a)
$\displaystyle\tilde{B}_{kc,ld}$	$\displaystyle=v^{kl}_{dc}$	(70b)

where $\tilde{\bm{A}}$ and $\tilde{\bm{B}}$ are closely related to the $\bm{A}$ and $\bm{B}$ matrices of particle-hole RPA, but differ by the sign of the two-electron integral. In terms of these newest quantities, the crossed-ring–CCD looks like

\bm{0}=-\tilde{\bm{B}}^{\star}+\tilde{\bm{T}}\,\tilde{\bm{A}}+\tilde{\bm{A}}^{\star}\,\tilde{\bm{T}}-\tilde{\bm{T}}\,\tilde{\bm{B}}\,\tilde{\bm{T}}.

(71)

From our discussion of particle-hole RPA, it should be clear that this is the Riccati equation corresponding to the symplectic eigenvalue problem

\begin{pmatrix}\tilde{\bm{A}}&-\tilde{\bm{B}}\\ \tilde{\bm{B}}^{\star}&-\tilde{\bm{A}}^{\star}\end{pmatrix}\begin{pmatrix}\bm{X}\\ \bm{Y}\end{pmatrix}=\begin{pmatrix}\bm{X}\\ \bm{Y}\end{pmatrix}\bm{\varpi}.

(72)

III.4 Combining the Three Channels

Returning to the Brueckner CCD amplitude equations, we see that we can write them

0=v+\left(\mathrm{Ring}-v\right)+\left(\mathrm{Ladder}-v\right)+\left(\textrm{Crossed-Ring}-v\right).

(73)

We emphasize that we cannot simply add the Riccati equations for the three channels, as the bosonic composite external indices are formed from the fermionic indices in three different ways. We also emphasize that including the mosaic terms in defining an effective one-body Hamiltonian is necessary to leave this simple form. Finally, note that the three Riccati equations, corresponding to three different RPA-like problems, are tied together by the requirement that the amplitudes $t$ are the same – these amplitudes, in other words, force the channels to interact. This is made possible by the inclusion of the crossed-ring terms, which provide the necessary antisymmetrization that the ring channel lacks.

The patching of these three problems as a conceptual tool for understanding CCD has never been discussed in the literature, to the best of our knowledge. From the analysis above, we see that the CC equations can be interpreted as the sum of three quadratic bosonic problems: rings, ladders, and crossed-rings with a renormalized one-body Hamiltonian ( $F$ ) and the regular Coulomb two-body interaction. Each of the channels contract the CC amplitudes in a different manner. In this sense, we could say that CC theory is as a correct bosonization of fermion excitations because it yields fully antisymmetrized excitation amplitudes fulfilling Pauli’s principle, together with a well-defined wavefunction and corresponding density matrices.

IV Results

In this section, we provide a few numerical results, showing the relative importance of the particle-particle and particle-hole channels of RPA, as well as the importance of the various contributions to the CCD equations. Variants of CCD will be identified by whether they include ladder terms (“l”), ring terms (“r”), and mosaic terms (“m”); thus, ring-CCD in this notation is rCCD and ladder-CCD is l-CCD. Because they do not appear in RPA, we do not include crossed-ring terms in this section except in the form of the full CCD. We shall see later that while the crossed-rings are vital for restoring the full fermionic character of the CCD wavefunction, they can also unbalance the CC equations when other classes of diagrams are omitted.

All results were generated using an in-house program. As this paper is not intended to generate benchmark data or even comparisons with experiment, but rather to discuss qualitative features, we will use small basis sets and not worry about basis set incompleteness error.

We begin our discussion with the dissociation of H₂, which is of course paradigmatic in quantum chemistry. As is well known, the restricted Hartree-Fock (RHF) solution goes to much too high an energy, as a result of contamination from ionic dissociation fragments. This is remedied (energetically) by unrestricted Hartree-Fock (UHF), which dissociates correctly at the cost of broken spatial and spin symmetries. There is thus an instability from RHF to UHF past the Coulson-Fischer point, and beyond this point particle-hole RPA based on the RHF reference yields unphysical complex correlation energies.

Refer to caption — Figure 1: Dissociation curve of H₂ in the cc-pvdz basis set, obtained with various flavors of RPA, compared to CCD, ring-CCD, and ladder-CCD. Top panel: RHF reference. Bottom panel: UHF reference.

The RPA dissociation of H₂ is shown in Fig. 1. There are several key features here we wish to point out. First, the ladder-CCD and the particle-particle RPA energies are indeed identical, bearing out our anayltic proof earlier. While particle-particle RPA undercorrelates, particle-hole RPA overcorrelates. Since quasiparticle RPA includes both, it overcorrelates even more. Near the Coulson-Fischer point, particle-hole RPA is particularly bad, with a cusp at the Coulson-Fischer point; this behavior is inherited by quasiparticle RPA. If one uses an RHF reference instead of a UHF reference, the particle-hole and quasiparticle RPA energies become complex. Though the particle-particle RPA remains well behaved, its undercorrelation is greatly exaggerated.

In Fig. 2, we show results for variants of CCD including selective terms. Let us start with results based on the RHF reference. In this case, as is well-known, the ring-CCD does not converge past the Coulson-Fischer point. Our results are nevertheless sufficient to illustrate that ladder-CCD undercorrelates while ring-CCD overcorrelates. Including both ladders and rings undercorrelates, but for a sufficiently stretched bond, we were unable to converge the equations. Adding the mosaic terms has, in general, a relatively small effect, but note that they cure the convergence difficulties we face due to the inclusion of the ring diagrams without corresponding inclusion of the crossed rings.

The story is qualitatively similar using the UHF reference. Again, the ladder-CCD undercorrelates while the ring-CCD overcorrelates. Including both ladders and rings qualitatively resembles ladder-CCD. In this case, the ring-CCD equations can be converged to dissociation, but the curve near the Coulson-Fischer point remains pathological. Again, the mosaic terms have a relatively small effect for the most part, but cure the worst of the pathologies of ring-CCD. We note finally that CCD itself has a shoulder at the Coulson-Fischer point, which would be cured by the inclusion of single excitations. Equivalently, we could iterate the BD equations, in which case the mosaic terms would be naturally taken care of and the CCD would become exact.

We next turn to the dissociation of LiH, for which we use the 6-31G** basis set. In Fig. 3 we show UHF-based RPA dissociation curves, while Fig. 4 shows both RHF- and UHF-based CCD dissociations. The main details remain the same as in H₂. Ring-CCD overcorrelates significantly, and particle-hole accordingly overcorrelates even more badly. In contrast, ladder-CCD = particle-particle RPA undercorrelates. Including both ladders and rings in quasiparticle RPA overcorrelates even more badly than does particle-hole RPA (not shown), while including them both in a coupled cluster approach improves but does not fully cure the undercorrelation from ladder-CCD. Starting from the RHF reference, ring-CCD does not converge past the Coulson-Fischer point and lr-CCD eventually stops converging as well. These convergence problems are cured by including the mosaic terms. Starting from the UHF reference, ring-CCD (and therefore particle-hole RPA) exhibits pathological behavior in the vicinity of the Coulson-Fischer point, which is again cured by adding the mosaic terms. In all other cases, the mosaic terms have essentially negligible effects.

Lastly, we consider the dissociation of N₂, which we have examined in the cc-pVDZ basis set with $d$ functions removed for ease of convergence, and since we remain interested only in the qualitative picture. Results are shown in Fig. 5, using a broken-symmetry reference. Yet again, the particle-hole RPA overcorrelates, as does the ring-CCD, while the particle-particle RPA undercorrelates. As usual, the particle-hole RPA and ring-CCD display especially severe problems near the Coulson-Fischer point. The mosaic terms are again fairly unimportant, except in the case of ring-CCD where they offer a large improvement. Including both ladders and rings in the CCD is superior to including just one or the other.

V Discussion

It should be clear from previous sections that there is a recurrying theme in all of these approaches, which is the attempt to treat the ground state with a fermionic wavefunction which is to be correlated via bosonic excitations. This seems like a natural way of using fermion and boson mean-field methods. The main difficulty is that composite fermions are not actually bosonic, so whatever sort of bosonization one undertakes can only be approximate when the composite indices are never broken (as in the various flavors of RPA). A secondary difficulty is that bosonizing the excitations, which are $\mathcal{O}(M^{2})$ in number where $M$ is the number of basis functions, leads to an $\mathcal{O}(M^{6})$ bosonic mean-field treatment, so that the computational scaling does not improve upon traditional CCD without further approximation. This is true both for particle-hole and for particle-particle RPA. Note that in practice, the scaling of RPA can often be considerably reduced through techniques such as Cholesky decomposition and variants.Eshuis et al. (2012)

The traditional CCD equations achieve a sort of correct bosonization by including not only the particle-particle (ladder) and particle-hole (ring) channels, but adding crossed-ring terms as well. Insisting, additionally, that the amplitudes obtained from the corresponding symplectic eigenvalue problems all be identical once the composite fermion indices are interpreted forces the three channels to interact, and in this way overcomes the overcorrelation of quasiparticle RPA while correctly enforcing fermionic antisymmetry.

From one perspective, the main function of the crossed-ring terms in the CCD equations is simply to guarantee this fermionic antisymmetry. Including the crossed-ring terms is not, however, the only way to impose this constraint. Second-order screened exchangeFreeman (1977) (SOSEX) does this for the case of ring-CCD. In SOSEX, one first solves the ring-CCD equations, then antisymmetrizes the resulting amplitudes before calculating the correlation energy (using, it should be noted, the correct coupled cluster factor of $\tfrac{1}{4}$ ). This reduces the ring-CCD correlation energy substantially, preventing the dramatic overcorrelation which we have seen. One might think that a better approach might be to include both rings and crossed-ring in a coupled cluster approach. This, however, is unsuccessful. Including only rings and crossed-rings in a CCD framework, we encounter severe convergence difficulties; when we can converge the results, they are quite poor. Adding the mosaic terms as well alleviates these convergence difficulties, but overcorrelates terribly, as seen in Fig. 6.

While the crossed-ring terms can be naturally incorporated into the CCD framework, their origin in RPA of all flavors is obscure. In fact, they seem entirely artificial from the RPA perspective. A proper RPA-style mechanism by which these antisymmetrizing pieces can be included would most likely be quite valuable.

An important point to bear in mind in quasiparticle RPA is that the 13 and 31 blocks of the Hamiltonian do not contribute at all. We can imagine a Brueckner-style renormalization such as

\tilde{H}_{pq}^{11}=H_{pq}^{11}+\left[H^{31}_{prst}\left(\bm{\mathcal{Y}}\,\bm{\mathcal{X}}^{-1}\right)_{stqr}+\textrm{h.c.}\right]

(74)

where, in analogy with the Brueckner effective Hamiltonian, the effects of $H^{31}$ and $H^{13}$ are incorporated iteratively by changing the reference quasiparticle determinants orbitals and energies. One might hope that the dressed Fock and pairing fields caused by this Brueckner-style renormalization would induce number fluctuations in the reference determinant. Were this to happen, the particle-particle and particle-hole channels would interact. This might cure some of the overcorrelation endemic to quasiparticle RPA. Recall that for a standard repulsive two-body interaction, the HFB method is entirely equivalent to Hartree-Fock.

Overall, we recommend a careful combination of particle-particle and particle-hole RPA if one wishes to work within an RPA framework. The two methods are in some sense complementary, and each describes different physics. Particle-hole RPA, for example, works quite well for long-range electron-electron interactions, and successfully incorporates van der Waals binding.Dobson et al. (2006) On the other hand, from the form of the bosonic excitation operators, we see that particle-particle RPA should be suitable for the description of charge fluctuations which are beyond the scope of conventional particle-hole RPA. These recommendations are particularly important in efforts to incorporate RPA correlation effects into density functional theory.¹¹1While this paper was about to be submitted, we learned about an effort to incorporate particle-particle RPA correlation in DFT. See H. van Aggelen, Y. Yang, and W. Yang, (2013) arXiv:1306.4957, and D. Peng, S. N. Steinmann, H. van Aggelen, and W. Yang, (2013) arXiv:1306.5638 We have seen, moreover, that the CCD result often lies between ladder-CCD and ring-CCD, and thus between particle-hole and particle-particle RPA. Finally, what we have christened the mosaic terms are also apparently quite important, and should perhaps be included in post–Hartree-Fock RPA schemes.

VI Acknowledgments

This work was supported by the National Science Foundation CHE-1110884 and the Welch Foundation (C-0036).

References

Scuseria et al. (2008) G. E. Scuseria, T. M. Henderson, and D. C. Sorensen, J. Chem. Phys. 129, 231101 (2008).
Ren et al. (2012) X. Ren, P. Rinke, J. Christian, and M. Scheffler, J. Mater. Sci. 47, 7447 (2012).
Eshuis et al. (2012) H. Eshuis, J. Bates, and F. Furche, Theor. Chem. Acc. 131, 1084 (2012).
Ring and Schuck (1980) P. Ring and P. Schuck, The Nuclear Many-Body Problem (Springer-Verlag, Berlin, 1980).
Blaizot and Ripka (1985) J.-P. Blaizot and G. Ripka, Quantum Theory of Finite Systems (The MIT Press, Cambridge, MA, 1985).
Bach et al. (1994) V. Bach, E. H. Lieb, and J. P. Solovej, J. Stat. Phys. 76, 3 (1994).
Scuseria and Tsuchimochi (2009) G. E. Scuseria and T. Tsuchimochi, J. Chem. Phys. 131, 164119 (2009).
Colpa (1986) J. H. P. Colpa, Physica A 134, 377 (1986).
Bartlett and Shavitt (2009) R. J. Bartlett and I. Shavitt, Many-Body Methods in Chemistry and Physics (Cambridge University Press, New York, 2009).
Scuseria (1995) G. E. Scuseria, Int. J. Quantum Chem. 55, 165 (1995).
Freeman (1977) D. L. Freeman, Phys. Rev. B 15, 5512 (1977).
Dobson et al. (2006) J. F. Dobson, A. White, and A. Rubio, Phys. Rev. Lett. 96, 073201 (2006).