# Pertubative corrections for Hartree-Fock-like Algebraic Bethe Ansatz
analogue.

Jean-David Moisset  Laurie Carrier  Paul A. Johnson
<EMAIL_ADDRESS>Département de chimie, Université Laval, Québec,
Québec, G1V 0A6, Canada

###### Abstract

Rayleigh-Schrödinger perturbation theory corrections are developed for an
algebraic Bethe ansatz of individual electrons. Numerical results are
ambiguous and would need either an orbital optimization or a configuration
interaction singles solution to be satisfactory. Two different expressions are
computed to be certain that results are correct.

##  I Introduction

Strongly-correlated systems present a challenge in quantum chemistry. The
dominant behaviour is not independent electrons and thus methods built in that
framework do not perform well. Powerful methods exist but are generally either
too difficult or too expensive, though substantial improvement is being
made.chan:2011 ; booth:2013 ; holmes:2017 ; li:2018  A promising alternative
is to consider wavefunctions built from weakly-correlated pairs of electrons,
or geminals. This idea is not new,hurley:1953 ; coleman:1965 ; silver:1969 ;
silver:1970 ; coleman:1997  though a strong renewed interestsurjan:1999 ;
kobayashi:2010 ; surjan:2012 ; neuscamman:2012 ; peter:2013 ; johnson:2013 ;
stein:2014 ; boguslawski:2014a ; boguslawski:2014b ; boguslawski:2014c ;
tecmer:2014 ; henderson:2014a ; henderson:2014b ; shepherd:2014 ; bulik:2015 ;
pastorczak:2015 ; henderson:2019 ; khamoshi:2019 ; henderson:2020 ; dutta:2020
; harsha:2020 ; khamoshi:2020 ; khamoshi:2021 ; dutta:2021 ; marie:2021  has
shown that these types of wavefunctions correctly describe many bond-breaking
processes. In particular, the antisymmetric product of 1-reference orbital
geminals (AP1roG), or equivalently pair coupled-cluster doubles (pCCD), scales
like $\mathcal{O}(N^{4})$ and dissociates Hydrogen chains almost perfectly.
Their principal drawback is that they must be solved projectively rather than
variationally to be feasible.

We are building approaches for strongly-correlated electrons from exactly
solvable models. In particular we are using the algebraic Bethe ansatz
(ABA)bethe ; faddeev:1981  solution to the reduced Bardeen-Cooper-Schrieffer
Hamiltonian,bardeen:1957a ; bardeen:1957b  which we refer to as Richardson-
Gaudin (RG)richardson:1963 ; richardson:1964 ; richardson:1965 ; gaudin:1976 ;
dukelsky:2004  states, as a mean-field geminal wavefunction.johnson:2020 ;
fecteau:2020 ; fecteau:2021 ; johnson:2021  The ABA is an approach capable of
solving a large class of models both in quantum mechanics,QISM_book  and in
2-dimensional classical statistical mechanics.baxter_sm  We are studying RG
states as the general mean-field for pairs of electrons, the so-called
antisymmetric product of interacting geminals (APIG), is intractable to
compute whereas RG states have polynomial cost and may be improved upon
systematically. In a previous contributioncarrier:2020  we studied an ABA for
individual electrons to gauge how well the general mean-field, in that case
Hartree-Fock (HF), was reproduced. We demonstrated numerically that there is
essentially no approximation at the mean-field level.

The purpose of this contribution is to develop perturbative corrections for
the ABA for electrons. The final equations look quite similar to the Møller-
Plesset (MP) corrections to Hartree-Fock, though the singles give a non-zero
contribution.MEST  Finally, the numerical results are ambiguous so we will not
further explore the ABA for individual electrons. To be certain, we derived
two different expressions that gave consistent numerical results. More
optimization parameters would be required, either in the form of an orbital
optimization or a configuration interaction (CI) singles solution.

In section II we briefly review the ABA for individual electrons in a general
spin-orbital basis before moving to an unrestricted basis and calculating
reduced density matrix (RDM) and transition density matrix (TDM) elements.
Section III presents second-order Rayleigh-Schrödinger perturbation theory
(RSPT) corrections computed in the original “primitive” basis, while section
IV presents RSPT expressions in the basis of ABA objects (the analogue of
molecular orbitals).

##  II ABA for individual electrons

The ABA provides the eigenvectors of the Hamiltonian

$\displaystyle\hat{H}_{ABA}$
$\displaystyle=\sum_{i}\varepsilon_{i}a^{\dagger}_{i}a_{i}+g\sum_{ij}a^{\dagger}_{i}a_{j}$
(1)

which describes an aufbau filling of the lowest single particle states
$\\{\varepsilon\\}$ along with an isotropic scattering $g$ of electrons
between each site. Obviously, this Hamiltonian may be solved exactly by
diagonalization. The point here is to try to understand the ABA in a simple
case so that we can use that information for RG, where the analogue of HF,
APIG, is intractable.

We employ spin-orbitals that have the usual structure,

$\displaystyle[a^{\dagger}_{i},a_{j}]_{+}=\delta_{ij}$ (2)

for which $i$ and $j$ are a complete set of spin and orbital indices. The ABA
is built from a set of objects in terms of a complex number $u$:

$\displaystyle a^{\dagger}(u)$
$\displaystyle=\sum_{i}\frac{a^{\dagger}_{i}}{u-\varepsilon_{i}}$ (3)
$\displaystyle a(u)$ $\displaystyle=\sum_{i}\frac{a_{i}}{u-\varepsilon_{i}}.$
(4)

These objects have the structure

$\displaystyle[a^{\dagger}(u),a(v)]_{+}=\frac{A(u)-A(v)}{u-v}$ (5)

where the scalar function $A(u)$ is

$\displaystyle A(u)=\frac{1}{g}-\sum_{i}\frac{1}{u-\varepsilon_{i}}.$ (6)

In the limit that $u$ and $v$ are the same, this is well defined, with the
result

$\displaystyle[a^{\dagger}(u),a(v)]_{+}=\frac{\partial A(u)}{\partial u}.$ (7)

###  II.1 ABA solution

The point of the ABA is that eigenvectors of (1) are products of
$a^{\dagger}(u)$ whose arguments satisfy a set of equations. Usually these
equations are coupled, but for this particular case they are not. We denote
the ABA states based on their _rapidities_ $\\{u\\}$ as

$\displaystyle\ket{\\{u\\}}=a^{\dagger}(u_{1})a^{\dagger}(u_{2})\dots
a^{\dagger}(u_{2M})\ket{\theta}$ (8)

where the vacuum state $\ket{\theta}$ is understood such that
$a_{j}\ket{\theta}=a(u)\ket{\theta}=0$. We emphasize that there are $2M$
electrons in $2N$ spin-orbitals. The state $\\{u\\}$ is a Slater determinant
of $2M$ electrons in the basis of $a^{\dagger}(u)$. To establish the action of
(1) on this state, the strategy is to move $\hat{H}_{ABA}$ to the right until
it destroys the vacuum. First, it is easy to confirm that

$\displaystyle[\hat{H}_{ABA},a^{\dagger}(u)]=ua^{\dagger}(u)-gA(u)\sum_{i}a^{\dagger}_{i}$
(9)

and after a little work

$\displaystyle\hat{H}_{ABA}\ket{\\{u\\}}=\sum_{a\in
occ}u_{a}\ket{\\{u\\}}-g\sum_{a\in
occ}A(u_{a})\sum_{i}a^{\dagger}_{i}\ket{\\{u\\}_{a}}.$ (10)

Thus we have an eigenvector provided that for each $a$

$\displaystyle
A(u_{a})=\frac{1}{g}-\sum_{i}\frac{1}{u_{a}-\varepsilon_{i}}=0.$ (11)

These are called the Bethe ansatz equations (BAE). Usually they are coupled,
and must be solved for each ABA state. In this particular case, the BAE are
decoupled, and hence one solution suffices _for all_ the ABA eigenvectors. The
BAE (11) have $2N$ solutions, but each eigenvector only contains $2M$. To be
in line with HF, we will refer to the rapidities that are present in the
ground state as _occupied_ and those that do not as _virtual_. The summations
over $a$ in equation (10) are of course over only the occupied rapidities. In
this particular case, the BAE may be solved by diagonalizing the matrix

$\displaystyle\begin{pmatrix}\varepsilon_{1}+g&g&\dots&g\\\
g&\varepsilon_{2}+g&\dots&g\\\ \vdots&\vdots&\ddots&\vdots\\\
g&g&\dots&\varepsilon_{2N}+g\end{pmatrix}\rightarrow\begin{pmatrix}u_{1}&0&\dots&0\\\
0&u_{2}&\dots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\dots&u_{2N}\end{pmatrix}$ (12)

which is substantially more stable numerically.

###  II.2 Scalar Products

Scalar products are computed by the form factor approach.zhou:2002 ;
faribault:2008 ; faribault:2010 ; GB:2011  Namely to evaluate
$\braket{\\{v\\}}{a^{\dagger}_{i}a_{j}}{\\{u\\}}$ move $a_{j}$ to the right
until it destroys the vacuum. This yields a sum over terms that are evaluated
as specific limits of a general scalar product formula. For arbitrary
$\\{u\\}$ and $\\{v\\}$ we can evaluate the scalar product easily as the
determinant of a $2M\times 2M$ matrix

$\displaystyle\braket{\\{v\\}}{\\{u\\}}=\det_{ab\in
occ}\left(\frac{A(u_{a})-A(v_{b})}{u_{a}-v_{b}}\right).$ (13)

Henceforth, we take one of the sets, say $\\{v\\}$ to be on-shell, i.e.
$A(v_{b})=0$ for each $b$ giving

$\displaystyle\braket{\\{v\\}}{\\{u\\}}=\prod_{c\in occ}A(v_{c})\det_{ab\in
occ}\left(\frac{1}{u_{a}-v_{b}}\right).$ (14)

In principal this expression could be simplified as the determinant of a
Cauchy matrix has a simple closed form expression, but this is not productive
as we will see briefly. All the required scalar products are specific limits
of (14). In particular, the norm is the limit when
$\\{u\\}\rightarrow\\{v\\}$:

$\displaystyle\braket{\\{v\\}}{\\{v\\}}=\prod_{c\in occ}\frac{\partial
A(v_{c})}{\partial v_{c}}$ (15)

since the only surviving terms come from indeterminate forms on the diagonal
of the determinant.

To evaluate the 1-electron reduced density matrix (1-RDM) we will need the
structure constants

$\displaystyle[a_{j},a^{\dagger}(u)]_{+}$
$\displaystyle=\frac{1}{u-\varepsilon_{j}}$ (16)
$\displaystyle[a^{\dagger}_{i},a^{\dagger}(u)]_{+}$ $\displaystyle=0$ (17)

to move $a_{j}$ to the right. The result is

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i}a_{j}}{\\{u\\}}$
$\displaystyle=\sum_{a\in
occ}(-1)^{a-1}[a_{j},a^{\dagger}(u_{a})]_{+}\braket{\\{v\\}}{a^{\dagger}_{i}}{\\{u\\}_{a}}$
(18) $\displaystyle=\sum_{a\in
occ}\frac{\braket{\\{v\\}}{\\{u\\}_{a\rightarrow i}}}{u_{a}-\varepsilon_{i}}$
(19)

where $\\{u\\}_{a}$ denotes the set $\\{u\\}$ without $u_{a}$, and the state
$\ket{\\{u\\}_{a\rightarrow i}}$ is $\ket{\\{u\\}}$ with $a^{\dagger}(u_{a})$
replaced with $a^{\dagger}_{i}$ _in the correct position_ which removes the
sign. Now, the original operators are the residues of the ABA operators

$\displaystyle
a^{\dagger}_{i}=\lim_{u\rightarrow\varepsilon_{i}}(u-\varepsilon_{i})a^{\dagger}(u)$
(20)

which carries through to the scalar product (14)

$\displaystyle\braket{\\{v\\}}{\\{u\\}_{a\rightarrow\varepsilon_{i}}}$
$\displaystyle=\lim_{u_{a}\rightarrow
i}(u_{a}-\varepsilon_{i})\braket{\\{v\\}}{\\{u\\}}$ (21) $\displaystyle=-\det
G^{a}_{i}\prod_{c(\neq a)\in occ}A(u_{c}).$ (22)

Here the matrix $G^{a}_{i}$ is the matrix from (14) with the $a$th column
replaced:

$\displaystyle\begin{pmatrix}\frac{1}{u_{a}-v_{1}}\\\ \frac{1}{u_{a}-v_{2}}\\\
\vdots\\\
\frac{1}{u_{a}-v_{M}}\end{pmatrix}\rightarrow\begin{pmatrix}\frac{1}{\varepsilon_{i}-v_{1}}\\\
\frac{1}{\varepsilon_{i}-v_{2}}\\\ \vdots\\\
\frac{1}{\varepsilon_{i}-v_{M}}\end{pmatrix}.$ (23)

Taking the limit $\\{u\\}\rightarrow\\{v\\}$ and normalizing we arrive at the
1-RDM elements

$\displaystyle\gamma_{ij}=\frac{\braket{\\{v\\}}{a^{\dagger}_{i}a_{j}}{\\{v\\}}}{\braket{\\{v\\}}{\\{v\\}}}=\sum_{a\in
occ}\frac{1}{(v_{a}-\varepsilon_{i})(v_{a}-\varepsilon_{j})}\frac{1}{\frac{\partial
A(v_{a})}{\partial v_{a}}}.$ (24)

The ABA operators represent electrons in orthogonal but not normal orbitals.
The quantity

$\displaystyle\frac{\partial A(v_{a})}{\partial
v_{a}}=\sum_{j}\frac{1}{(v_{a}-\varepsilon_{j})^{2}}$ (25)

plays the role of the (diagonal) overlap matrix.

Similarly for two electron operators,

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}}{\\{u\\}}=\sum_{a<b\in
occ}\begin{vmatrix}\frac{1}{u_{a}-\varepsilon_{k}}&\frac{1}{u_{a}-\varepsilon_{l}}\\\
\frac{1}{u_{b}-\varepsilon_{k}}&\frac{1}{u_{b}-\varepsilon_{l}}\end{vmatrix}\braket{\\{v\\}}{\\{u\\}_{ab\rightarrow
ij}}.$ (26)

and

$\displaystyle\braket{\\{v\\}}{\\{u\\}_{ab\rightarrow ij}}=\det
G^{ab}_{ij}\prod_{c(\neq a,b)\in occ}A(u_{c})$ (27)

leads to the 2-electron reduced density matrix (2-RDM) elements

$\displaystyle\Gamma_{ijkl}$
$\displaystyle=\frac{\braket{\\{v\\}}{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}}{\\{v\\}}}{\braket{\\{v\\}}{\\{v\\}}}$
(28) $\displaystyle=\sum_{a<b\in
occ}\begin{vmatrix}\frac{1}{v_{a}-\varepsilon_{k}}&\frac{1}{v_{a}-\varepsilon_{l}}\\\
\frac{1}{v_{b}-\varepsilon_{k}}&\frac{1}{v_{b}-\varepsilon_{l}}\end{vmatrix}\begin{vmatrix}\frac{1}{v_{a}-\varepsilon_{i}}&\frac{1}{v_{a}-\varepsilon_{j}}\\\
\frac{1}{v_{b}-\varepsilon_{i}}&\frac{1}{v_{b}-\varepsilon_{j}}\end{vmatrix}\frac{1}{\frac{\partial
A(v_{a})}{\partial v_{a}}\frac{\partial A(v_{b})}{\partial v_{b}}}$ (29)
$\displaystyle=\gamma_{ik}\gamma_{jl}-\gamma_{il}\gamma_{jk}.$ (30)

The 2-RDM factors into 1-RDM information as the wavefunction is factored into
1-electron wavefunctions. This extends to any order

$\displaystyle\Gamma^{(N)}_{i_{1}\dots i_{N}j_{1}\dots j_{N}}$
$\displaystyle=\frac{\braket{\\{v\\}}{a^{\dagger}_{i_{1}}\dots
a^{\dagger}_{i_{N}}a_{j_{N}}\dots
a_{j_{1}}}{\\{v\\}}}{\braket{\\{v\\}}{\\{v\\}}}$ (31)
$\displaystyle=\sum_{a_{1}<\dots a_{N}\in
occ}\begin{vmatrix}\frac{1}{v_{a_{1}}-\varepsilon_{j_{1}}}&\dots&\frac{1}{v_{a_{1}}-\varepsilon_{j_{N}}}\\\
\vdots&\ddots&\vdots\\\
\frac{1}{v_{a_{N}}-\varepsilon_{j_{1}}}&\dots&\frac{1}{v_{a_{N}}-\varepsilon_{j_{N}}}\end{vmatrix}\begin{vmatrix}\frac{1}{v_{a_{1}}-\varepsilon_{i_{1}}}&\dots&\frac{1}{v_{a_{1}}-\varepsilon_{i_{N}}}\\\
\vdots&\ddots&\vdots\\\
\frac{1}{v_{a_{N}}-\varepsilon_{i_{1}}}&\dots&\frac{1}{v_{a_{N}}-\varepsilon_{i_{N}}}\end{vmatrix}\frac{1}{\frac{\partial
A(v_{a_{1}})}{\partial v_{a_{1}}}\dots\frac{\partial A(v_{a_{N}})}{\partial
v_{a_{N}}}}$ (32)
$\displaystyle=\begin{vmatrix}\gamma_{i_{1}j_{1}}&\dots&\gamma_{i_{1}j_{N}}\\\
\vdots&\ddots&\vdots\\\
\gamma_{i_{N}j_{1}}&\dots&\gamma_{i_{N}j_{N}}\end{vmatrix}.$ (33)

Transition density matrix (TDM) elements are computed with the same approach
and the results are even simpler. For single excitations, i.e. replacing one
$v_{a}$ with $v_{b}$, denoted set-wise as $\\{v\\}^{b}_{a}$, there is only one
non-zero contribution. All others are proportional to $A(v_{b})$ without a
vanishing denominator and hence are identically zero. The non-zero
contributions from single excitations are

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i}a_{j}}{\\{v\\}^{b}_{a}}=\frac{1}{(v_{b}-\varepsilon_{j})(v_{a}-\varepsilon_{i})}\prod_{c(\neq
a)\in occ}\frac{\partial A(v_{c})}{\partial v_{c}}$ (34)

and

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}}{\\{v\\}^{b}_{a}}=\sum_{c(\neq
a)\in
occ}\begin{vmatrix}\frac{1}{v_{b}-\varepsilon_{k}}&\frac{1}{v_{b}-\varepsilon_{l}}\\\
\frac{1}{v_{c}-\varepsilon_{k}}&\frac{1}{v_{c}-\varepsilon_{l}}\end{vmatrix}\begin{vmatrix}\frac{1}{v_{a}-\varepsilon_{k}}&\frac{1}{v_{a}-\varepsilon_{l}}\\\
\frac{1}{v_{c}-\varepsilon_{k}}&\frac{1}{v_{c}-\varepsilon_{l}}\end{vmatrix}\prod_{d(\neq
a,c)\in occ}\frac{\partial A(v_{d})}{\partial v_{d}}.$ (35)

Double excitations have the only non-zero contributions

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i}a^{\dagger}_{j}a_{l}a_{k}}{\\{v\\}^{cd}_{ab}}=\begin{vmatrix}\frac{1}{v_{c}-\varepsilon_{k}}&\frac{1}{v_{c}-\varepsilon_{l}}\\\
\frac{1}{v_{d}-\varepsilon_{k}}&\frac{1}{v_{d}-\varepsilon_{l}}\end{vmatrix}\begin{vmatrix}\frac{1}{v_{a}-\varepsilon_{i}}&\frac{1}{v_{a}-\varepsilon_{j}}\\\
\frac{1}{v_{b}-\varepsilon_{i}}&\frac{1}{v_{b}-\varepsilon_{j}}\end{vmatrix}\prod_{c(\neq
a,b)\in occ}\frac{\partial A(v_{c})}{\partial v_{c}}.$ (36)

###  II.3 Unrestricted Orbitals

To arrive at restricted results, we work out the result in unrestricted
orbitals before taking the limit where the $\alpha$ and $\beta$ elements are
identical. The elements can be separated into two sets. It’s equivalent to
having two separate ABA set-ups:

$\displaystyle\hat{H}^{\alpha}_{ABA}$
$\displaystyle=\sum_{i}\varepsilon^{\alpha}_{i}a^{\dagger}_{i}a_{i}+g^{\alpha}\sum_{ij}a^{\dagger}_{i}a_{j}$
(37) $\displaystyle\hat{H}^{\beta}_{ABA}$
$\displaystyle=\sum_{i}\varepsilon^{\beta}_{i}a^{\dagger}_{i}a_{i}+g^{\beta}\sum_{ij}a^{\dagger}_{i}a_{j}$
(38)

with solutions the roots of the equations:

$\displaystyle
A(v^{\alpha}):=\frac{1}{g^{\alpha}}-\sum_{i}\frac{1}{v^{\alpha}-\varepsilon^{\alpha}_{i}}$
$\displaystyle=0$ (39) $\displaystyle
A(v^{\beta}):=\frac{1}{g^{\beta}}-\sum_{i}\frac{1}{v^{\beta}-\varepsilon^{\beta}_{i}}$
$\displaystyle=0.$ (40)

The scalar products simplify as the matrix entering in (14) becomes block
diagonal. Specifically, it being understood that
$\\{v\\}=\\{v^{\alpha}\\}\cup\\{v^{\beta}\\}$, we have

$\displaystyle\left(\hat{H}^{\alpha}_{ABA}+\hat{H}^{\beta}_{ABA}\right)\ket{\\{v\\}}=\left(\sum_{a\in
occ}v^{\alpha}_{a}+v^{\beta}_{a}\right)\ket{\\{v\\}}$ (41)

and the scalar product with an arbitrary ABA vector $\ket{\\{u\\}}$ is

$\displaystyle\braket{\\{v\\}}{\\{u\\}}\prod_{c\in
occ}A(u^{\alpha}_{c})A(u^{\beta}_{c})\det_{ab\in
occ}\left(\frac{1}{u^{\alpha}_{a}-v^{\alpha}_{b}}\right)\det_{a^{\prime}b^{\prime}\in
occ}\left(\frac{1}{u^{\beta}_{a^{\prime}}-v^{\beta}_{b^{\prime}}}\right).$
(42)

The norm of the on-shell state becomes

$\displaystyle\braket{\\{v\\}}{\\{v\\}}=\prod_{a\in occ}\frac{\partial
A(v^{\alpha}_{a})}{\partial v^{\alpha}_{a}}\frac{\partial
A(v^{\beta}_{a})}{\partial v^{\beta}_{a}}.$ (43)

The 1-RDM elements are:

$\displaystyle\gamma^{\sigma}_{ij}=\frac{\braket{\\{v\\}}{a^{\dagger}_{i\sigma}a_{j\sigma}}{\\{v\\}}}{\braket{\\{v\\}}{\\{v\\}}}=\sum_{a\in
occ}\frac{1}{(v^{\sigma}_{a}-\varepsilon^{\sigma}_{i})(v^{\sigma}_{a}-\varepsilon^{\sigma}_{j})}\frac{1}{\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}}$ (44)

and the 2-RDM elements reduce to 1-RDM elements

$\displaystyle\Gamma^{\sigma\tau}_{ijkl}=\frac{\braket{\\{v\\}}{a^{\dagger}_{i\sigma}a^{\dagger}_{j\tau}a_{l\tau}a_{k\tau}}{\\{v\\}}}{\braket{\\{v\\}}{\\{v\\}}}=\gamma^{\sigma}_{ik}\gamma^{\tau}_{jl}-\delta_{\sigma\tau}\gamma^{\sigma}_{il}\gamma^{\sigma}_{jk}.$
(45)

As in the previous section, excited states are labelled by which occupieds are
removed, and which virtuals are added, i.e. the state
$\ket{\\{v\\}^{p\sigma}_{a\sigma}}$ is the state with $a\sigma$ occupied
orbital replaced by the $p\sigma$ virtual orbital.

We will require transition density matrix elements, for which $\sigma$ denotes
either spin, while $\bar{\sigma}$ denotes the opposite spin of $\sigma$. They
are obtained in the same manner as for the spin-orbital basis:

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i\sigma}a_{j\sigma}}{\\{v\\}^{p\sigma}_{a\sigma}}=\frac{\braket{\\{v\\}}{\\{v\\}}}{(v^{\sigma}_{p}-\varepsilon^{\sigma}_{j})(v^{\sigma}_{a}-\varepsilon^{\sigma}_{i})\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}}$ (46)

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i\sigma}a^{\dagger}_{j\sigma}a_{l\sigma}a_{k\sigma}}{\\{v\\}^{p\sigma}_{a\sigma}}=\sum_{c(\neq
a)\in
occ}\begin{vmatrix}\frac{1}{v^{\sigma}_{p}-\varepsilon^{\sigma}_{k}}&\frac{1}{v^{\sigma}_{p}-\varepsilon^{\sigma}_{l}}\\\
\frac{1}{v^{\sigma}_{c}-\varepsilon^{\sigma}_{k}}&\frac{1}{v^{\sigma}_{c}-\varepsilon^{\sigma}_{l}}\end{vmatrix}\begin{vmatrix}\frac{1}{v^{\sigma}_{a}-\varepsilon^{\sigma}_{i}}&\frac{1}{v^{\sigma}_{a}-\varepsilon^{\sigma}_{j}}\\\
\frac{1}{v^{\sigma}_{c}-\varepsilon^{\sigma}_{i}}&\frac{1}{v^{\sigma}_{c}-\varepsilon^{\sigma}_{j}}\end{vmatrix}\frac{\braket{\\{v\\}}{\\{v\\}}}{\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}\frac{\partial
A(v^{\sigma}_{c})}{\partial v^{\sigma}_{c}}}.$ (47)

Strictly speaking, the summation excludes the term $c=a$, but this term would
give zero contribution so will be included to simplify later summations. The
last TDM element for single excitations is

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i\sigma}a^{\dagger}_{j\bar{\sigma}}a_{l\bar{\sigma}}a_{k\sigma}}{\\{v\\}^{p\sigma}_{a\sigma}}=\frac{\braket{\\{v\\}}{\\{v\\}}}{(v^{\sigma}_{p}-\varepsilon^{\sigma}_{k})(v^{\sigma}_{a}-\varepsilon^{\sigma}_{i})\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}}\gamma^{\bar{\sigma}}_{jl}.$ (48)

Double excitations only couple with the ground state through two-electron
operators:

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i\sigma}a^{\dagger}_{j\sigma}a_{l\sigma}a_{k\sigma}}{\\{v\\}^{p\sigma
q\sigma}_{a\sigma
b\sigma}}=\begin{vmatrix}\frac{1}{v^{\sigma}_{p}-\varepsilon^{\sigma}_{k}}&\frac{1}{v^{\sigma}_{p}-\varepsilon^{\sigma}_{l}}\\\
\frac{1}{v^{\sigma}_{q}-\varepsilon^{\sigma}_{k}}&\frac{1}{v^{\sigma}_{q}-\varepsilon^{\sigma}_{l}}\end{vmatrix}\begin{vmatrix}\frac{1}{v^{\sigma}_{a}-\varepsilon^{\sigma}_{i}}&\frac{1}{v^{\sigma}_{a}-\varepsilon^{\sigma}_{j}}\\\
\frac{1}{v^{\sigma}_{b}-\varepsilon^{\sigma}_{i}}&\frac{1}{v^{\sigma}_{b}-\varepsilon^{\sigma}_{j}}\end{vmatrix}\frac{\braket{\\{v\\}}{\\{v\\}}}{\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}\frac{\partial
A(v^{\sigma}_{b})}{\partial v^{\sigma}_{b}}}$ (49)

and

$\displaystyle\braket{\\{v\\}}{a^{\dagger}_{i\sigma}a^{\dagger}_{j\bar{\sigma}}a_{l\bar{\sigma}}a_{k\sigma}}{\\{v\\}^{p\sigma
q\bar{\sigma}}_{a\sigma
b\bar{\sigma}}}=\frac{\braket{\\{v\\}}{\\{v\\}}}{(v^{\sigma}_{a}-\varepsilon^{\sigma}_{i})(v^{\bar{\sigma}}_{b}-\varepsilon^{\bar{\sigma}}_{j})(v^{\sigma}_{p}-\varepsilon^{\sigma}_{k})(v^{\bar{\sigma}}_{q}-\varepsilon^{\bar{\sigma}}_{l})\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}\frac{\partial
A(v^{\bar{\sigma}}_{b})}{\partial v^{\bar{\sigma}}_{b}}}.$ (50)

All other TDM elements either vanish or will not connect through the
Hamiltonian.

##  III PT expressions: Primitive basis

We wish to solve the Coulomb Hamiltonian for molecules

$\displaystyle\hat{H}_{C}=\sum_{ij}h_{ij}\sum_{\sigma}a^{\dagger}_{i\sigma}a_{\sigma}+\frac{1}{2}\sum_{ijkl}V_{ijkl}\sum_{\sigma\tau}a^{\dagger}_{i\sigma}a^{\dagger}_{j\tau}a_{l\tau}a_{k\sigma}$
(51)

where the 1- and 2-electron integrals are expressed in a basis $\\{\phi\\}$

$\displaystyle h_{ij}$ $\displaystyle=\int
d\mathbf{r}\phi^{*}_{i}(\mathbf{r})\left(-\frac{1}{2}\nabla^{2}-\sum_{I}\frac{Z_{I}}{|\mathbf{r}-\mathbf{R}_{I}|}\right)\phi_{j}(\mathbf{r})$
(52) $\displaystyle V_{ijkl}$ $\displaystyle=\int
d\mathbf{r}_{1}d\mathbf{r}_{2}\frac{\phi^{*}_{i}(\mathbf{r}_{1})\phi^{*}_{j}(\mathbf{r}_{2})\phi_{k}(\mathbf{r}_{1})\phi_{l}(\mathbf{r}_{2})}{|\mathbf{r}_{1}-\mathbf{r}_{2}|}.$
(53)

Using the ABA for electrons as a wavefunction ansatz, and expanding

$\displaystyle\hat{H}_{C}\ket{\Psi_{k}}=E_{k}\ket{\Psi_{k}}$ (54)

in Rayleigh-Schrödinger perturbation theory (RSPT), we obtain

$\displaystyle\hat{H}_{C}$ $\displaystyle=\hat{H}_{0}+\lambda\hat{H}_{1}$ (55)
$\displaystyle E_{k}$ $\displaystyle=E^{(0)}_{k}+\lambda
E^{(1)}_{k}+\lambda^{2}E^{(2)}_{k}+\dots$ (56) $\displaystyle\ket{\Psi_{k}}$
$\displaystyle=\ket{\Psi^{(0)}_{k}}+\lambda\ket{\Psi^{(1)}_{k}}+\lambda^{2}\ket{\Psi^{(2)}_{k}}+\dots.$
(57)

Collecting powers of $\lambda$, we highlight the zeroth order problem

$\displaystyle\hat{H}_{0}\ket{\Psi^{(0)}_{k}}=E^{(0)}_{k}\ket{\Psi^{(0)}_{k}}$
(58)

which is identically the problem (41). The RSPT first-order correction to the
wavefunctions are

$\displaystyle\ket{\Psi^{(1)}_{k}}=\sum_{l\neq
k}\frac{\braket{\Psi^{(0)}_{l}}{\hat{H}_{1}}{\Psi^{(0)}_{k}}}{E^{(0)}_{k}-E^{(0)}_{l}}\ket{\Psi^{(0)}_{l}}$
(59)

and since our first order Hamiltonian is

$\displaystyle\hat{H}_{1}=\hat{H}_{C}-\left(\hat{H}^{\alpha}_{ABA}+\hat{H}^{\beta}_{ABA}\right)$
(60)

and for $k\neq l$,

$\displaystyle\braket{\Psi^{(0)}_{l}}{\left(\hat{H}^{\alpha}_{ABA}+\hat{H}^{\beta}_{ABA}\right)}{\Psi^{(0)}_{k}}=0$
(61)

we can instead write the first order correction as

$\displaystyle\ket{\Psi^{(1)}_{k}}=\sum_{l\neq
k}\frac{\braket{\Psi^{(0)}_{l}}{\hat{H}_{C}}{\Psi^{(0)}_{k}}}{E^{(0)}_{k}-E^{(0)}_{l}}\ket{\Psi^{(0)}_{l}}.$
(62)

The 2nd order energy is then

$\displaystyle E^{(2)}_{k}=\sum_{l\neq
k}\frac{|\braket{\Psi^{(0)}_{l}}{\hat{H}_{C}}{\Psi^{(0)}_{k}}|^{2}}{E^{(0)}_{k}-E^{(0)}_{l}}.$
(63)

The ABA excited states are simply ABA states with different “occupied”
rapidities since the BAE (11) do not couple rapidities. Like the case for
Hartree-Fock, the only excited states that couple with the ground state
through the Hamiltonian are single and double excitations.

The parameters defining the ABA in an unrestricted basis are
$\\{\varepsilon^{\alpha}\\},\\{\varepsilon^{\beta}\\},g^{\alpha},g^{\beta}$,
and since the 2-RDM factors into the 1-RDM contributions
$\gamma^{\alpha}_{ij},\gamma^{\beta}_{ij}$, as shown previouslycarrier:2020
the variational ground state energy expression is

$\displaystyle E^{ABA}_{U}$ $\displaystyle=E^{(0)}_{0}+E^{(1)}_{0}$ (64)
$\displaystyle=\min_{\\{\varepsilon^{\alpha}\\},\\{\varepsilon^{\beta}\\},g^{\alpha},g^{\beta}}\sum_{ij}h_{ij}\left(\gamma^{\alpha}_{ij}+\gamma^{\beta}_{ij}\right)$
$\displaystyle+\frac{1}{2}\sum_{ijkl}V_{ijkl}\left(\gamma^{\alpha}_{ik}\gamma^{\alpha}_{jl}-\gamma^{\alpha}_{il}\gamma^{\alpha}_{jk}+\gamma^{\alpha}_{ik}\gamma^{\beta}_{jl}+\gamma^{\beta}_{ij}\gamma^{\alpha}_{kl}+\gamma^{\beta}_{ik}\gamma^{\beta}_{jl}-\gamma^{\beta}_{il}\gamma^{\beta}_{jk}\right).$
(65)

To arrive at a restricted expression, we take the $\alpha$ and $\beta$
parameters to be identical. This amounts to solving for one ABA Hamiltonian,
and using states that use each rapidity twice (once for each spin projection).
The resulting energy expression is

$\displaystyle E^{ABA}_{R}$
$\displaystyle=\min_{\\{\varepsilon\\},g}2\sum_{ij}h_{ij}\gamma_{ij}+\sum_{ijkl}V_{ijkl}(2\gamma_{ik}\gamma_{jl}-\gamma_{il}\gamma_{jk})$
(66)

with the 1-RDM elements:

$\displaystyle\gamma_{ij}=\sum_{a\in
occ}\frac{1}{(v_{a}-\varepsilon_{i})(v_{a}-\varepsilon_{j})}\frac{1}{\frac{\partial
A(v_{a})}{\partial v_{a}}}.$ (67)

We will now work out the 2nd-order energy correction. Our strategy is to work
out the unrestricted expression first, then take the limit that the two spin
projections are identical to get the restricted version.

The ABA ground state couples with single and double excitations through the
Hamiltonian (51), so (63) becomes

$\displaystyle E^{(2)}_{U,0}$ $\displaystyle=\sum_{\begin{subarray}{c}a\in
occ\\\ p\in
virt\end{subarray}}\sum_{\sigma}\frac{1}{v^{\sigma}_{a}-v^{\sigma}_{p}}\frac{|\braket{\\{v\\}}{\hat{H}_{C}}{\\{v\\}^{p\sigma}_{a\sigma}}|^{2}}{\braket{\\{v\\}}{\\{v\\}}\braket{\\{v\\}^{p\sigma}_{a\sigma}}{\\{v\\}^{p\sigma}_{a\sigma}}}$
$\displaystyle+\sum_{\begin{subarray}{c}a<b\in occ\\\ p<q\in
virt\end{subarray}}\sum_{\sigma}\frac{1}{v^{\sigma}_{a}+v^{\sigma}_{b}-v^{\sigma}_{p}-v^{\sigma}_{q}}\frac{|\braket{\\{v\\}}{\hat{H}_{C}}{\\{v\\}^{p\sigma
q\sigma}_{a\sigma
b\sigma}}|^{2}}{\braket{\\{v\\}}{\\{v\\}}\braket{\\{v\\}^{p\sigma
q\sigma}_{a\sigma b\sigma}}{\\{v\\}^{p\sigma q\sigma}_{a\sigma b\sigma}}}$
$\displaystyle+\sum_{\begin{subarray}{c}ab\in occ\\\ pq\in
virt\end{subarray}}\frac{1}{v^{\alpha}_{a}+v^{\beta}_{b}-v^{\alpha}_{p}-v^{\beta}_{q}}\frac{|\braket{\\{v\\}}{\hat{H}_{C}}{\\{v\\}^{p\alpha
q\beta}_{a\alpha
b\beta}}|^{2}}{\braket{\\{v\\}}{\\{v\\}}\braket{\\{v\\}^{p\alpha
q\beta}_{a\alpha b\beta}}{\\{v\\}^{p\alpha q\beta}_{a\alpha b\beta}}}.$ (68)

With the TDM elements computed in section (II.3), this becomes

$\displaystyle E^{(2)}_{U,0}$ $\displaystyle=\sum_{\begin{subarray}{c}a\in
occ\\\ p\in
virt\end{subarray}}\sum_{\sigma}\frac{1}{v^{\sigma}_{a}-v^{\sigma}_{p}}\frac{(D^{p\sigma}_{a\sigma})^{2}}{\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}\frac{\partial
A(v^{\sigma}_{p})}{\partial v^{\sigma}_{p}}}$
$\displaystyle+\frac{1}{2}\sum_{\begin{subarray}{c}a<b\in occ\\\ p<q\in
virt\end{subarray}}\sum_{\sigma}\frac{1}{v^{\sigma}_{a}+v^{\sigma}_{b}-v^{\sigma}_{p}-v^{\sigma}_{q}}\frac{(W^{\sigma\sigma}_{abpq}-W^{\sigma\sigma}_{abqp})^{2}}{\frac{\partial
A(v^{\sigma}_{a})}{\partial v^{\sigma}_{a}}\frac{\partial
A(v^{\sigma}_{b})}{\partial v^{\sigma}_{b}}\frac{\partial
A(v^{\sigma}_{p})}{\partial v^{\sigma}_{p}}\frac{\partial
A(v^{\sigma}_{q})}{\partial v^{\sigma}_{q}}}$
$\displaystyle+\sum_{\begin{subarray}{c}ab\in occ\\\ pq\in
virt\end{subarray}}\frac{1}{v^{\alpha}_{a}+v^{\beta}_{b}-v^{\alpha}_{p}-v^{\beta}_{q}}\frac{(W^{\alpha\beta}_{abpq})^{2}}{\frac{\partial
A(v^{\alpha}_{a})}{\partial v^{\alpha}_{a}}\frac{\partial
A(v^{\beta}_{b})}{\partial v^{\beta}_{b}}\frac{\partial
A(v^{\alpha}_{p})}{\partial v^{\alpha}_{p}}\frac{\partial
A(v^{\beta}_{q})}{\partial v^{\beta}_{q}}}$ (69)

where the coupling matrix elements for the singles are

$\displaystyle D^{p\sigma}_{a\sigma}=t^{\sigma}_{ap}+\sum_{c\in
occ}\left(\frac{W^{\sigma\sigma}_{acpc}-W^{\sigma\sigma}_{accp}}{\frac{\partial
A(v^{\sigma}_{c})}{\partial
v^{\sigma}_{c}}}+\frac{W^{\sigma\bar{\sigma}}_{acpc}}{\frac{\partial
A(v^{\bar{\sigma}}_{c})}{v^{\bar{\sigma}}_{c}}}\right)$ (70)

and the integrals have been transformed directly:

$\displaystyle t^{\sigma}_{ap}$
$\displaystyle=\sum_{ij}\frac{h_{ij}}{(v^{\sigma}_{a}-\varepsilon^{\sigma}_{i})(v^{\sigma}_{p}-\varepsilon^{\sigma}_{j})}$
(71) $\displaystyle W^{\sigma\tau}_{abpq}$
$\displaystyle=\sum_{ijkl}\frac{V_{ijkl}}{(v^{\sigma}_{a}-\varepsilon^{\sigma}_{i})(v^{\tau}_{v}-\varepsilon^{\tau}_{j})(v^{\sigma}_{p}-\varepsilon^{\sigma}_{k})(v^{\tau}_{q}-\varepsilon^{\tau}_{l})}.$
(72)

The restricted expression is obtained in the limit that the two sets of
parameters are identical. Specifically,

$\displaystyle E^{(2)}_{R,0}$ $\displaystyle=2\sum_{a\in occ}\sum_{p\in
virt}\frac{1}{v_{a}-v_{p}}\frac{(D^{p}_{a})^{2}}{\frac{\partial
A(v_{a})}{\partial v_{a}}\frac{\partial A(v_{p})}{\partial v_{p}}}$
$\displaystyle+\sum_{ab\in occ}\sum_{pq\in
virt}\frac{1}{v_{a}+v_{b}-v_{p}-v_{q}}\frac{W_{abpq}(2W_{abpq}-W_{abqp})}{\frac{\partial
A(v_{a})}{\partial v_{a}}\frac{\partial A(v_{b})}{\partial
v_{b}}\frac{\partial A(v_{p})}{\partial v_{p}}\frac{\partial
A(v_{q})}{\partial v_{q}}}$ (73)

where the matrix element for the singles contribution is

$\displaystyle D^{p}_{a}=t_{ap}+\sum_{b\in
occ}\frac{2W_{abpb}-W_{abbp}}{\frac{\partial A(v_{b})}{\partial v_{b}}}$ (74)

and the integrals are

$\displaystyle t_{ap}$
$\displaystyle=\sum_{ij}\frac{h_{ij}}{(v_{a}-\varepsilon_{i})(v_{p}-\varepsilon_{j})}$
(75) $\displaystyle W_{abpq}$
$\displaystyle=\sum_{ijkl}\frac{V_{ijkl}}{(v_{a}-\varepsilon_{i})(v_{b}-\varepsilon_{j})(v_{p}-\varepsilon_{k})(v_{q}-\varepsilon_{l})}.$
(76)

The expression (III) is identical to the restricted MP2 expression MEST
except that the orbitals are not normalized and there are non-zero
contributions from single excitations.

##  IV PT expressions: ABA basis

Expression (III) was computed using the Hamiltonian (51) in the primitive
basis, using transition density matrix elements computed with the ABA. The
integrals $h_{ij}$ and $V_{ijkl}$ end up being transformed directly. To be
certain that (III) is correct, we can compute the 2nd order perturbation in
another manner, by transforming (51) to the basis of ABA quasiparticles (the
equivalent of molecular orbitals) and using particle-hole excitations to
generate the relevant excited states. Notice that with the _complete_ set of
solutions to the BAE we get a linear transformation of the creation operators:

$\displaystyle\begin{pmatrix}a^{\dagger}(v^{\sigma}_{1})\\\ \vdots\\\
a^{\dagger}(v^{\sigma}_{N})\end{pmatrix}=\begin{pmatrix}\frac{1}{v^{\sigma}_{1}-\varepsilon^{\sigma}_{1}}&\dots&\frac{1}{v^{\sigma}_{1}-\varepsilon^{\sigma}_{N}}\\\
\vdots&\ddots&\vdots\\\
\frac{1}{v^{\sigma}_{N}-\varepsilon^{\sigma}_{1}}&\dots&\frac{1}{v^{\sigma}_{N}-\varepsilon^{\sigma}_{N}}\end{pmatrix}\begin{pmatrix}a^{\dagger}_{1\sigma}\\\
\vdots\\\ a^{\dagger}_{N\sigma}\end{pmatrix}$ (77)

or

$\displaystyle
a^{\dagger}(\textbf{v}^{\sigma})=C^{\sigma}\textbf{a}^{\dagger}_{\sigma}.$
(78)

The matrix $C$ has a known explicit inverse,

$\displaystyle
C^{-1}_{Ii}=(v_{I}-\varepsilon_{i})\prod_{\begin{subarray}{c}K\neq I\\\ k\neq
i\end{subarray}}\frac{(\varepsilon_{i}-v_{K})(v_{I}-\varepsilon_{k})}{(v_{I}-v_{K})(\varepsilon_{i}-\varepsilon_{k})}$
(79)

but this expression is not optimal numerically as floating point precision
will be lost rather quickly. As a result, $C$ will be inverted numerically.
Upper case letters have been used to label parameters in the ABA basis.
Inverting (78) gives

$\displaystyle\textbf{a}^{\dagger}_{\sigma}=(C^{\sigma})^{-1}a^{\dagger}(\textbf{v}^{\sigma}),$
(80)

so that the Coulomb Hamiltonian can then be written in the ABA parameters

$\displaystyle\hat{H}_{C}=\sum_{IJ}\sum_{\sigma}\tilde{h}^{\sigma}_{IJ}a^{\dagger}(v^{\sigma}_{I})a(v^{\sigma}_{J})+\frac{1}{2}\sum_{IJKL}\sum_{\sigma\tau}\tilde{V}^{\sigma\tau}_{IJKL}a^{\dagger}(v^{\sigma}_{I})a^{\dagger}(v^{\tau}_{J})a(v^{\tau}_{L})a(v^{\sigma}_{K})$
(81)

where the integrals have been transformed

$\displaystyle\tilde{h}^{\sigma}_{IJ}$
$\displaystyle=\sum_{ij}h_{ij}(C^{\sigma}_{Ii})^{-1}(C^{\sigma}_{Jj})^{-1}$
(82) $\displaystyle\tilde{V}^{\sigma\tau}_{IJKL}$
$\displaystyle=\sum_{ijkl}V_{ijkl}(C^{\sigma}_{Ii})^{-1}(C^{\tau}_{Jj})^{-1}(C^{\sigma}_{Kk})^{-1}(C^{\tau}_{Ll})^{-1}.$
(83)

Excited ABA states are generated by acting on the ground state with the ABA
operators. Specifically the single excitation from the $a$th occupied orbital
with spin $\mu$ to the $p$th virtual with spin $\nu$ is written

$\displaystyle
a^{\dagger}(v^{\nu}_{p})a(v^{\mu}_{a})\ket{\\{v\\}}=\frac{\partial
A(v^{\mu}_{a})}{\partial v^{\mu}_{a}}\ket{\\{v\\}^{p\nu}_{a\mu}}$ (84)

and likewise double excitations are generated

$\displaystyle
a^{\dagger}(v^{\kappa}_{p})a^{\dagger}(v^{\eta}_{q})a(v^{\nu}_{b})a(v^{\mu}_{a})\ket{\\{v\\}}=\frac{\partial
A(v^{\mu}_{a})}{\partial v^{\mu}_{a}}\frac{\partial A(v^{\nu}_{b})}{\partial
v^{\nu}_{b}}\ket{\\{v\\}^{p\kappa q\eta}_{a\mu b\nu}}.$ (85)

Again, the 2nd order energy correction will have contributions from singles
and doubles

$\displaystyle E^{(2)}_{U,0}$ $\displaystyle=\sum_{\begin{subarray}{c}a\in
occ\\\ p\in
virt\end{subarray}}\sum_{\mu\nu}\frac{1}{v^{\mu}_{a}-v^{\nu}_{p}}\frac{|\braket{\\{v\\}}{\hat{H}_{C}}{\\{v\\}^{p\nu}_{a\mu}}|^{2}}{\braket{\\{v\\}}{\\{v\\}}\braket{\\{v\\}^{p\nu}_{a\mu}}{\\{v\\}^{p\nu}_{a\mu}}}$
$\displaystyle+\frac{1}{4}\sum_{\begin{subarray}{c}ab\in occ\\\ pq\in
virt\end{subarray}}\sum_{\mu\nu\kappa\eta}\frac{1}{v^{\mu}_{a}+v^{\nu}_{b}-v^{\kappa}_{p}-v^{\eta}_{q}}\frac{|\braket{\\{v\\}}{\hat{H}_{C}}{\\{v\\}^{p\kappa
q\eta}_{a\mu b\nu}}|^{2}}{\braket{\\{v\\}}{\\{v\\}}\braket{\\{v\\}^{p\kappa
q\eta}_{a\mu b\nu}}{\\{v\\}^{p\kappa q\eta}_{a\mu b\nu}}}$ (86)

for which we will have to work out TDM elements. This may be done by Wick’s
theorem as the ABA operators yield orbitals that are orthogonal but not
normal. The coupling between the ABA ground state and single excitations has a
contribution from 1-electron operators:

$\displaystyle\braket{\\{v\\}}{a^{\dagger}(v^{\sigma}_{I})a(v^{\sigma}_{J})}{\\{v\\}^{p\nu}_{a\mu}}=\delta_{aI}\delta_{pJ}\delta_{\nu\sigma}\delta_{\mu\sigma}\frac{\partial
A(v^{\nu}_{p})}{\partial v^{\nu}_{p}}\braket{\\{v\\}}{\\{v\\}}$ (87)

and a more complicated contribution from 2-electron operators

$\displaystyle\braket{\\{v\\}}{a^{\dagger}(v^{\sigma}_{I})a^{\dagger}(v^{\tau}_{J})a(v^{\tau}_{L})a(v^{\sigma}_{K})}{\\{v\\}^{p\nu}_{a\mu}}$
$\displaystyle=\frac{\partial A(v^{\nu}_{p})}{\partial
v^{\nu}_{p}}\frac{\partial A(v^{\tau}_{L})}{\partial
v^{\tau}_{L}}\braket{\\{v\\}}{\\{v\\}}\delta_{pK}\delta_{JL}\delta_{aI}\delta_{\nu\sigma}\delta_{\mu\sigma}\sum_{c\in
occ}\delta_{cJ}$ $\displaystyle+\frac{\partial A(v^{\nu}_{p})}{\partial
v^{\nu}_{p}}\frac{\partial A(v^{\sigma}_{K})}{\partial
v^{\sigma}_{K}}\braket{\\{v\\}}{\\{v\\}}\delta_{pL}\delta_{IK}\delta_{aJ}\delta_{\nu\tau}\delta_{\mu\tau}\sum_{c\in
occ}\delta_{cI}$ $\displaystyle-\frac{\partial A(v^{\nu}_{p})}{\partial
v^{\nu}_{p}}\frac{\partial A(v^{\sigma}_{L})}{\partial
v^{\sigma}_{L}}\braket{\\{v\\}}{\\{v\\}}\delta_{pK}\delta_{IL}\delta_{aJ}\delta_{\nu\sigma}\delta_{\mu\tau}\delta_{\sigma\tau}\sum_{c\in
occ}\delta_{cI}$ $\displaystyle-\frac{\partial A(v^{\nu}_{p})}{\partial
v^{\nu}_{p}}\frac{\partial A(v^{\sigma}_{K})}{\partial
v^{\sigma}_{K}}\braket{\\{v\\}}{\\{v\\}}\delta_{pL}\delta_{JK}\delta_{aI}\delta_{\nu\tau}\delta_{\mu\sigma}\delta_{\sigma\tau}\sum_{c\in
occ}\delta_{cJ}.$ (88)

The only non-vanishing coupling between the ABA ground state and double
excitations is

$\displaystyle\braket{\\{v\\}}{a^{\dagger}(v^{\sigma}_{I})a^{\dagger}(v^{\tau}_{J})a(v^{\tau}_{L})a(v^{\sigma}_{K})}{\\{v\\}^{p\kappa
q\eta}_{a\mu b\nu}}=\frac{\partial A(v^{\kappa}_{p})}{\partial
v^{\kappa}_{p}}\frac{\partial A(v^{\eta}_{q})}{\partial
v^{\eta}_{q}}\braket{\\{v\\}}{\\{v\\}}\times$
$\displaystyle\times(\delta_{aI}\delta_{bJ}\delta_{\mu\sigma}\delta_{\nu\tau}-\delta_{aJ}\delta_{bI}\delta_{\mu\tau}\delta_{\nu\sigma})(\delta_{pK}\delta_{qL}\delta_{\kappa\sigma}\delta_{\eta\tau}-\delta_{pL}\delta_{qK}\delta_{\kappa\tau}\delta_{\eta\sigma}).$
(89)

The final expression for the unrestricted 2nd order correction is thus

$\displaystyle E^{(2)}_{U,0}$ $\displaystyle=\sum_{\begin{subarray}{c}a\in
occ\\\ p\in
virt\end{subarray}}\sum_{\mu}\frac{1}{v^{\mu}_{a}-v^{\mu}_{p}}\frac{\partial
A(v^{\mu}_{p})}{\partial v^{\mu}_{p}}\frac{\partial A(v^{\mu}_{a})}{\partial
v^{\mu}_{a}}(\tilde{D}^{p\mu}_{a\mu})^{2}$
$\displaystyle+\frac{1}{2}\sum_{\begin{subarray}{c}ab\in occ\\\ pq\in
virt\end{subarray}}\sum_{\mu\nu}\frac{\tilde{V}^{\mu\nu}_{pqab}\tilde{V}^{\mu\nu}_{pqab}}{v^{\mu}_{a}+v^{\nu}_{b}-v^{\mu}_{p}-v^{\nu}_{q}}\frac{\partial
A(v^{\mu}_{p})}{\partial v^{\mu}_{p}}\frac{\partial A(v^{\nu}_{q})}{\partial
v^{\nu}_{q}}\frac{\partial A(v^{\mu}_{a})}{\partial v^{\mu}_{a}}\frac{\partial
A(v^{\nu}_{b})}{\partial v^{\nu}_{b}}$
$\displaystyle-\frac{1}{2}\sum_{\begin{subarray}{c}ab\in occ\\\ pq\in
virt\end{subarray}}\sum_{\mu}\frac{\tilde{V}^{\mu\mu}_{pqab}\tilde{V}^{\mu\mu}_{pqba}}{v^{\mu}_{a}+v^{\mu}_{b}-v^{\mu}_{p}-v^{\mu}_{q}}\frac{\partial
A(v^{\mu}_{p})}{\partial v^{\mu}_{p}}\frac{\partial A(v^{\mu}_{q})}{\partial
v^{\mu}_{q}}\frac{\partial A(v^{\mu}_{a})}{\partial v^{\mu}_{a}}\frac{\partial
A(v^{\mu}_{b})}{\partial v^{\mu}_{b}}$ (90)

where the singles contribution is

$\displaystyle\tilde{D}^{p\mu}_{a\mu}=\tilde{h}^{\mu}_{pa}+\sum_{c\in
occ}\left(\sum_{\sigma}\tilde{V}^{\sigma\mu}_{cpca}\frac{\partial
A(v^{\sigma}_{c})}{\partial
v^{\sigma}_{c}}-\tilde{V}^{\mu\mu}_{cpac}\frac{\partial
A(v^{\mu}_{c})}{\partial v^{\mu}_{c}}\right).$ (91)

Taking the restricted limit, we arrive at

$\displaystyle E^{(2)}_{R,0}$ $\displaystyle=2\sum_{\begin{subarray}{c}a\in
occ\\\ p\in virt\end{subarray}}\frac{1}{v_{a}-v_{p}}\frac{\partial
A(v_{p})}{\partial v_{p}}\frac{\partial A(v_{a})}{\partial
v_{a}}\left(\tilde{h}_{pa}+\sum_{c\in
occ}(2\tilde{V}_{cpca}-\tilde{V}_{cpac})\frac{\partial A(v_{c})}{\partial
v_{c}}\right)^{2}$ $\displaystyle+\sum_{\begin{subarray}{c}ab\in occ\\\ pq\in
virt\end{subarray}}\frac{(2\tilde{V}_{pqab}-\tilde{V}_{pqba})\tilde{V}_{pqab}}{v_{a}+v_{b}-v_{p}-v_{q}}\frac{\partial
A(v_{p})}{\partial v_{p}}\frac{\partial A(v_{q})}{\partial
v_{q}}\frac{\partial A(v_{a})}{\partial v_{a}}\frac{\partial
A(v_{b})}{\partial v_{b}}.$ (92)

Equation (IV) is quite similar to (III), except that the integrals have been
transformed in the opposite manner.

##  V Discussion

We derived two expressions, (III) and (IV), for the restricted 2nd order RSPT
energy correction. Both were implemented, and found to give numerically the
same results. Unfortunately these results were ambiguous. The ABA parameters
$\\{\varepsilon\\},g$ may be changed without affecting the mean-field energy
(66), though they have a large effect on the 2nd-order energy. In addition,
unlike the case for MP2, the singles couple to the ground state through
$\hat{H}_{C}$. To study these effects we tried minimizing the energy
functional

$\displaystyle F=E^{ABA}_{R}+\sum_{\begin{subarray}{c}a\in occ\\\ p\in
virt\end{subarray}}|\braket{\\{v\\}}{\hat{H}_{C}}{\\{v\\}^{p}_{a}}|^{2}$ (93)

which leads to numerically the same mean-field energy, and the singles do
vanish. _But so do the doubles_ and thus the 2nd order RSPT correction is
numerically zero. The $N+1$ variables $\\{\varepsilon\\},g$ are enough to
reproduce the ground state RHF energy as all that is required is an idempotent
1-RDM in the correct basis. However, they are not enough to reproduce the
entire RHF spectrum and hence do not reproduce the MP2 perturbative
correction. To fix this, one could employ a Hylleraas functionalhylleraas  to
solve for the perturbative correction variationally. Solving the equations
explicitly on paper would of course lead to the same expressions (III) and
(IV). Instead one should consider variationally solving for the coefficients
$\\{c_{ki}\\}$

$\displaystyle\ket{\Psi^{(1)}_{k}}=\sum_{i}c_{ki}\ket{\Psi^{(0)}_{i}}$ (94)

but this amounts to the same number of parameters as an orbital optimization.
Including an orbital optimization would necessarily reproduce MP2 as the RHF
solution is defined by the optimal orbitals. Similarly, we could solve the CI
singles problem for ABA states as a reference wavefunction, but again this
amounts to employing $N^{2}$ parameters. Therefore we do not consider it
productive to demonstrate numerically.

##  VI Conclusion

Perturbative corrections for the ABA for individual electrons were calculated.
The first expression (III) was computed in the primitive basis using TDM
elements obtained from the ABA. The second expression (IV) was computed in the
ABA basis using particle-hole excitations. Both expressions are in close
analogy with MP2 corrections to RHF. Numerically computing these corrections
lead to ambiguous results that could, in principle, be corrected with an
orbital optimization or a CI singles solution. Going forward to pairs of
electrons, we understand that an orbital optimization is required along with
the ABA wavefunction. Care will be required when computing perturbative
corrections.

###### Acknowledgements.

We are grateful for support form NSERC. This research was enabled by support
from CalculQuébec and Compute Canada.

## References

  * (1) G. K. L. Chan and S. Sharma.  Annual Review of Physical Chemistry, 62:465, 2011.
  * (2) G. H. Booth, A. Gruneis, G. Kresse, and A. Alavi.  Nature, 493:431, 2013.
  * (3) A. A. Holmes, C. J. Umrigar, and S. Sharma.  The Journal of Chemical Physics, 147:164111, 2017.
  * (4) J. Li, M. Otten, A. A. Holmes, S. Sharma, and C. J. Umrigar.  The Journal of Chemical Physics, 149:214110, 2018.
  * (5) A. C. Hurley, J. Lennard-Jones, and J. A. Pople.  Proceedings of the Royal Society, A220:446, 1953.
  * (6) A. J. Coleman.  Journal of Mathematical Physics, 6:1425, 1965.
  * (7) D. M. Silver.  The Journal of Chemical Physics, 50:5108, 1969.
  * (8) D. M. Silver, E. L. Mehler, and K. Ruedenberg.  The Journal of Chemical Physics, 52:1174, 1970.
  * (9) A. J. Coleman.  International Journal of Quantum Chemistry, 63:23, 1997.
  * (10) P. R. Surján.  An Introduction to the Theory of Geminals.  Springer, Berlin, 1999.
  * (11) M. Kobayashi, A. Szabados, H. Nakai, and P. Surjan.  Journal of Chemical Theory and Computation, 6:2024, 2010.
  * (12) P. R. Surján, Á Szabados, P. Jeszenski, and T. Zoboki.  Journal of Mathematical Chemistry, 50:534, 2012.
  * (13) E. Neuscamman.  Physical Review Letters, 109:203001, 2012.
  * (14) P. A. Limacher, P. W. Ayers, P. A. Johnson, S. De Baerdemacker, D. Van Neck, and P. Bultinck.  Journal of Chemical Theory and Computation, 9:1394, 2013.
  * (15) P. A. Johnson, P. W. Ayers, P. A. Limacher, S. De Baerdemacker, D. Van Neck, and P. Bultinck.  Computational and Theoretical Chemistry, 1003:101, 2013.
  * (16) T. Stein, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 140:214113, 2014.
  * (17) K. Boguslawski, P. Tecmer, P. W. Ayers, P. Bultinck, S. De Baerdemacker, and D. Van Neck.  Physical Review B, 98:201106(R), 2014.
  * (18) K. Boguslawski, P. Tecmer, P. Bultinck, S. De Baerdemacker, D. Van Neck, and P. W. Ayers.  Journal of Chemical Theory and Computation, 10:4873, 2014.
  * (19) K. Boguslawski, P. Tecmer, P. A. Limacher, P. A. Johnson, P. W. Ayers, P. Bultinck, S. De Baerdemacker, and D. Van Neck.  Journal of Chemical Theory and Computation, 140:214114, 2014.
  * (20) P. Tecmer, K. Boguslawski, P. A. Johnson, M. Chan, T. Verstraelen, and P. W. Ayers.  Journal of Physical Chemistry, A118:9058, 2014.
  * (21) T. M. Henderson, G. E. Scuseria, J. Dukelsky, A. Signoracci, and T. Duguet.  Physical Review C, 89:054305, 2014.
  * (22) T. M. Henderson, I. W. Bulik, T. Stein, and G. E. Scuseria.  The Journal of Chemical Physics, 141:244104, 2014.
  * (23) J. J. Shepherd, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 144:094112, 2014.
  * (24) I. W. Bulik, T. M. Henderson, and G. E. Scuseria.  Journal of Chemical Theory and Computation, 11:3171, 2015.
  * (25) E. Pastorczak and K. Pernal.  Physical Chemistry Chemical Physics, 17:8622, 2015.
  * (26) T. M. Henderson and G. E. Scuseria.  The Journal of Chemical Physics, 151:051101, 2019.
  * (27) A. Khamoshi, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 151:184103, 2019.
  * (28) T. M. Henderson and G. E. Scuseria.  The Journal of Chemical Physics, 153:084111, 2020.
  * (29) R. Dutta, T. M. Henderson, and G. E. Scuseria.  Journal of Chemical Theory and Computation, 16:6358, 2020.
  * (30) G. Harsha, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 153:124115, 2020.
  * (31) A. Khamoshi, F. A. Evangelista, and G. E. Scuseria.  Quantum Science and Technology, 6:014004, 2020.
  * (32) A. Khamoshi, G. P. Chen, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 154:074113, 2021.
  * (33) R. Dutta, G. P. Chen, T. M. Henderson, and G. E. Scuseria.  The Journal of Chemical Physics, 154:114112, 2021.
  * (34) A. Marie, F. Kossoki, and P.-F. Loos.  arXiv:2106.11305.
  * (35) H. Bethe.  Zeitschrift für Physik, 71:205, 1931.
  * (36) L. D. Faddeev and L. A. Takhtajan.  Soviet Scientific Reviews, C 1:107, 1981.
  * (37) J. Bardeen, L. N. Cooper, and J. R. Schrieffer.  Physical Review, 106:162, 1957.
  * (38) J. Bardeen, L. N. Cooper, and J. R. Schrieffer.  Physical Review, 108:1175, 1957.
  * (39) R. W. Richardson.  Physics Letters, 3:277, 1963.
  * (40) R. W. Richardson and N. Sherman.  Nuclear Physics, 52:221, 1964.
  * (41) R. W. Richardson.  Journal of Mathematical Physics, 6:1034, 1965.
  * (42) M. Gaudin.  Journal de Physique, 37:1087, 1976.
  * (43) J. Dukelsky, S. Pittel, and G. Sierra.  Reviews of Modern Physics, 76:643, 2004.
  * (44) P. A. Johnson, C.-É. Fecteau, F. Berthiaume, S. Cloutier, L. Carrier, M. Gratton, P. Bultinck, S. De Baerdemacker, D. Van Neck, P. Limacher, and P. W. Ayers.  The Journal of Chemical Physics, 153:104110, 2020.
  * (45) C.-É. Fecteau, H. Fortin, S. Cloutier, and P. A. Johnson.  The Journal of Chemical Physics, 153:164117, 2020.
  * (46) C.-É. Fecteau, F. Berthiaume, M. Khalfoun, and P. A. Johnson.  Journal of Mathematical Chemistry, 59:289, 2021.
  * (47) P. A. Johnson, H. Fortin, S. Cloutier, and C.-É. Fecteau.  The Journal of Chemical Physics, 154:124125, 2021.
  * (48) V. E. Korepin, N. M. Bogoliubov, and A. G. Izergin.  Quantum Inverse Scattering Method and Correlation Functions.  Cambrige University Press, Cambridge, 1993.
  * (49) R. Baxter.  Exactly Solved Models in Statistical Mechanics.  Dover, New York, 2007.
  * (50) L. Carrier, C.-É. Fecteau, and P. A. Johnson.  International Journal of Quantum Chemistry, page e26255, 2020.
  * (51) T. Helgaker, P. Jørgensen, and J. Olsen.  Molecular Electronic Structure Theory.  Wiley, Chichester, 2004.
  * (52) H.-Q. Zhou, J. Links, R. H. McKenzie, and M. D. Gould.  Physical Review B, 65:060502, 2002.
  * (53) A. Faribault, P. Calabrese, and J.-S. Caux.  Physical Review B, 77:064503, 2008.
  * (54) A. Faribault, P. Calabrese, and J.-S. Caux.  Physical Review B, 81:174507, 2010.
  * (55) G. Gorohovsky and E. Bettelheim.  Physical Review B, 84:224503, 2011.
  * (56) E. A. Hylleraas.  Zeitschrift für Physik, 65:209, 1930.