Beyond a Richardson-Gaudin mean-field: Slater-Condon rules and perturbation theory

Weakly correlated electrons are well described by Kohn-Sham density functional theory (KS-DFT) and coupled-cluster (CC) theory.^{1, 2, 3, 4, 5, 6} The qualitative behaviour is approximated by a Slater determinant and corrected by an approximate functional (KS-DFT) or by developing the wavefunction (CC). Both approaches are succesful as the corrections are small. In either case the approach is well understood, systematic improvement is possible, and results are easily computable by the non-expert with software packages.

For strongly correlated electrons there is generally no obvious Slater determinant reference. Many Slater determinants are required so that a qualitative mean-field description is difficult. The complete active space self-consistent field (CASSCF)^{7, 8, 9, 10} is successful if the important orbitals are few enough in number and easily identified. Approximate CASSCF solvers^{11, 12, 13, 14, 15, 16, 17, 18, 19} have been developed to treat larger systems, for which the density matrix renormalization group (DMRG)^{20, 21, 22, 23, 24, 25} has also been applied with great success. These methods target the wavefunction directly in a computationally efficient manner without making reference to a mean-field picture.

It is understood that pairs of electrons are a more appropriate starting point for strongly correlated electrons.²⁶ The antisymmetrized product of interacting geminals (APIG) provides an excellent description of paired electrons for an intractable cost.^{27, 28, 29, 30, 31, 32, 33} These results are obtainable at mean-field cost through the antisymmetric product of 1-reference orbital geminals (AP1roG),³⁴ or eqiuvalently, pair-coupled-cluster doubles (pCCD).^{35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53} AP1roG/pCCD shows excellent results for the ground state energy of systems with no unpaired electrons. Generalizations to systems with unpaired electrons, or broken pairs of electrons, is however not obvious.

The reduced Bardeen-Cooper-Schrieffer (BCS)^{54, 55, 56} Hamiltonian is an exactly solvable model whose eigenvectors are products of pairs of electrons. These states were first obtained by Richardson^{57, 58, 59} and generalized by Gaudin⁶⁰ and are thus referred to as Richardson-Gaudin (RG) states. As eigenvectors of a physical model, RG states form a basis for the Hilbert space. It has been shown for strongly correlated model systems that a single RG state is a good starting point, much like a Slater determinant for weakly correlated systems.^{61, 62, 63}

RG states are feasible as their reduced density matrix (RDM) elements are cheap to compute,^{64, 65, 66} and their couplings with excited states decay rapidly with excitation level.⁶⁷ Previous constructions of the transition density matrix (TDM) elements to perform a configuration interaction with singles and doubles (CISD) of RG states scaled very poorly, which is a situation remedied herein. The purpose of this contribution is to nail down the details of using RG states for fully paired electrons, the so-called seniority-zero sector. This is meant as a stepping stone to RG states with unpaired electrons, which are much more complicated and will be treated in a future contribution.

This manuscript is organized as follows. First, a brief overview of RG states is given, discussing in particular their density matrix elements and low-lying excited states. Optimal expressions for the TDM elements are obtained using an approach developed by Chen and Scuseria for Wick’s theorem applied to Hartree-Fock-Bogoliubov states.⁶⁸ Second, contributions of the individual RG states to the RGCISD wavefunction are studied, leading to analogues of the Slater-Condon rules for RG states. Finally, perturbation theory is shown to give comparable results to RGCISD at a fraction of the cost.

This section provides a brief outline of Richardson-Gaudin states, keeping the details to a minimum. For a complete description see ref.⁶⁹ There are always $N$ sites or spatial orbitals labelled with indices $i,j,k,l$, and $M$ pairs labelled with indices $a,b$.

RG states are built from pairs of electrons, which are elementary representations of the Lie algebra su(2). For each spatial orbital there are three operators

$S^{+}_{i}$ creates a pair in spatial orbital $i$, $S^{-}_{i}$ removes a pair in spatial orbital $i$, while $S^{z}_{i}$ gives $+\frac{1}{2}$ for a pair-occupied spatial orbital, $-\frac{1}{2}$ for an empty spatial orbital, and 0 for a singly-occupied orbital. These operators have the usual structure of the spin operators of su(2)

It is convenient to use $\hat{n}_{i}=2S^{z}_{i}+1$, which counts the number of electrons in spatial orbital $i$. RG states are the eigenvectors of the reduced BCS Hamiltonian

They are products of RG pairs

which depend on a set of complex numbers usually called rapidities. The RG pair (5) is a linear combination of all one-pair states weighted by the difference of the rapidity and the energy $\varepsilon_{i}$ of the spatial orbital $i$. Rapidities are not free variables. In particular for the RG state

to be an eigenvector of the reduced BCS Hamiltonian, the rapidities $\{u\}$ must be a solution of the nonlinear equations

first derived by Richardson, and are hence known as Richardson’s equations. Each RG state is defined by a complete set of rapidities that together solve Richardson’s equations. Thus, strictly speaking there are no rapidities shared between eigenvectors of the reduced BCS Hamiltonian.

Numerically solving Richardson’s equations is possible, but extra care must be taken near the critical points where a rapidity coincides with a single-particle energy.^{70, 71, 72, 73} It is far easier to work with so-called eigenvalue-based variables (EBV)

which are solutions of the non-linear equations

Taking the sum of Richardson’s equations, for each $a$, one finds that the EBV satisfy the sum rule

These equations are easy to solve numerically and do not have troublesome critical points. For the complete details of solving the EBV equations (9) see refs.^{74, 75, 62} It is convenient to use rapidities when looking at the geminal coefficients in each pair (5), but numerically they are not productive to compute. RG pairs and states have clear definitions in terms of rapidities but do not in terms of EBV. Thus, we continue to label RG states with rapidities even if we do not compute them numerically.

When $g=0$, the reduced BCS Hamiltonian is non-interacting and the EBV equations (9) decouple completely

whose solutions are $M$ EBV equal to 2 and $N-M$ EBV equal to 0. These RG states are Slater determinants whose pair-occupied orbitals correspond to the particular EBV equal to 2. The remarkable fact is that as $g$ grows in either direction, the RG states are connected uniquely to one solution at $g=0$.^{76, 77} RG states can thus be labelled by the bitstring of the Slater determinant they represent at $g=0$ even though at finite $g$ they are not Slater determinants. The ground state of the reduced BCS Hamiltonian is always labelled by $M$ 1s followed by $N-M$ 0s, while the highest excited state is always labelled by $N-M$ 0s followed by $M$ 1s. Other RG states will cross, but these are always strict crossings and not avoided crossings.^{78, 79, 80}

In ref.⁶⁷ it was seen that the RG ground state $1^{M}0^{N-M}$ coupled appreciably only with its single and double pair excitations. This will now be established and justified for the RG state $\ket{M}$. Previous studies of small strongly correlated systems of Hydrogen atoms were very well described with a configuration interaction (CI) of RG states.^{63, 69} These systems do not exist in reality: a linear chain of equidistant Hydrogen atoms would undergo a Peirls distortion and become a set of independent H₂ molecules. However, these systems are studied as they maximize the effects of strong correlation while being small enough to treat exactly.

RG states have no unpaired electrons: they have zero seniority. They are thus approximations to a CI of all Slater determinants with no unpaired electrons, the so-called doubly-occupied configuration interaction (DOCI). DOCI is not exact though it has been shown that classifying Slater determinants by seniority leads provides a systematic hierarchy for strongly correlated systems, even multiple bond-breaking processes like N₂.⁹² DOCI depends on the choice of orbitals: a Slater determinant that has no unpaired electrons in a given set of orbitals will have unpaired electrons in a different set of orbitals. It is therefore necessary to use orbital-optimized (OO-)DOCI. The OO-DOCI results for equidistant linear H₄ and H₈ were computed in the basis STO-6G in ref.⁶¹ For these systems, a CI of the RG state $\ket{M}$ along with its singles and doubles, RGCISD, is numerically indistinguishable from OO-DOCI.⁶⁹ The CI coefficients of each RG state contributing to RGCISD will now be studied to demonstrate that the expansion is short, and that the couplings between the RG states follow simple rules.

CI coefficients of the RG states contributing to the ground state of equidistant linear H₄ are shown as a function of the H – H distance in Figure 2. For this system, RGCISD and DOCI have the same cardinality and thus the two treatments are equivalent. Raw energy curves are not informative, as RGCISD is indistinguishable from OO-DOCI, but are shown in ref.⁶⁹ As expected, the coefficient of the reference $\ket{M}$ completely dominates the expansion, only differing from 1 near the “equilibrium” geometry at $r=1.65$ bohr. The transfers $\ket{M^{2}_{1}}$ and $\ket{M^{1}_{2}}$ provide the interpair weak correlation and drop off very quickly once $r$ is greater than $r_{C}=3.1$ bohr. The contributions from the swaps $\ket{M^{1}_{1}}$ and $\ket{M^{2}_{2}}$ are weak and centred at zero. Finally, the double swap $\ket{M^{12}_{12}}$ gives a weak contribution near equilibrium, but grows and peaks at $r=2.95$ bohr, which is just before $r_{C}$, before decaying back to zero.

The same behaviour is seen in the CI coefficients of the RG states for linear H₈. Here RGCISD is not DOCI, but it has been seen previously that the difference is on the order of $10^{-10}\;E_{h}$.⁶⁹ Sums of absolute values of CI coefficients are plotted rather than individual coefficients for two reasons. First, to keep the plots clean and legible. Second, it is near-impossible to get clean continuous curves as even a very small discontinuity in any of the parameters, in particular the orbital coefficients, will cause a large difference in the individual CI coeffcients but no difference in the energy. Sums of absolute values of CI coefficients are invariant to such discontinuities while retaining the physical meaning. The expansion is completely dominated by the reference $\ket{M}$ with important contributions from the single transfers $\ket{M^{a}_{b}}$ near the equilibrium geometry at $r=1.75$ bohr. These contributions drop off very quickly, reaching roughly one tenth of their maximal value by $r=3.3$ bohr. The single swaps $\ket{M^{a}_{a}}$ are small, and on either side of zero. Individual coefficients for the swaps are plotted to convey the message. Contributions from the double swaps $\ket{M^{ab}_{ab}}$ are maximal at $r=3.0$ bohr, again just before the critical points $r_{C}$.

The coupling of the reference $\ket{M}$ to each state goes to zero as $r$ becomes large, but in different ways. One might expect that all the TDM elements go to zero, but this is not the case. Couplings between the reference and swaps have non-zero $\gamma_{k}$ and $P_{kl}$ elements within a VBS, but with opposite sign. In the large $r$ limit, the intra-VBS 1-body integrals $h_{a_{1}a_{1}}=h_{a_{2}a_{2}}$ are degenerate, as are the 2-body direct integrals

while the exchange integrals are always the same by symmetry

so the sum of $\gamma_{k}$ and $P_{kl}$ contributions will vanish. There are non-zero $D_{kl}$ elements between VBS, which at dissociation become for the particular VBS $a$ and all other VBS $b$

The exchange integrals $V_{ijji}$ between VBS are zero as the orbitals are localized, and the direct elements are degenerate

so the couplings $\braket{M}{\hat{H}_{C}}{M^{a}_{a}}=0$ when $r$ becomes large.

Couplings between the reference and single transfers have no contribution from one-body elements past $r_{C}$. The TDM elements $\gamma_{k}$ themselves vanish as $J$ develops a second small singular value. While small, it is orders of magnitude larger than the fundamental singular value $\sigma_{N}$. At large $r$ there are small, but non-zero, $D_{kl}$ elements between VBS, i.e. for the coupling of the reference with the single transfer state $\ket{M^{b}_{a}}$ describing the pair transfer $a\rightarrow b$

The direct integrals between VBS are degenerate

while the exchange integrals between VBS are zero, so the $D_{kl}$ couplings between the reference and the single transfer give no contribution. There are non-zero $P_{kl}$ elements between VBS, in particular there are the expected forward scattering

in addition to much smaller backward scattering

elements. However, as these will be weighted by exchange (real pair-transfer) integrals

when $r$ is large these give no contribution either and hence $\braket{M}{\hat{H}_{C}}{M^{b}_{a}}=0$.

For couplings of the reference $\ket{M}$ with doubles, the matrix $J$ always has a second small singular value, though again it is generally much larger than the fundamental $\sigma_{N}$. As a result, none of the doubles couple to the reference through 1-body elements $\gamma_{k}$. At large $r$, the coupling between the reference $\ket{M}$ and a double swap $\ket{M^{ab}_{ab}}$ has non-zero elements

which give zero contribution as the corresponding integrals are again symmetric. There are also very small

elements between VBS, which also give zero contribution as the exchange integrals between VBS go to zero for large $r$. Couplings of the reference $\ket{M}$ with pair plus transfers $\ket{M^{ac}_{ab}}$ and double transfers $\ket{M^{cd}_{ab}}$ only have non-zero contributions from $P_{kl}$ type elements, though the largest of these is on the order of $10^{-6}$. Thus the reference does not couple at all with doubles at large $r$.

For H₈, there are a limited set of triple excitations, as well as a quadruple swap. The triples always have at least 3 small singular values, with an additional small singular value developed for those involving a transfer. The quadruple swap always has 4 small singular values. As a result, the first and second cofactors of $J$ are always near-zero and there is no coupling to the reference at any value of $r$.

Following these observations, analogues of the Slater-Condon rules for the RG reference $\ket{M}$ may be stated:

1. 1.

   For a $k$-fold excitation, $J$ always has $k$ small singular values. Additional small singular values appear for each transfer once $\Delta\varepsilon_{a}<|g|$: a single transfer acquires one extra small singular value, double transfers acquire two extra small singular values, etc.

2. 2.

   If $J$ has one small singular value (in particular $\sigma_{N}$), there are one- and two-body couplings.

3. 3.

   If $J$ has two small singular values, only two-body couplings are present.

4. 4.

   If $J$ has more than two small singular values, there is no coupling for a two-body operator.

Having established that a given RG reference $\ket{M}$ has no couplings beyond doubles, and thus that RGCISD $\approx$ DOCI, it would be nice to reduce the calculation to something actually feasible. Building the RGCISD matrix means computing the complete TDM ($\mathcal{O}(N^{4})$) for the $M(N-M)$ single excitations, and $\binom{M}{2}\binom{N-M}{2}$ double excitations, which scales like $\mathcal{O}(N^{4}M^{8})$. Even building the RGCIS matrix has a cost on the order of $\mathcal{O}(N^{4}M^{4})$, which is more expensive than the resulting matrix diagonalization. As the RGCISD wavefunction is dominated by the reference $\ket{M}$, single reference perturbation theory (PT) would seem to be an attractive alternative.

The standard Rayleigh-Schrödinger (RS) PT construction decomposes the Hamiltonian we wish to solve, in this case $\hat{H}_{C}$ the Coulomb Hamiltonian (12), as an exactly solvable reference $\hat{H}_{0}$ plus a perturbation

and builds energetic corrections order by order with the eigenvectors of $\hat{H}_{0}$

From a given reference $\ket{\psi^{(0)}_{0}}$, the 2nd order correction to the energy is

In the present case, the chosen reference would be $\hat{H}_{0}=\hat{H}_{VB}$ (45) which is far from $\hat{H}_{C}$. Worse, as $\ket{\psi^{(0)}_{0}}=\ket{M}$ is not the ground state of (45), the 2nd order correction (61) can be positive. For H₈ this does happen, and is so disastrous the curve won’t be presented. Such failure was also observed for RSPT corrections to RG states for individual electrons.^{93, 94}

A much better treatment is obtained by defining $\hat{H}_{0}$ in terms of the expected values of the target $\hat{H}_{C}$

which is the Epstein⁹⁵-Nesbet⁹⁶ (EN) partitioning. The state $\ket{M}$ is the now the ground state of $\hat{H}_{0}$, and the 2nd order energy correction is

As only TDM elements involving the reference are required, this correction is computable with $\mathcal{O}(N^{4}M^{2})$ cost. ENPT is known to have issues with size-consistency. In the present case, the results appear to be size-consistent so long as the orbitals are localized. If the orbitals are not localized, then RG states are not size-consitent either.

ENPT2 corrections employing single and double excitations for H₈ are shown in Figure (4). Raw energies are not plotted as it is impossible to discern the curves. Notice that the difference between ENPT2 using singles is very close to RGCISD.

The same is generally true for the Stair-Evangelista isomers of H₁₀⁹⁷: the chain of 10 equidistant H atoms, the circular ring of 10 equidistant H atoms, the 3-4-3 planar sheet of 10 equidistant H atoms and the pyramid of 10 equidistant H atoms. The pyramid is in fact a tetrahedron with an H atom at each vertex and an H atom at the midpoint of each edge. OO-DOCI is a reasonable description for the chain and the ring, but not for the sheet nor the pyramid. In all cases, RGCISD $\approx$ DOCI, with an agreement on the order of $10^{-9}$ for the chain and the ring, and $10^{-6}$ for the sheet and the pyramid.⁶³ Figure 5 shows the errors of ENPT2 with respect to RGCISD. In all cases the disagreement between the two is less than 1 mE_h. Note that for the first three points of the sheet, and the first four points of the pyramid, there is a curve crossing and the refence RG state is different.⁶³ Even so, the agreement between ENPT2 with singles and RGCISD is quite good for a fraction of the cost.

Finally, a variational RG treatment and ENPT2 with singles of linear equidistant H₅₀ is shown in figure 6. This system is too large to treat with DOCI so the main point is that ENPT2 can be computed in a reasonable time. By computing the TDM elements in parallel, this computation requires roughly an hour on a modest desktop machine. As expected, the reference RG state is $\ket{M}$, and the ENPT2 correction is largest near the equilibrium geometry before decaying to zero at dissociation. As OO-DOCI is unfeasible, the orbitals employed are those optimal for the Piris natural orbital functional PNOF7^{98, 99, 100} from ref.¹⁰¹ Natural orbital functionals are close cousins of geminal wavefunctions: the 2-RDM elements are explicit functions of the 1-RDM elements, which is the case for the antisymmetrized product of strongly orthogonal geminals (APSG)^{102, 103, 104} and the antisymmetrized geminal power (AGP).^{105, 106, 107} PNOF7 is not N-representable, but can be loosely understood as having intrapair elements like APSG and interpair elements like AGP.³³