# Generalized Slater-Jastrow trial functional: Application to high order
correlations in the two and three dimensional electron gas

James C. Porter<EMAIL_ADDRESS>Department of Physics and Astronomy,
Eastern Michigan University,
Ypsilanti, MI 48197

###### Abstract

The GSJ trial functional is a modification of the Slater-Jastrow functional
where, effectively, the argument of the Jastrow factor can be momentum
dependent. The associated Euler equations, which provide an alternative to the
usual approach via perturbation theory, are solved to yield the first three
terms of the electron correlation energies for the two and three dimensional
electron gas at high density. In two dimensions we find
$\epsilon_{c}=-0.4039-0.50\,r_{s}\,\log r_{s}$ plus an estimated
$-0.87\,r_{s}$, and for three dimensions
$\epsilon_{c}=0.0622\,\log\,r_{s}-0.082+0.059\,r_{s}\,\log r_{s}$ plus an
estimated $-0.015\,r_{s}$ where $r_{s}$ is the applicable Seitz radius in each
case. The present paper extends the findings of the similarly titled 2009 work
(Phys. Rev. B 80) in two respects. First the two dimensional as well as the
three dimensional electron gas is treated and second high order coulomb
correlations are obtained in each case, furnishing an alternative to the
findings of Gell-Mann and Breuckner (Phys. Rev. 106 (1957)), Carr and
Maradudin (Phys. Rev. 133 (1964)) and Rajagopal and Kimball (Phys. Rev. B 15
(1977)). The latter development is due to the discovery of an analytic
solution for the lowest order two electron functions, obtained previously only
numerically, and a procedure which isolates the logarithmic divergence in the
three dimensional case without the need to set the Seitz radius to zero. Other
points of interest include a more complete accounting of background energy
cancelation in the second section, and the evaluation of a 12-dimensional
integral contributing to the high order coulomb correlation energy in Appendix
B.

Hartree-Fock, jellium, high orderr correlation energy

###### pacs:

71.10.Ca 71.15.Nc 31.15.xt 02.30Rz

## Introduction

Variational approximations to the wave function satisfying the many body
Schrodinger equation can be obtained from a great variety of comparison
functionalsFoulkes et al. (2001). A subset of these, consists of a series
which begins with the Hartree Fock (HF) functional, an antisymmetrized product
of $N$ one electron functions often written as a determinant,Slater (1960) and
continues with the GSJ functional,Porter (2009) constructed from both one and
two electron functions. The direct terms of the functionals in this series
approximate the many electron wave function in a manner suggested by the Mayer
cluster expansion,Mayer and Mayer (1940) in that the Hartree two electron
density matrix is approximated as the product of two one electron matrices
while the direct terms of the GSJ functional approximate the three electron
density matrix by terms representing the three ways of combining one and two
electron matrices less twice the product of three one electron matrices.Porter
(1965a)

Like the one-electron functions of the HF functional, the one and two-electron
functions of the GSJ functional are a priory unconstrained, save for simple
orthogonality and normalization requirements, and are determined by associated
integro-differential Euler equations. Also, like the HF functional, the
coordinates (position and momentum) of a given electron appear only once in
each additive term of the GSJ functional.

The Euler equations associated with the GSJ functional are tractable in the
case of the homogeneous electron gas in a background of compensating positive
charge (possibly with the addition of an externally imposed electric
field.Oberman et al. (1962)) Here plane waves satisfy the one electron Euler
equation, as in the HF theory,Fetter and Walecka (1971) and the two electron
Euler equation has an approximate solution in closed form which allows
evaluation of the order $r_{s}\,\log\,r_{s}$ term of the correlation energy,
$r_{s}$ being the applicable Seitz radius. It is instructive to compare these
results with those of perturbation theory. The latter assumes that the system
energy is expressible as a power series in some small parameter measuring the
ratio of potential to kinetic energies. This assumption does not hold in the
case of the electron gas because of the energies logarithmic dependence on the
assumed small Seitz radius so that complicatedGell-Mann and Breuckner (1957);
Carr and Maradudin (1964) and possibly problematicFarmelo (2009) work-arounds
are required.

The following sections cover the construction of the GSJ functional, Euler
equations for the electron gas Hamiltonian, their solution and correlation
energy calculations through the $r_{s}\,\log\,r_{s}$ term, with an estimate
for the order $r_{s}$ term. Results for both the two and three dimensional gas
are presented in parallel prior to the evaluation of the individual
correlation energies.

## GSJ functional and energy integral

Denoting the one and two electron functions by $U(j_{1};i_{1})$ and
$f(j_{1},j_{2};i_{1},i_{2})$, where indices $j_{1},\,j_{2}\cdots$ appearing
prior to the semicolons label momenta and spin, and integers
$i_{1},\,i_{2}\cdots$ following the semicolons show position dependences
$\mathbf{r}_{i_{1}},\,\mathbf{r}_{i_{2}}\cdots$, the GSJ functional is
constructed from the Laplace developmentKorn and Korn (1968) of the N electron
HF determinant, $\Psi_{HF}=\det((U(j;i)))$, by 2n-rowed minors,

$\Psi_{HF}=\sum_{1\leq j_{1}<j_{2}<\cdots j_{2n}\leq N}M(j_{1},j_{2},\cdots
j_{2n};i_{1},i_{2}\cdots i_{2n})D(j_{1},j_{2},\cdots j_{2n};i_{1},i_{2}\cdots
i_{2n}),$ (1)

where $M$ is the minor obtained by retaining only those elements of
$\Psi_{HF}$ found in columns $j_{1},j_{2}.\cdots j_{2n}$ and rows
$i_{1},i_{2},\cdots i_{2n}$ and where the $N\\!\\!-\\!2n$ rowed determinant
$D(j_{1},j_{2},\cdots j_{2n};i_{1},i_{2}\cdots i_{2n})$ is the corresponding
cofactor.

$\Psi_{GSJ}$ is obtained from Eq.(1) upon replacing the minors there by $1/n!$
times antisymmetrized products of $n$ two electron functions and summing over
all distinguishable combinations of the row numbers chosen for the Laplace
developments. A final summation over $n=0,1,\,2,\cdots[(N-1)/2]$, where the
$n=0$ term is just the HF determinant itself, gives the GSJ seriesPorter
(1965b); not

$\displaystyle\Psi_{GSJ}=\sum_{n=0}^{[(N-1)/2]}\frac{1}{n!}\sum_{\stackrel{{\scriptstyle\scriptstyle
1\leq j_{1}<j_{2}<\cdots<j_{2n}\leq N}}{{1\leq i_{1}<i_{2}<\cdots<i_{2n}\leq
N}}}(\hat{A}^{j}\hat{A}_{i}\prod_{k=1}^{n}f(j_{2k-1},j_{2k};i_{2k-1},i_{2k})\frac{}{})$
$\displaystyle\times D(j_{1},\cdots,j_{2n};i_{1},\cdots,i_{2n}),\qquad$ (2)

where $\hat{A}^{j}$ and $\hat{A}_{i}$ are commuting antisymmetrizing operators
acting on momenta and positions respectively.

The Coulomb correlation energy is given by $1/N$ times the difference between
$E_{GSJ}$, the energy integral, and the corresponding HF estimate, where

$\displaystyle E_{GSJ}=(\Psi_{GSJ},\left[-\frac{\hbar^{2}}{2m}\sum_{1\leq
i\leq N}\mathbf{\nabla}_{i}^{2}+\frac{1}{2}\sum_{1\leq i_{1}\neq i_{2}\leq
N}\phi(i_{1},i_{2})\right]\Psi_{GSJ})/(\Psi_{GSJ},\Psi_{GSJ})$
$\displaystyle+E_{background}.$ (3)

As in Slater’s derivation of the HF equations,Slater (1960) the evaluation of
$E_{GSJ}$ is facilitated by the use of Lagrange multipliers
$\\{\epsilon(j)\\}$ to constrain the one and two electron functions. The
constraints imposed are that the one electron functions for an orthonormal set
and that

$\int
d^{D}r_{1}d^{D}r_{2}U^{\dagger}(j_{1};1)U^{\dagger}(j_{2};2)\hat{A}^{j}\hat{A}_{i}f(j_{1},j_{2};1,2)=0,$
(4a) where the spatial integrals may be over either two or three dimensions,
$D=2$ or $D=3$. Since $f(j_{1},j_{2};1,2)\equiv f(j_{2},j_{1};2,1)$, Eq.(3a)
can be written more simply in terms of $F(j_{1},j_{2};1,2)\equiv
2f(j_{1},j_{2};1,2)$ as $\int
d^{D}r_{1}d^{D}r_{2}U^{\dagger}(j_{1};1)U^{\dagger}(j_{2};2)\hat{A}^{j}F(j_{1},j_{2};1,2)=0.$
(4b)

Terms, where $n$ in Eq.(2) is much less than $N$ but also much greater than
unity, dominate in Eq.(3), where the large numeric factors which appear in
both numerator and denominator reach broad maxima. It can be shownPorter
(1965a, b) that the quotient of these multipliers can be dropped, and the GSJ
series replaced by just its first two terms, provided that one inserts a
weighting factor $(1+\lambda(j_{i}))^{-1}$ under every sum over a momentum
$j_{i}$, where

$\lambda(j_{1})\equiv\sum_{j_{2}}(F(j_{1},j_{2};1,2),\hat{A}^{j}F(j_{1},j_{2};1,2))$
(5)

is assumed small compared to unity. Because we shall evaluate the correlation
series exactly only through its $r_{s}\,\log r_{s}$ term, where $r_{s}$ is the
applicable Seitz radius, weighting factors may be dropped entirely from the
Euler equations. Additionally it is only necessary to include their
contributions to the nominally infinite terms of $E_{GSJ}$ (which must sum to
zero) and to the HF portion of the kinetic energy, so that Eq.(3) givesPorter
(1965a, b)

$\displaystyle
E_{GSJ}=\sum_{j_{1}}(1-\lambda(j_{1}))\left(U(j_{1};1),(-\frac{\hbar^{2}}{2m}\mathbf{\nabla}_{1}^{2})\,U(j_{1};1)\right)$
$\displaystyle+\frac{1}{2}\sum_{j_{1},\,j_{2}}(1-\lambda(j_{1})-\lambda(j_{2}))\left(U(j_{1};1)U(j_{2};2),\phi(1,2)U(j_{1};1)U(j_{2};2)\frac{}{}\right)\frac{}{}$
$\displaystyle+\sum_{j_{1},\,j_{2},\,j_{3}}\left(\frac{}{}U(j_{3};3)F(j_{1},j_{2};1,2),(\phi(1,3)+\phi(2,3))U(j_{3};3)\hat{A}^{j}F(j_{1},j_{2};1,2)\frac{}{}\right)$
$\displaystyle+\frac{1}{2}\sum_{j_{1},\,j_{2}}\left(F(j_{1},j_{2};1,2),-\frac{\hbar^{2}}{2m}(\mathbf{\nabla}_{1}^{2}+\mathbf{\nabla}^{2}_{2})\hat{A}^{j}F(j_{1},j_{2};1,2)\frac{}{}\right)$
$\displaystyle+\sum_{j_{1},\,j_{2}}\mathsf{Re}\left(U(j_{1};1)U(j_{2};2),\phi(1,2)\hat{A}^{j}F(j_{1},j_{2};1,2)\frac{}{}\right)$
$\displaystyle+\sum_{j_{1},\,j_{2},\,j_{3}}\left(\frac{}{}U(j_{3};3)F(j_{1},j_{2};1,2),\phi(2,3)U(j_{2};2)F(j_{1},j_{3};1,3)\right)$
$\displaystyle-\sum_{j_{1},\,j_{2},\,j_{3}}\left(\frac{}{}U(j_{3};3)F(j_{1},j_{2};1,2),\phi(2,3)(U(j_{3};2)F(j_{1},j_{2};1,3)+U(j_{2};3)F(j_{1},j_{3};1,2))\frac{}{}\right)$
$\displaystyle-\frac{1}{2}\sum_{j_{1},\,j_{2}}\left(U(j_{1};1)U(j_{2};2),\phi(1,2)U(j_{2};1)U(j_{1};2)\frac{}{}\right)\frac{}{}$
$\displaystyle+\frac{1}{2}\sum_{j_{1},\,j_{2}}\left(F(j_{1},j_{2};1,2),\phi(1,2)F(j_{1},j_{2};1,2)\right)+E_{background}.$
(6)

The convention here is that sums run from $1$ to $N$ and that any summand,
considered as a function of the indices $j_{1},\,j_{2},\cdots$, is to be
omitted when any two of these indices become equal.

In Eq.(6) the nominally infinite energy associated with the positively charged
background $E_{b}$ reduces the factor weighting the HF Hartree term to just
$(-\lambda(j_{1})-\lambda(j_{2}))$ so that this and the first triple sum (also
nominally infinite) in Eq.(6) cancel. When one considers infinite terms of
next higher order, where the corresponding weighting for the Hartree term is
now more accurately
$(1-\lambda(j_{1})-\lambda(j_{2})+\lambda(j_{1})\lambda(j_{2})-\lambda(j_{1})^{2}-\lambda(j_{2})^{2})$,
with weighting $(1-\lambda(j_{1})-\lambda(j_{2})-\lambda(j_{3}))$ applied to
the initial triple sum of Eq.(6), and then adding the higher order analogue of
this triple summation, (this with a pair of the two electron functions both
pre- and post-multiplying a sum of four two body potentials), cancelation
again occurs. We shall assume that cancelations continue through all higher
orders.

Finally we note that Eq.(6) has been simplified, in comparison to the form
which must be used to calculate variations,Porter (1965a, b) by dropping terms
proportionate to $\int d^{D}r_{1}U^{\dagger}(j_{1};1)F(j_{1},j_{2};1,2)$. It
is readily shown that putting this latter integral to zero is consistent with
solutions of the following Euler equations.

## Euler equations

Variation with respect to $U^{\dagger}(j_{1};1)$ of $E_{GSJ}$ when just the
initial two terms of the GSJ series are inserted into Eq.(3) gives the
complicated one electron Euler equation discussed in connection with an order
$r_{s}$ term below. Because we do not attempt to evaluate the entirety of
contributions to the correlation series order $r_{s}$ term, it is sufficient
to approximate the one electron equation by the applicable form of the HF
relation,

$(\epsilon(j_{1})+\frac{\hbar^{2}}{2m}\mathbf{\nabla}_{1}^{2})\,U(j_{1};1)=-\sum_{j_{3}}\int
d^{D}r_{3}U^{\dagger}(j_{3};3)\phi(1.3)U(j_{3};1)U(j_{1};3).$ (7)

Variation with respect to $f^{\dagger}(j_{1},j_{2};1,2)$ of the initial two
terms of the GSJ series in Eq.(3) leads to an integro-differential equation
whose differential operators act on the combination
$(U(j_{1};1)U(j_{2};2)+f(j_{1},j_{2};1,2))$ of one and two electron functions.
This, together with Eq.(6), can be manipulated to give the two electron Euler
equationPorter (1965a)

$\displaystyle(\epsilon(j_{1})+\epsilon(j_{2})-\frac{\hbar^{2}}{2m}(\mathbf{\nabla}_{1}^{2}+\mathbf{\nabla}_{2}^{2}))F(j_{1},j_{2};1,2)$
$\displaystyle-\sum_{j_{3}}\int
d^{D}r_{3}U^{\dagger}(j_{3};3)(\tilde{P}(1)\phi(1,3)U(j_{1};1)F(j_{2},j_{3};2,3)+\tilde{P}(2)\phi(2,3)U(j_{2};2)F(j_{1},j_{3};1,3))$
$\displaystyle-\tilde{P}(1)\tilde{P}(2)U(j_{1};1)U(j_{2};2)=\tilde{P}(1)\tilde{P}(2)F(j_{1},j_{2};1,2)$
$\displaystyle-\sum_{j_{3}}\int
d^{D}r_{3}U^{\dagger}(j_{3};3)\left(\phi(1,3)U(j_{3};1)F(j_{2},j_{1};2,3)+\tilde{P}(1)\phi(1,3)U(j_{1};3)F(j_{3},j_{2};1,2)\right.$
$\displaystyle+\left.\phi(2,3)U(j_{3};2)F(j_{1},j_{2};1,3)+\tilde{P}(2)\phi(2,3)U(j_{2};3)F(j_{1},j_{3};1,2)\right),\quad$
(8)

where the idempotent projection operators $\tilde{P}(1)\equiv
1-\sum_{j_{1}^{\prime}}\int
d^{D}r^{\prime}_{1}U^{\dagger}(j^{\prime}_{1};1^{\prime})\tilde{\bigcirc}(1\\!\\!\rightarrow\\!1^{\prime})$,
with $\tilde{\bigcirc}(1\\!\\!\rightarrow\\!1^{\prime})$ operating to replace
functions of $\mathbf{r}_{1}$ by the same function of
$\mathbf{r}^{\prime}_{1}$, subtract from an arbitrary function that part which
is orthogonal to each of the one electron functions.

Eq.(8) has been arranged with direct terms on the left, while on the right are
terms diminished by an additional projection operator or by an exchange of
indices $j_{2}\leftarrow\\!\\!\\!\\!\rightarrow j_{3}$. Still smaller terms,
with exchanges involving three indices, are negligible.

## Approximating the two electron functions

For the spatially homogeneous electron gas
$\phi(1,2)=e^{2}/|\mathbf{r}_{1}-\mathbf{r}_{2}|$, and the one electron
functions are plane waves satisfying periodic boundary conditions over a $D=2$
or $D=3$ dimensional region of side $L$,

$U(j_{1};1)=\frac{1}{\sqrt{L^{D}}}\exp{\left(\imath
K_{D}\mathbf{k}_{1}\cdot\mathbf{r}_{1}\right)}\sigma(j_{1}).$ (9)

Here the Fermi wave number $K_{D}$ is $K_{2}=\sqrt{2\pi N/L^{2}}$ in two
dimensions and $K_{3}=\sqrt[3]{3\pi^{2}N/L^{3}}$ in three dimensions and
$\mathbf{k}_{1}$ denotes the dimensionless wavenumber $\mathbf{k}(j_{1})$. It
is convenient to take $\sigma(j)=\uparrow$ for $1\leq j\leq[N/2]$ and
$\sigma(j)=\downarrow$ for $[N/2]<j\leq[N]$.

Writing the two electron functions in terms of the real-valued Fourier-
transforms ${\cal F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})$,

$F(j_{1},j_{2};1,2)=U(j_{1};1)U(j_{2};2)\int d^{D}q\exp(\imath
K_{D}\,\mathbf{q}\cdot(\mathbf{r}_{1}-\mathbf{r}_{2})){\cal
F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q}),$ (10)

and substituting from Eq.(7) for the Lagrange multipliers, Eq.(8) becomes

$\displaystyle((\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\hat{\mathbf{e}}_{\mathbf{q}}+q){\cal
F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})$ $\displaystyle+\int
d^{D}k_{3}\Theta(1-k_{3})(\Theta(|\,\mathbf{q}+\mathbf{k}_{1}|-1){\cal
F}(\mathbf{k}_{3},\mathbf{k}_{2};\mathbf{q})+\Theta(|\,\mathbf{q}-\mathbf{k}_{2}|-1){\cal
F}(\mathbf{k}_{1},\mathbf{k}_{3};\mathbf{q}))$
$\displaystyle=-\frac{\alpha_{D}R_{D}}{\pi^{D-1}q^{D}}\frac{1}{2}\Theta(|\,\mathbf{q}+\mathbf{k}_{1}|-1)\Theta(|\,\mathbf{q}-\mathbf{k}_{2}|-1)$
$\displaystyle+\frac{\alpha_{D}R_{D}}{\pi^{D-1}q^{D}}\left[-\frac{1}{2}\Theta(|\,\mathbf{q}+\mathbf{k}_{1}|-1)\Theta(|\,\mathbf{q}-\mathbf{k}_{2}|-1)\int\frac{d^{D}q^{\prime}}{|\mathbf{q}-\mathbf{q}^{\prime}|^{D-1}}{\cal
F}(\mathbf{k_{1}},\mathbf{k}_{2};\mathbf{q}^{\prime})\right.$
$\displaystyle+\left.\int
d^{D}k_{3}\Theta(1-k_{3})\left((\frac{1}{|\mathbf{q}+\mathbf{k}_{1}+\mathbf{k}_{3}|^{D-1}}-\frac{1}{|\mathbf{k}_{1}+\mathbf{k}_{3}|^{D-1}}){\cal
F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})\right.\right.$
$\displaystyle+\left.\left.\frac{\Theta(|\mathbf{q}+\mathbf{k}_{1}|-1)}{|\mathbf{k}_{1}-\mathbf{k}_{3}|^{D-1}}{\cal
F}(\mathbf{k}_{3},\mathbf{k}_{2};\mathbf{q})\right)\right.$
$\displaystyle+\left.\int
d^{D}k_{3}\Theta(1-k_{3})\left((\frac{1}{|\mathbf{q}-\mathbf{k}_{2}+\mathbf{k}_{3}|^{D-1}}-\frac{1}{|\mathbf{k}_{2}-\mathbf{k}_{3}|^{D-1}}){\cal
F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})\right.\right.$
$\displaystyle+\left.\left.\frac{\Theta(|\mathbf{q}-\mathbf{k}_{2}|-1)}{|\mathbf{k}_{2}+\mathbf{k}_{3}|^{D-1}}{\cal
F}(\mathbf{k}_{1},\mathbf{k}_{3};\mathbf{q})\right)\right].\quad$ (11)

Here $\hat{\mathbf{e}}_{\mathbf{q}}=\mathbf{q}/q$, $\Theta$ is the unit step
function, $\alpha_{2}=1/\sqrt{2}$, $\alpha_{3}=\sqrt[3]{4/9\pi}$, and the
Seitz radius $r_{s}$ is denoted by
$R_{2}=(\alpha_{2}K_{2}(\hbar^{2}/2m))^{-1}$ in two dimensions and
$R_{3}=(\alpha_{3}K_{3}(\hbar^{2}/2m))^{-1}$ in three dimensions.

As in Eq.(8) the direct terms on the left in Eq.(11) are large in comparison
with the homogeneous terms on the right. Neglecting these latter terms gives a
useful first approximation ${\cal F}\approx{\cal F}^{(1)}$ which can be found
by setting just the initial terms on left and right in Eq.(11) equal, so that

${\cal
F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})\approx\;\approx-\frac{\alpha_{D}R_{D}}{2(\pi
q)^{D-1}}\frac{\Theta(|\,\mathbf{q}+\mathbf{k}_{1}|-1)\Theta(|\,\mathbf{q}-\mathbf{k}_{2}|-1)}{(\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{q}+q^{2}}.$

This can be written as

$\displaystyle{\cal
F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})\approx\;\approx-\frac{\alpha_{D}R_{D}}{2\pi^{D-1}q^{D}}\Theta(|\,\mathbf{q}+\mathbf{k}_{1}|-1)\Theta(|\,\mathbf{q}-\mathbf{k}_{2}|-1)$
$\displaystyle\times\mathsf{Re}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi\imath}\frac{1}{-\nu+\mathbf{k}_{1}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}\,\frac{1}{\nu-\mathbf{k}_{2}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}\,{\cal
G}(\nu,q),$ (12)

where $\eta$ is a positive infinitesimal and where ${\cal G}(\nu,q)=1$. Next
substituting Eq.(12), now with ${\cal G}(\nu,q)$ undetermined, into the left
hand member of Eq.(11), where the coefficient of ${\cal
F}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})$ has been written as
$(-\nu+\mathbf{k}_{1}\cdot\mathbf{e}_{\mathbf{q}}+q/2)+(\nu-\mathbf{k}_{2}\cdot\mathbf{e}_{\mathbf{q}}+q/2)$,
and equating the result to zero leads to

$\mathsf{Re}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi\imath}\left(\frac{1}{-\nu+\mathbf{k}_{1}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}{\cal
D}(-\nu,q)+\frac{1}{\nu-\mathbf{k}_{2}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}{\cal
D}(\nu,q)\right){\cal G}(\nu,q)=1$ (13)

for the determination of ${\cal G}$. Here it is convenient to put

${\cal D}(\nu,q)=1+\frac{R_{D}}{q^{D-1}}{\cal U}(\nu,q),$ (14)

with ${\cal U}(\nu,q)$ given by

${\cal U}(\nu,q)=\frac{\alpha_{D}}{\pi^{D-1}q}\int
d^{D}k_{3}\frac{\Theta(1-k_{3})\Theta(|\mathbf{q}+\mathbf{k}_{3}|-1)}{-\nu+\mathbf{k}_{3}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}.$
(15)

The locations of the poles in the integrands of each of the terms on the left
of Eq.(13) can be restricted to just one of the half $\nu$-planes by choosing

${\cal G}(\nu,q)=({\cal D}(\nu,q){\cal D}(-\nu,q))^{-1},$ (16)

so that, on account of its symmetry, Eq.(13) will be satisfied provided that

$\mathsf{Re}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi\imath}\,\frac{1}{-\nu+\mathbf{k}_{1}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}\,\frac{1}{{\cal
D}(\nu,q)}=\frac{1}{2}.$ (17)

The validity of Eq.(17) rests on assumptions that the reciprocal of ${\cal
D}(\nu,q)$ can be expanded in powers of $R_{D}$ and that, in each term of the
resulting series, the integration over $\nu$ can be performed first. Then the
first term of the series integrates to one half while the remaining terms
vanish since their integrands decrease at least as rapidly as $|\nu|^{-2}$,
for large $|\nu|$, and are free of poles in the upper half plane.

Combining Eqs.(12) and (16) gives

$\displaystyle{\cal
F}^{(1)}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})=-\frac{\alpha_{D}R_{D}}{2\pi^{D-1}q^{D}}\Theta(|\,\mathbf{q}+\mathbf{k}_{1}|-1)\Theta(|\,\mathbf{q}-\mathbf{k}_{2}|-1)$
$\displaystyle\times\mathsf{Re}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi\imath}\frac{1}{-\nu+\mathbf{k}_{1}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}\,\frac{1}{\nu-\mathbf{k}_{2}\cdot\mathbf{e}_{\mathbf{q}}+q/2-\imath\eta}\frac{1}{{\cal
D}(\nu,q){\cal D}(-\nu,q)}.$ (18)

Replacing ${\cal F}$ on the right in Eq.(11) by ${\cal F}^{(1)}$ and on the
left by ${\cal F}^{(1)}+{\cal F}^{(2)}$, determines a second order
approximation which is discussed below in connection with $R_{3}\,\log\,R_{3}$
term of the correlation energy.

## Correlation energy

The correlation energy $\epsilon_{D}$, again in both two and three dimensions,
is obtained from Eq.(6) upon subtracting the HF contributions and dividing by
$N$. Similar to the relative importance of terms in Eq.(8), we expect that the
direct terms of $\epsilon_{D}$ will give contributions to the correlation
energy of all orders, while terms diminished by exchange, or in the case of
the last double sum of Eq.(6), by the implicit presence of an extra pair of
projection operators, will give contributions to orders $R_{D}\,\log\,R_{D}$
and higher. Identifying these contributions as $\epsilon_{D}^{(1)}$ and
$\epsilon_{D}^{(2)}$ respectively, and expressing the results in terms of the
transformed two electron functions Eq.(6) givesN (3)

$\displaystyle\epsilon_{D}^{(1)}({\cal
F})=\frac{D}{2^{D-2}\pi(\alpha_{D}R_{D})^{2}}\int
d^{D}k_{1}\Theta(1-k_{1})\int d^{D}k_{2}\Theta(1-k_{2})\int
d^{D}q\,((\mathbf{k}_{1}+\mathbf{k}_{2})\cdot\mathbf{q}+q^{2})$
$\displaystyle\times\left({\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q})^{2}-\frac{1}{2}{\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q}){\cal
F}(-\mathbf{k}_{2},\mathbf{k}_{1};\mathbf{q}+\mathbf{k}_{1}+\mathbf{k}_{2})\right)$
$\displaystyle+\frac{D}{2^{D-2}\pi^{D}\alpha_{D}R_{D}}\int
d^{D}k_{1}\Theta(1-k_{1})\int d^{D}k_{2}\Theta(1-k_{2})$
$\displaystyle\times\int d^{D}q\left(\frac{{\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q})}{q^{D-1}}-\frac{1}{2}\frac{{\cal
F}(-\mathbf{k}_{2},\mathbf{k}_{1};\mathbf{q}+\mathbf{k}_{1}+\mathbf{k}_{2})}{|\mathbf{q}+\mathbf{k}_{1}+\mathbf{k}_{2}|^{D-1}}\right.$
$\displaystyle+\left.\frac{2}{q^{D-1}}\int d^{D}k_{3}\Theta(1-k_{3}){\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q}){\cal
F}(\mathbf{k}_{1},\mathbf{k}_{3};\mathbf{q})\right),\qquad$ (19)

and

$\displaystyle\epsilon_{D}^{(2)}({\cal
F})=\frac{D}{2^{D-2}\pi^{D}\alpha_{D}R_{D}}\int d^{D}k_{1}\Theta(1-k_{1})\int
d^{D}k_{2}\Theta(1-k_{2})\int d^{D}q_{1}{\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q}_{1})$
$\displaystyle\times\left[-\int
d^{D}k_{3}\Theta(1-k_{3})\left((\frac{1}{|\mathbf{q}_{1}+\mathbf{k}_{1}+\mathbf{k}_{3}|^{D-1}}-\frac{1}{|\mathbf{k}_{1}+\mathbf{k}_{3}|^{D-1}}){\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q}_{1})\right.\right.$
$\displaystyle\left.\left.+\frac{{\cal
F}(\mathbf{k}_{2},-\mathbf{k}_{3};\mathbf{q}_{1})}{|\mathbf{k}_{1}-\mathbf{k}_{3}|^{D-1}}\right)+\frac{1}{2}\int
d^{D}q_{2}\frac{{\cal
F}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q}_{2})}{|\mathbf{q}_{1}-\mathbf{q}_{2}|^{D-1}}\right],\quad$
(20)

where the units are Ry rather than ergs. Note that the $\lambda(j_{1})$
dependent reduction of the HF kinetic energy term in Eq.(6) contributes to the
kinetic energy in Eq.(19).

Substituting Eq.(18) in (19), and using the assumptions which underlie Eq.(17)
to evaluate one of the two contour integrations in the two direct bilinear
terms of Eq.(19), gives

$\displaystyle\epsilon_{D}^{(1)}({\cal
F}^{(1)})=\frac{D}{2\pi\alpha_{D}^{2}}\int_{0}^{\infty}d\nu\int_{0}^{\infty}\frac{dq}{q^{D-2}}\mathsf{Im}\left(\frac{{\cal
U}(\nu,q){\cal U}(-\nu,q)}{|{\cal D}(\nu,q)|^{2}{\cal
D}(-\nu,q)^{2}}-2\frac{{\cal U}(\nu,q){\cal U}(-\nu,q)}{{\cal D}(\nu,q){\cal
D}(-\nu,q)}\right.$ $\displaystyle\left.+\frac{R_{D}}{q^{D-1}}\frac{{\cal
U}(\nu,q){\cal U}(-\nu,q)^{2}}{|{\cal D}(\nu,q)|^{2}{\cal
D}(-\nu,q)^{2}}\right)+\left\\{\begin{array}[]{c}0.2039\\\
0.0484\end{array}\right\\}$ (23)

which is equivalent to

$\epsilon_{D}^{(1)}({\cal
F}^{(1)})=-\frac{D}{2\pi\alpha^{2}}\int_{0}^{\infty}d\nu\int_{0}^{\infty}\frac{dq}{q^{D-2}}\mathsf{Im}\left(\frac{{\cal
U}(\nu,q){\cal U}(-\nu,q)}{{\cal D}(\nu,q){\cal
D}(-\nu,q)}\right)+\left\\{\begin{array}[]{c}0.2039\\\
0.0484\end{array}\right\\},$ (24)

where we have used the fact that the imaginary part of the function ${\cal
U}(\nu,q)$, evaluated in Appendix A, vanishes for $\nu<0$.

The additive constants (second order Rayleigh-Schrodinger exchange energies)
in Eqs.(21) and (22) are the sum of the kinetic and potential energy exchange
terms of Eq.(19) in the limit $R_{D}\rightarrow 0$. These were evaluated in
cylindrical coordinates in two dimensions (or Cartesian coordinates in three
dimensions) with symmetry axis $z$ along $\mathbf{q}$ and with the common
origin of the wavenumber vectors at $z=q/2$. A 45 degree rotation of axes
allows reductions to three-fold quadratures in both cases with results in
agreement with those of previous authors.Onsager et al. (1966); Rajagopal and
Kimball (1977)

In two dimensions, development of Eq.(22) in partial fractions now yields a
constant, $-0.4039$ Ry, plus terms proportional to $R_{2}\,\log R_{2}$ and to
$R_{2}$. Numerical fitting at 100 points
$R_{2}^{(k)}=10^{0.05k-8};\,k=1,\cdots 100$, so that $R_{2}<10^{-3}$ where the
contribution from order $R_{2}^{2}$ terms is expected to be small, gives

$\epsilon_{2}^{(1)}({\cal F}^{(1)})=-0.4039-0.50\,R_{2}\,\log
R_{2}-0.87\,R_{2},$ (25)

with an RMS deviation of $1.1\,10^{-5}$ from computed values. The result here
is in agreement with that found by Rajagopal and KimballRajagopal and Kimball
(1977) with the exception of the $R_{2}\,\log\,R_{2}$ term which is
approximately three times larger.

In three dimensions the presence of a logarithmic divergence with small
$R_{3}$ makes direct fitting unfeasible. Instead means must be found to
subtract the logarithmic dependence prior to evaluation. One way to do this is
to integrate in Eq.(21) with respect to only its explicit dependence on the
momentum $q$. This renders the integrand in Eq.(22) in the form of a partial
derivative with respect to $q$, of a quantity proportionate to

$G(\nu,q)\equiv\frac{{\cal U}(\nu,q){\cal U}(-\nu,q)}{{\cal U}(-\nu,q)-{\cal
U}(\nu,q)}({\cal U}(\nu,q)\log(q^{2}+R_{3}\,{\cal U}(\nu,q))-{\cal
U}(-\nu,q)\log(q^{2}+R_{3}\,{\cal U}(-\nu,q))).$ (26)

Expressing this partial derivative in terms of its related total derivative
isolates the divergence at small $R_{3}$.

A complication of this method is that the imaginary part of the integrand in
Eq.(22) picks up contributions from a new area in the first quadrant of the
$\nu-q$ plane where the argument of $\log(q^{2}+R_{3}\,{\cal U}(\nu,q))$ in
Eq.(24) becomes negative. This happens on and to the right of the line
$q=2(\nu-1)$ and below the curve $q=q_{c}(\nu)$, where
$q_{c}(\nu)^{2}+R_{3}\,{\cal U}(\nu,q_{c}(\nu))=0$. Here $q_{c}(\nu)$
originates at a point $0<\nu<1$ on the $\nu$ axis, reaches a maximum for
$\nu\stackrel{{{}_{>}}}{{{}_{\approx}}}1$, and, with increasing $\nu$, turns
downward, crossing the line $q=2(\nu-1)$ at the point
$(\nu_{c},2(\nu_{c}-1))$, where $\nu_{c}\geq 1$ for $R_{3}\geq 0$. Finally
$q_{c}(\nu)$ approaches zero for large $\nu$.

Summing contributions in Eq.(22) from the area between the lines $q=2(\nu\pm
1)$, where we have $\mathsf{Im}\,{\cal U}(\nu,q)\neq 0$, and then from the
area where $\mathsf{Im}\,\log(q^{2}+R_{3}\,{\cal U}(\nu,q))=\pi$, gives

$\displaystyle\epsilon_{3}^{(1)}=\frac{3}{\pi\alpha^{2}}\left[\left(\int_{0}^{1}d\nu\int_{0}^{2(\nu+1)}dq+\int_{1}^{\infty}d\nu\int_{2(\nu-1)}^{2(\nu+1)}dq\right)\right.\
$ $\displaystyle\left.\times\left(\frac{d}{dq}-\frac{\partial\,{\cal
U}(\nu,q)}{\partial q}\frac{\partial}{\partial\,{\cal
U}(\nu,q)}-\frac{\partial\,{\cal U}(-\nu,q)}{\partial
q}\frac{\partial}{\partial\,{\cal
U}(-\nu,q)}\right)\mathsf{Im}\,G(\nu,q)\right.$
$\displaystyle+\left.\int_{1}^{\infty}d\nu\int_{0}^{2(\nu-1)}dq\frac{{\cal
U}(\nu,q)^{2}{\cal U}(-\nu,q)}{{\cal U}(-\nu,q)-{\cal
U}(\nu,q)}\frac{\partial}{\partial q}\pi\Theta(-(q^{2}+R_{3}\,{\cal
U}(\nu,q)))\right].\quad$ (27)

Here upper limits of the $q$ integrations over the total derivative do not
contribute, while the lower limit of the first of these integrations is
proportionate to $\log\,R_{3}$. The lower limit of the second $q$ integration
over the total derivative leaves an integration over $\nu$ running from one to
infinity with $\mathsf{Im}\,G(\nu,q)$ evaluated on the right most bounding
line $q=2(\nu-1)$. However only along the initial portion of this bounding
line, where $\nu<\nu_{c}$ so that $q^{2}+R_{3}\,{\cal U}(\nu,q)<0$, is
$\mathsf{Im}\,G(\nu,q)$ non-zero. The sum of the double integrals over the
partial derivatives of $\mathsf{Im}\,G(\nu,q)$ is bounded.

In the final double integral in Eq.(25) the derivative of the step function is
proportional to $\delta(q-q_{c}(\nu))$ so the integral over $q$ vanishes where
$2(\nu-1)<q_{c}(\nu)$, that is where $\nu<\nu_{c}$. This leaves for the last
line of Eq.(25)

$-\frac{3}{4\alpha_{3}^{2}}\int^{\infty}_{\nu_{c}}d\nu\frac{{\cal
U}(\nu,q_{c}(\nu))^{2}{\cal U}(-\nu,q_{c}(\nu))}{{\cal
U}(-\nu,q_{c}(\nu))-{\cal
U}(\nu,q_{c}(\nu))}\,\frac{1}{1+(R_{3}/2q_{c}(\nu))(\partial/\partial
q_{c}(\nu)){\cal U}(\nu,q_{c}(\nu))},$

where the final fraction approaches unity as $R_{3}$ approaches zero since the
$q$ derivative of $\mathsf{Re}\;{\cal U}(\nu,q)$ is itself proportional to
$q$.

With the divergent term, $0.0622\,\log\,R_{3}$ Ry, set aside, the remainder
was examined numerically at 60 points $R_{3}^{(k)}=10^{0.05k-6};\,k=1,\cdots
60$, so that $0<R_{3}<0.001$, with the result, inclusive of the second order
exchange energy $0.0484$ Ry,

$\epsilon_{3}^{(1)}({\cal F}^{(1)})=0.0622\,\log R_{3}-0.082-0.015\,R_{3},$
(28)

with an RMS deviation of $1.2\,10^{-5}$ from computed values. No $\,R_{3}\log
R_{3}$ dependence was found. We note that if the direct potential energy term
quadratic in the two electron functions in Eq.(8), (this is responsible for
the term proportionate to $R_{3}$ on the right in Eq.(21)), were not included
in the subsequent calculations, the constant term of Eq.(26) becomes $-0.093$
Ry, approximately the result found from perturbation theory.Gell-Mann and
Breuckner (1957)

## Order $R_{D}\,\log\,R_{D}$

Logarithmic contributions are determined in the transformed domain by behavior
at small $q$, where the functions ${\cal U}$ in the denominator of Eq.(18)
approach finite positive limits. Accordingly we may now simplify ${\cal
F}^{(1)}$ by replacing each of the functions ${\cal D}$ in the denominator of
Eq.(18) by $(1+R_{D}/q^{D-1})$ and performing the $\nu$-integration to get

${\cal
F}^{(1)}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})\approx-\frac{\alpha_{D}R_{D}}{2\pi^{D-1}q^{D}}\frac{\Theta(|\,\mathbf{q}+\mathbf{k}_{1}|-1)\Theta(|\,\mathbf{q}-\mathbf{k}_{2}|-1)}{(1+R_{D}/q^{D-1})^{2}((\mathbf{k}_{1}-\mathbf{k}_{2})\cdot\mathbf{\hat{e}}_{\mathbf{q}}+q)}.$
(29)

In terms of this simplified version of ${\cal F}^{(1)}$, we have

$\displaystyle\epsilon_{D}^{(2)}({\cal
F}^{(1)})=\frac{D\alpha_{D}R_{D}}{2^{D}\pi^{3D-2}}\left\\{-\int\frac{d^{D}q}{q^{2D}}\frac{1}{(1+R_{D}/q^{D-1})^{4}}\int
d^{D}k_{1}\Theta(1-k_{1})\Theta(|\mathbf{q}+\mathbf{k}_{1}|-1)\right.$
$\displaystyle\left.\times\int
d^{D}k_{2}\Theta(1-k_{2})\frac{\Theta(|\mathbf{q}+\mathbf{k}_{2}|-1)}{((\mathbf{k}_{1}+\mathbf{k}_{2})\cdot\hat{\mathbf{e}}_{\mathbf{q}}+q)}\left[\frac{1}{(\mathbf{k}_{1}+\mathbf{k}_{2})\cdot\hat{\mathbf{e}}_{\mathbf{q}}+q}\int
d^{D}k_{3}\,\Theta(1-k_{3})\right.\right.$
$\displaystyle\times\left.\left.\left(\frac{1}{|\mathbf{q}+\mathbf{k}_{1}+\mathbf{k}_{3}|^{D-1}}-\frac{1}{|\mathbf{k}_{1}+\mathbf{k}_{3}|^{D-1}}\right)+\int
d^{D}k_{3}\frac{\Theta(1-k_{3})\,\Theta(|\mathbf{q}+\mathbf{k}_{3}|-1)}{((\mathbf{k}_{2}+\mathbf{k}_{3})\cdot\hat{\mathbf{e}}_{\mathbf{q}}+q)\;|\mathbf{k}_{1}-\mathbf{k}_{3}|^{D-1}}\right]\right.$
$\displaystyle+\left.\frac{1}{2}\int\frac{d^{D}q_{1}}{q^{D}_{1}}\frac{1}{(1+R_{D}/q_{1}^{D-1})^{2}}\int\frac{d^{D}q_{2}}{q^{D}_{2}}\frac{1}{(1+R_{D}/q^{D-1}_{2})^{2}}\frac{1}{|\mathbf{q}_{1}-\mathbf{q}_{2}|^{D-1}}\int
d^{D}k_{1}d^{D}k_{2}\Theta(1-k_{1})\right.$
$\displaystyle\left.\times\Theta(1-k_{2})\frac{\Theta(|\mathbf{q}_{1}+\mathbf{k}_{1}|-1)\Theta(|\mathbf{q}_{1}+\mathbf{k}_{2}|-1)}{(\mathbf{k}_{1}+\mathbf{k}_{2})\cdot\mathbf{e}_{\mathbf{q}_{1}}+q_{1}}\,\frac{\Theta(|\mathbf{q}_{2}+\mathbf{k}_{1}|-1)\Theta(|\mathbf{q}_{2}+\mathbf{k}_{2}|-1)}{(\mathbf{k}_{1}+\mathbf{k}_{2})\cdot\mathbf{e}_{\mathbf{q}_{2}}+q_{2}}\right\\}.\qquad$
(30)

In two dimensions the curly-bracketed term of Eq.(28) remains finite in the
limit $R_{2}\rightarrow 0$ so that there is no order $R_{2}\,\log R_{2}$
contribution. In three dimensions the $q$-integrations over the square
bracketed integrals in Eq.(28), considered in isolation, each diverge
logarithmically. However, as shown in Appendix B, these divergences cancel (as
they must) with the result that

$\epsilon_{3}^{(2)}=0.059\,R_{3}\,\log\,R_{3}.$ (31)

A second potential source of order $R_{D}\,\log\,R_{D}$ contributions is

$\displaystyle\epsilon_{D}^{(1)}({\cal F}^{(2)}+{\cal
F}^{(1)})-\epsilon_{D}^{(1)}({\cal
F}^{(1)})=\frac{D}{(2\pi)^{D-2}(\alpha_{D}R_{D})^{2}}\int
d^{D}k_{1}\Theta(1-k_{1})\int d^{D}k_{2}\Theta(1-k_{2})$
$\displaystyle\times\int
d^{D}q((\mathbf{k}_{1}+\mathbf{k}_{2})\cdot\mathbf{q}+q^{2}){\cal
F}^{(1)}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q}){\cal
F}^{(2)}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q})$
$\displaystyle+\frac{D}{\pi^{D}\alpha_{D}R_{D}}\int
d^{D}k_{1}\Theta(1-k_{1})\int
d^{D}k_{2}\Theta(1-k_{2})\int\frac{d^{D}q}{q^{D-1}}{\cal
F}^{(2)}(\mathbf{k}_{1},-\mathbf{k}_{2};\mathbf{q}),\qquad$ (32)

where the small exchange and final terms of Eq.(19) have been dropped. Here we
expect ${\cal F}^{(2)}$ to be well defined as a result of cancelations of the
small $q$ divergences on the right in Eq.(11) analogous to the cancelation
observed in Eq.(28).

Expressing the left hand member of Eq.(11) in the form of an integral operator
$\hat{B}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})$ acting on ${\cal F}$,
then, analogous with the classical case,Montgomery and Tidman (1964) the
inverse operator

$\hat{B}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})^{-1}=\mathsf{Re}\int_{-\infty}^{\infty}\frac{d\nu}{2\pi\imath}\hat{V}(\mathbf{k}_{1},\mathbf{q},\nu)\hat{V}(\mathbf{k}_{2},-\mathbf{q},-\nu),$
(33)

where

$\displaystyle\hat{V}(\mathbf{k}_{1},\mathbf{q},\nu)=\frac{1}{-\nu+\mathbf{k}\cdot\mathbf{\hat{e}}_{\mathbf{q}}+q/2-\imath\eta}$
$\displaystyle\times\left(1-\frac{\alpha_{D}R_{D}}{\pi^{D-1}q^{D}}\frac{\Theta(|\mathbf{q}+\mathbf{k}|-1)}{{\cal
D}(\nu,q)}\int
d\mathbf{k}^{\prime}\frac{\Theta(1-k^{\prime})}{-\nu+\mathbf{k}^{\prime}\cdot\mathbf{\hat{e}}_{\mathbf{q}}+q/2-\imath\eta}\hat{\bigcirc}(\mathbf{k}\rightarrow\mathbf{k}^{\prime})\right),$

allows one to obtain an explicit expression for ${\cal
F}^{(2)}(\mathbf{k}_{1},\mathbf{k}_{2};\mathbf{q})$. This, however, turns out
not to be needed since substitution from Eq.(27) shows that the order
$R_{D}\,\log\,R_{D}$ terms on the right in Eq.(30) sum to zero in both two and
three dimensions.

Summing Eqs.(26) and (29) gives finally

$\epsilon_{3}=0.0622\,\log\,R_{3}-0.082+0.059\,R_{3}\,\log
R_{3}-0.015\,R_{3}.$ (34)

## Order $R_{3}$

An accounting of contributions to order $R_{3}$ of $\epsilon_{3}$ requires use
of the initial two terms of the GSJ series in Eq.(3) when deriving the one
electron Euler equation with the consequent inclusion on the right in Eq.(7)
of the additional termsPorter (1965b)

$\displaystyle+\sum_{j_{3}}\int
d^{3}r_{3}U^{\dagger}(j_{3};3)\phi(1.3)F(j_{1},j_{2};1,2)$
$\displaystyle-\sum_{j_{3},\,j_{4}}\int
d^{3}r_{3}d^{3}r_{4}U^{\dagger}(j_{3};3)U^{\dagger}(j_{4};4)\phi(3,4)U(j_{1};3)F(j_{4},j_{3};1,4).$
(35)

(It is still unnecessary to include weighting functions.) Like the HF exchange
term, these latter two terms, when imported into the two electron equation via
the Lagrange multipliers, diverge with small $q$. As happened with the HF
exchange contributions, the divergent portions of these contributions are
canceled by four additional bilinear terms on the right in Eq.(8),

$\displaystyle\sum_{j_{3},\,j_{4}}\int
d^{3}r_{3}d^{3}r_{4}F^{\dagger}(j_{3},j_{4};3,4)\phi(3,4)U(j_{4};4)\left(U(j_{3};1)F(j_{2},j_{1};2,3)+U(j_{3};2)F(j_{1},j_{2};1,3)\right)$
$\displaystyle-\sum_{j_{3},\,j_{4},\,j_{5}}\int
d^{3}r_{3}d^{3}r_{4}d^{3}r_{5}F^{\dagger}(j_{3},j_{4};3,4)U^{\dagger}(j_{5};5)\phi(4,5)U(j_{3};5)U(j_{4};4)$
$\displaystyle\times\left(U(j_{5};1)F(j_{2},j_{1};2,3)+U(j_{5};2)F(j_{1},j_{2};1,3)\right),\qquad$
(36)

which arise when one extracts, without approximation, a relation for
$F(j_{1},j_{2};1,2)$ alone from the original result obtained from variations,
with respect to $f^{\dagger}(j_{1},j_{2};1,2)$, as discussed following Eq.(7).

Numerical evaluation of all order $R_{3}$ contributions is problematic since
the simplification in Eq.(27) cannot be used. In the transformed domain, order
$R_{3}$ integrands contain two or three additional powers of $q$ in comparison
to those of order $R_{3}\,\log\,R_{3}$ terms so that short range effects
dominate. The alternating signs of these integrals, and the fact that each
contains at least one exchange term, leaves $-0.015$ Ry, Eq.(26), from the
direct terms of Eq.(8), as a reasonable guess for the order of magnitude of
the order $R_{3}$ contribution to $\epsilon_{3}$. Thus the correlation
energy’s order $R_{3}$ term is found to be roughly equal, and the order
$R_{3}\,\log\,R_{3}$ term some three times larger, than the respective results
obtained from perturbation theory.Carr and Maradudin (1964)

*

## Appendix A

In the following inequalities are arranged with the expectation that the $\nu$
integrations required by Eq.(25) be performed first.

In two dimensions it is convenient to define the auxiliary functions

$\displaystyle\begin{array}[]{l}L(z,\nu,q)=\sqrt{1-(\nu-q/2)^{2}}\,\log\left(\frac{f(z,\nu,q)-\sqrt{(1-(\nu-q/2)^{2})(1-(z-q/2)^{2})}}{z-\nu}\right)+g(z,\nu,q),\\\
A(z,\nu,q)=\sqrt{(\nu-q/2)^{2}-1}\,\arcsin\left(\frac{f(z,\nu,q)}{z-\nu}\right)+g(z,\nu,q),\\\
C(z,\nu,q)=\sqrt{(\nu-q/2)^{2}-1}\,\arccos\left(\frac{f(z,\nu,q)}{z-\nu}\right)+g(z,\nu,q),\end{array}$

where

$\displaystyle\begin{array}[]{l}f(z,\nu,q)=1-q^{2}/4+\nu q/2-\nu z+qz/2,\\\
g(z,\nu,q)=(\nu-q/2)\arcsin(q/2-z)+\sqrt{1-(z-q/2)^{2}}.\end{array}$

Next we define

$\displaystyle\begin{array}[]{l}Llt2(\nu,q)=L(1+q/2,\nu,q)-L(0,\nu,q),\\\
Lgt2(\nu,q)=L(q/2+1,\nu,q)-L(q/2-1,\nu,q),\\\
Alt2(\nu,q)=A(1+q/2,\nu,q)-A(0,\nu,q),\\\
Agt2(\nu,q)=A(q/2+1,\nu,q)-A(q/2-1,\nu,q),\\\
Clt2(\nu,q)=C(1+q/2,\nu,q)-C(0,\nu,q),\\\
Cgt2(\nu,q)=C(q/2+1,\nu,q)-C(q/2-1,\nu,q).\end{array}$

In terms of this last set of definitions

$\displaystyle\mathsf{Re}\,{\cal U}(q,\nu)=\frac{2\alpha_{2}}{\pi
q}\mathsf{Re}\left\\{\begin{array}[]{l}Llt2(\nu,q)-Llt2(\nu,-q);\quad-1+q/2<\nu<1+q/2\text{
and }0<q<2,\\\ Slt2(\nu,q)-Slt2(\nu,-q);\quad 1+q/2<\nu\text{ and }0<q<2,\\\
Clt2(\nu,q)-Clt2(\nu,-q);\quad\nu<-1+q/2\text{ and }0<q<2,\\\
Lgt2(\nu,q);\quad-1+q/2<\nu<1+q/2\text{ and }q>2,\\\ Sgt2(\nu,q);\quad
1+q/2<\nu\text{ and }q>2,\\\ Cgt2(\nu,q);\quad\nu<-1+q/2\text{ and }q>2,\\\
\end{array}\right.$ (A5)

and

$\displaystyle\mathsf{Im}\,{\cal
U}(q,\nu)={2\alpha_{2}\left\\{\begin{array}[]{l}\vspace{.1in}\frac{2\nu}{\sqrt{1-(\nu-q/2)^{2}}+\sqrt{1-(\nu+q/2)^{2}}};\quad
0<\nu<1-q/2\text{ and }0<q<2,\\\
\frac{\sqrt{1-(\nu-q/2)^{2}}}{q};\quad(1-q/2<\nu<1+q/2\text{ and }0<q<2)\\\
\qquad\text{ or }(q/2-1<\nu<q/2+1\text{ and }q<2),\\\
0;\quad\text{otherwise}.\end{array}\right.}$ (A8) (A9)

In three dimensions

$\displaystyle\mathsf{Re}\,{\cal
U}(q,\nu)=\frac{\alpha_{3}}{\pi}\mathsf{Re}{\left\\{\begin{array}[]{l}1-\nu-2\nu\log(\nu)+\frac{1-(\nu-q/2)^{2}}{q}\log(1-\nu+q/2)\\\
\qquad-\frac{1-(\nu+q/2)^{2}}{q}\log(1-\nu-q/2);\quad 0<q<2,\\\
1-\frac{2\nu}{q}+\frac{1-(\nu-q/2)^{2}}{q}\log\left(\frac{1-\nu+q/2}{1+\nu-q/2}\right);\quad
q>2,\end{array}\right.}$ (A11) (A12)

and

$\displaystyle\mathsf{Im}\,{\cal
U}(q,\nu)={\alpha_{3}\left\\{\begin{array}[]{lcl}2\nu;\quad 0<\nu<1-q/2\text{
and }0<q<2,\\\ \frac{(1-(\nu-q/2)^{2})}{q};\quad(1-q/2<\nu<1+q/2\text{ and
}0<q<2)\\\ \qquad\text{ or }(q/2-1<\nu<q/2+1\text{ and }q>2),\\\
0;\quad\text{otherwise}.\end{array}\right.}$ (A15) (A16)

*

## Appendix B

We adopt, for the $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$ integrations in the
first term within the square brackets of Eq.(28), spherical coordinates with
$z$ axis parallel to $\mathbf{q}$. For small $q$ the step functions require
that $1>k_{1}>1-q\cos\,\theta_{1}$, where $\cos\,\theta_{1}>0$, with a
corresponding relation for $k_{2}$. Carrying out the trivial integrations over
$\mathbf{k}_{3}$, $\phi_{1}$ and $\phi_{2}$ gives

$\displaystyle-\frac{3\alpha
R_{3}}{8\pi^{7}}(2\pi)^{3}\int_{0}\frac{d\mathbf{q}}{q^{6}}\int_{0}^{\pi/2}\sin\theta_{1}d\theta_{1}\int_{1-q\cos\,\theta_{1}}^{1}dk_{1}\,k_{1}^{2}\int_{0}^{\pi/2}\sin\theta_{2}d\theta_{2}\int_{1-q\cos\,\theta_{2}}^{1}dk_{2}\,k_{2}^{2}$
$\displaystyle\times\frac{1}{(1+R_{3}/q^{2})^{4}}\frac{1}{(k_{1}\cos\theta_{1}+k_{2}\cos\theta_{2}+q)^{2}}$
$\displaystyle\times\left(\frac{1-|\mathbf{k}_{1}+\mathbf{q}|^{2}}{2|\mathbf{k}_{1}+\mathbf{q}|}\log\left|\frac{1+|\mathbf{k}_{1}+\mathbf{q}|}{1-|\mathbf{k}_{1}+\mathbf{q}|}\right|-\frac{1-k_{1}^{2}}{2k_{1}}\log\left|\frac{1+k_{1}}{1-k_{1}}\right|\right).$
(B1)

Putting $|\mathbf{k}_{1}+\mathbf{q}|\approx k_{1}+q\cos\theta_{1}$ and
introducing dummy variables of integration $y_{i}$ such that
$k_{i}=1-qy_{i}\cos\theta_{i}$ for $i=1,2$ gives, again for small $q$,

$\displaystyle-\frac{12\alpha
R_{3}}{\pi^{3}}\int_{0}dq\frac{q^{6}}{(q^{2}+R_{3})^{4}}\int_{0}^{1}dz_{1}\int_{0}^{1}dy_{1}\,z_{1}\int_{0}^{1}dz_{2}\int_{0}^{1}dy_{2}\,z_{2}\frac{1}{(z_{1}+z_{2})^{2}}$
$\displaystyle\times\left(\frac{1-(1-qy_{1}z_{1}+qz_{1})^{2}}{2}\log\left|\frac{2}{1-(1-qy_{1}z_{1}+qz_{1})}\right|\right.$
$\displaystyle-\left.\frac{1-(1-qy_{1}z_{1})^{2}}{2}\log\left|\frac{2}{1-(1-qy_{1}z_{1})}\right|\right),$
(B2)

where we have denoted the cosines $\cos\theta_{i}$ by $z_{i}$ for $i=1,2$.

For the second term within the square brackets of Eq.(28) we adopt cylindrical
coordinates, with symmetry axis along $\mathbf{q}$, for the integrations over
$\mathbf{k}_{i}$, $i=1,2,3$. Here the origins of the cylindrical systems are
located at a point $-\mathbf{q}/2$ from the origins of the $\mathbf{k}_{i}$ on
the $z$-axis. Thus, if $\mathbf{k}_{i}$ locates a point
$(\rho_{i},\,z_{i},\,\phi_{i})$ then
$\mathbf{k}_{i}\cdot\hat{\mathbf{e}}_{\mathbf{q}}=z_{i}-q/2$. The small $q$
approximations of the integral operators $\int
d\mathbf{k}_{i}\Theta(1-k_{i})\,\Theta(|\mathbf{q}+\mathbf{k}_{i}|-1)$ now
appear as

$\int_{0}^{1}dz_{i}\int_{\sqrt{1-(z_{i}+q/2)^{2}}}^{\sqrt{1-(z_{i}-q/2)^{2}}}\rho_{i}d\rho_{i}\int_{0}^{2\pi}d\phi_{i};\;i=1,\,2,\,3.$

Because $\int_{0}^{2\pi}d\phi_{3}(\mathbf{k}_{1}-\mathbf{k}_{3})^{-2}$ is a
function of $\rho_{1}^{2}$ and $\rho_{3}^{2}$ only, it is convenient to define
$Y_{i}=\rho_{i}^{2}$ for $i=1,\,2,\,3$. Thus for small $q$ this second term
becomes, upon performing the $\phi$-integrations

$\displaystyle-\frac{3\alpha
R_{3}}{\pi^{3}}\int_{0}dq\frac{q^{4}}{(q^{2}+R_{3})^{4}}\int_{0}^{1}dz_{1}\int_{0}^{1}dz_{2}\int_{0}^{1}dz_{3}$
$\displaystyle\times\int_{1-(z_{1}+q/2)^{2}}^{1-(z_{1}-q/2)^{2}}dY_{1}\int_{1-(z_{2}+q/2)^{2}}^{1-(z_{2}-q/2)^{2}}dY_{2}\int_{1-(z_{3}+q/2)^{2}}^{1-(z_{3}-q/2)^{2}}dY_{3}$
$\displaystyle\times\frac{1}{\sqrt{\left(\frac{}{}Y_{1}+(z_{1}-q/2)^{2}+Y_{3}+(z_{3}-q/2)^{2}-2(z_{1}-q/2)(z_{3}-q/2)\right)^{2}-4Y_{1}Y_{3}}}.$
(B3)

A change of dummy variables from the $Y_{i}$ to $y_{i}$, where
$Y_{i}=1-z_{i}^{2}+z_{i}y_{i}q-q^{2}/4$ for $i=1,\,2,\,3$, allows the
interchange of the integrals over $z_{3}$ and $y_{3}$, while reducing the
argument of the square root in (B3) to a quadratic in $z_{3}$, with the
reciprocal of the square root being multiplied by
$z_{1}z_{2}z_{3}/((z_{1}+z_{2})(z_{2}+z_{3}))$, so that the $z_{3}$ integral
in (B3) can be evaluated analytically. These changes of variable give a
multiplicative factor of $q^{3}$ so that the lower limit of the initial
integration over $q$ yields (for at least the $q$-independent part of its
integrand) the expected the $R_{3}\,\log R_{3}$ dependence.

The $z_{3}$-integration reveals a divergent term proportional to $\log\,q$
which cancels the $\log\,q$ divergence found in (B2). The remaining terms form
a power series in $q$ whose constant term can be evaluated numerically with
the result $0.0172\,R_{3}\,\log\,R_{3}$ Ry.

The wave number integrations found in the last term of Eq.(28) are best
evaluated in spherical coordinates. Here we observe that, for small $q_{1}$
and $q_{2}$, the arguments of the four step functions present can be replaced
by $k_{i}-(1-q_{j}\cos\theta_{i\,j})$, with $i=1,\,2$ and $j=1,\,2$ as
appropriate. Here, for a non-zero result,
$\cos\theta_{i\,j}=\mathbf{k}_{i}\cdot\mathbf{q}_{j}/(k_{i}\,q_{j})$ must be
positive. We next introduce spherical coordinates with $z$-axis along
$\mathbf{k}_{1}$ and $x$-axis in the $\mathbf{k}_{1}-\mathbf{k}_{2}$ plane
such that the $x$-component of $\mathbf{k}_{2}$ is positive. The angle
$\theta_{21}$ between $\mathbf{k}_{2}$ and $\mathbf{q}_{1}$ can be expressed
in terms of the polar and azimuthal angles of $\mathbf{q}_{1}$ by means of the
relation
$\cos\theta_{21}=\cos\beta\,\cos\theta_{11}+\sin\beta\,\sin\theta_{11}\,\cos\phi_{1}$,
where $\beta$ is the angle separating $\mathbf{k}_{1}$ and $\mathbf{k}_{2}$. A
corresponding expression gives $\cos\theta_{22}$ in terms of $\beta$,
$\theta_{12}$ and $\phi_{2}$.

The final term of Eq.(28) now becomes

$\displaystyle\frac{3\alpha_{3}R_{3}}{4\pi^{5}}\int^{\pi}_{0}d\beta\,\sin\beta\int_{0}d\mathbf{q_{2}}\frac{q_{2}}{(q^{2}_{2}+R_{3})^{2}}\int_{0}d\mathbf{q}_{1}\frac{q_{1}}{(q^{2}_{1}+r-3)^{2}}\frac{1}{q_{1}^{2}-2\mathbf{q}_{1}\cdot\mathbf{q}_{2}+q_{2}^{2}}$
$\displaystyle\times\int^{1}_{0}dk_{1}\,k_{1}^{2}\int^{1}_{0}dk_{2}\,k_{2}^{2}\,\frac{\Theta(k_{1}-(1-q_{1}\cos\theta_{11}))\,\Theta(\cos\theta_{11})\,\Theta(k_{2}-(1-q_{1}\cos\theta_{21}))\,\Theta(\cos\theta_{21})}{k_{1}\cos\,\theta_{11}+k_{2}\cos\,\theta_{21}}$
$\displaystyle\times\frac{\Theta(k_{1}-(1-q_{2}\cos\theta_{12}))\,\Theta(\cos\theta_{12})\,\Theta(k_{2}-(1-q_{2}\cos\theta_{22}))\,\Theta(\cos\theta_{22})}{k_{1}\cos\,\theta_{12}+k_{2}\cos\,\theta_{22}},\qquad$
(B4)

where the zero subscript marking the $\mathbf{q}$ integrations indicates that
only the $q\rightarrow 0$ lower limit is expected to contribute.

The step functions in (B4) require $k_{1}$ and $k_{2}$ to differ from unity
only by terms of order $q$ so that, with the exception of their appearance in
the arguments of the step functions, $k_{1}$ and $k_{2}$ in (B4) may each be
replaced by unity. The resulting trivial integrals over $k_{1}$ and $k_{2}$
give

$\displaystyle\frac{3\alpha
R_{3}}{4\pi^{5}}\int^{\pi}_{0}d\beta\,\sin\beta\int_{0}d\mathbf{q}_{2}\frac{q_{2}}{(q_{2}^{2}+R_{3})^{2}}\int_{0}d\mathbf{q}_{1}\frac{q_{1}}{(q_{1}^{2}+R_{3})^{2}}\frac{1}{q_{1}^{2}-2\mathbf{q}_{1}\cdot\mathbf{q}_{2}+q_{2}^{2}}$
$\displaystyle\times\frac{\Theta(\cos\,\theta_{11})\,\Theta(\cos\,\theta_{12})\Theta(\cos\,\theta_{21})\,\Theta(\cos\,\theta_{22})}{(\cos\theta_{11}+\cos\theta_{21})(\cos\theta_{12}+\cos\theta_{22})}\cos\theta_{12}\;\Theta(q_{1}-q_{2}b_{1})$
$\displaystyle\times\left(q_{1}q_{2}\,\cos\theta_{21}\,(1-\Theta(q_{1}-q_{2}b_{2}))+q_{2}^{2}\cos\theta_{22}\,\Theta(q_{1}-q_{2}b_{2})\frac{}{}\right)+\cdots,\quad$
(B5)

where we have defined $b_{1}=\cos\theta_{12}/\cos\theta_{11}>0$ and
$b_{2}=\cos\theta_{22}/\cos\theta_{21}>0$, and where the trailing dots
indicate terms duplicating those explicitly given with the exception that the
roles of $q_{1}$ and $q_{2}$ are interchanged.

Limits defining the angular portions of the $\mathbf{q}$ integrations depend
in a complicated way on the angle $\beta$ separating $\mathbf{k}_{1}$ and
$\mathbf{k}_{2}$. In the case where the angle $\beta$ equals zero, the angular
portions of the $q$-integrations in (B5) range over the $2\pi$ solid angle
above a plane perpendicular to both wave vectors $\mathbf{k}_{1}$ and
$\mathbf{k}_{2}$ and tangent to the unit Fermi sphere at a point an
infinitesimal distance above their coincident tips. As $\beta$ increases from
zero, the two wave vectors separate so that the solid angle over which the $q$
integrations are to range is reduced to $2\pi-2\beta$. This is because the
half of the $\mathbf{k}_{1}$ tangent plane for which $x<0$ is now replaced by
an upturned half plane which has been rotated counterclockwise (as viewed from
a great distance out along the positive $y$-axis) through an angle $\beta$
about a fold line parallel to the $y$-axis. This guarantees that
$\mathbf{k}_{2}$ plus an infinitesimally small $\mathbf{q}_{2}$ will still lie
outside the unit sphere.

The extent of the azimuthal integration for $\mathbf{q}_{1}$ is curtailed when
$\beta+\theta_{11}$ is greater than $\pi/2$. In this case the circle traced by
the tip of the unit vector $\mathbf{q}_{1}/q_{1}$, with the polar angle
$\theta_{11}$ fixed, intersects the above up-turned plane at the azimuthal
positions $\phi_{1}=\pi\pm\xi_{1}$, where
$\xi_{1}=\arccos(\tan(\pi/2-\beta)/\tan\,\theta_{11})$, so that the
integration over the azimuthal angle of $\mathbf{q}_{1}$ is restricted to the
range $-(\pi-\xi_{1})<\phi_{1}<\pi-\xi_{1}$. A corresponding restriction
applies to $\mathbf{q}_{2}$ where
$\xi_{2}=\arccos(\tan(\pi/2-\beta)/\tan\,\theta_{12})$.

An integration by parts, $vdu=uv-udv$, with $u$ the indefinite integral of
either $q_{1}^{4}$ or $q_{1}^{5}$ times the reciprocal of
$(q_{1}^{2}+R_{3})^{2}(q_{1}^{2}-2q_{1}\,q_{2}\,\cos\,\psi+q_{2}^{2})$, where
$\cos\,\psi=\cos\,\theta_{11}\cos\,\theta_{12}+\sin\,\theta_{11}\sin\,\theta_{12}\,\cos(\phi_{1}-\phi_{2})$,
and with $v$ a step function, (or a product of two step functions) which
vanishes at the lower limit $q_{1}=0$, gives a $udv$ term which includes a
delta function: $\delta(q_{1}-q_{2}\,b_{1})$ or $\delta(q_{1}-q_{2}\,b_{2})$.
The subsequent integration over $q_{2}$ of the explicitly written out terms in
(B5) gives

$\displaystyle\frac{3\alpha
R_{3}\log\,R_{3}}{4\pi^{5}}\left\\{\int^{\pi/2}_{0}d\beta\sin\,\beta\left[\left(\int^{\pi/2}_{\pi/2-\beta}d\theta_{11}\sin\,\theta_{11}\int^{\pi-\xi_{1}}_{-(\pi-\xi_{1})}d\phi_{1}+\int^{\pi/2-\beta}_{0}d\theta_{11}\sin\,\theta_{11}\int^{\pi}_{-\pi}d\phi_{1}\right)\right.\right.$
$\displaystyle\times\left(\int^{\pi/2}_{\pi/2-\beta}d\theta_{12}\sin\,\theta_{12}\int^{\pi-\xi_{2}}_{-(\pi-\xi_{2})}d\phi_{2}+\left.\left.\int^{\pi/2-\beta}_{0}d\theta_{12}\sin\,\theta_{12}\int^{\pi}_{-\pi}d\phi_{2}\right)\right]\right.$
$\displaystyle+\left.\int^{\pi}_{\pi/2}d\beta\,\sin\,\beta\int^{\pi/2}_{\beta-\pi/2}d\theta_{11}\sin\,\theta_{11}\int^{\xi_{1}}_{-\xi_{1}}d\phi_{1}\int^{\pi/2}_{\beta-\pi/2}d\theta_{12}\,\sin\,\theta_{12}\int^{\xi_{2}}_{-\xi_{2}}d\phi_{2}\right\\}$
$\displaystyle\times\Theta(\cos\,\theta_{21})\;\Theta(\cos\,\theta_{22})\cos\,\theta_{12}\left[\cos\,\theta_{21}\Theta(b_{2}-b_{1})\left(\frac{}{}I4(\psi,b_{1})-I4(\psi,b_{2})\right)\right.$
$\displaystyle+\left.\cos\,\theta_{22}\left(\Theta(b_{1}-b_{2})I5(\psi,b_{1})+\Theta(b_{2}-b_{1})I5(\psi,b_{2})\frac{}{}\right)\right].\qquad$
(B6)

In (B6) the functions $I4$ and $I5$ are defined by

$I4(\psi,b)=\csc\,\psi\,\arctan\,(b\csc\,\psi-\cot\,\psi)$

and

$I5(\psi,b)=\log\,(b^{2}-2b\cos\,\psi+1)/2+\cot\,\psi\arctan\,(b\csc\,\psi-\cot\,\psi)-1/2.$

To evaluate the five-fold integration over angles in (B6), ten values of
$\beta$, linearly distributed for $0\leq\beta\leq 0.014$, plus thirty
additional values with logarithmically increasing separations over the range
$0.014<\beta<\pi$, were considered. For each $\beta$ the polar angles
$\theta_{11}$ and $\theta_{12}$ were each allowed to assume thirty distinct
values chosen so as to implement an approximation, via a product of Chebyshev
polynomials over the planar area
$(0<\theta_{11}<\pi/2)\times(0<\theta_{12}<\pi/2)$. The initial double
integral over $\phi_{1}$ and $\phi_{2}$ in (B6) was performed numerically for
each combination of $\theta_{11}$ and $\theta_{12}$ filling a $30\times 30$
matrix for each value of $\beta$ considered. Adding each of these matrices to
its transpose then accounted for the contributions of the trailing terms in
(B5) for which $q_{1}$ and $q_{2}$ were interchanged. With the Chebyshev
coefficients available, the double integrals over $\theta_{11}$ and
$\theta_{12}$ could be evaluated with little additional computational effort.
Finally cubic spline interpolation over forty $\beta$ values facilitated the
final $\beta$ integration giving $0.0414\,R_{3}\,\log\,R_{3}\text{ Ry}$ for
the numeric contribution of the final term of Eq.(28).

## References

  * Foulkes et al. (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73 (2001).
  * Slater (1960) J. C. Slater, _Quantum Theory of Atomic Structure_ , vol. 2 (McGraw Hill, New York, 1960), pages 1-30 and Appendix 22.
  * Porter (2009) J. C. Porter, Phys. Rev. B 80, 205102 (2009).
  * Mayer and Mayer (1940) J. E. Mayer and M. G. Mayer, _Statistical Mechanics of Fluids_ (Wiley and Sons, New York, 1940).
  * Porter (1965a) J. C. Porter, Phys. Rev. 140, A732 (1965a).
  * Oberman et al. (1962) C. Oberman, A. Ron, and J. Dawson, Phys. Fluids 5 (1962).
  * Fetter and Walecka (1971) A. L. Fetter and J. D. Walecka, _Quantum Theory of Many-Particle Systems_ (McGraw-Hill, New York, 1971), page 156.
  * Gell-Mann and Breuckner (1957) M. Gell-Mann and K. A. Breuckner, Phys. Rev. 106 (1957).
  * Carr and Maradudin (1964) J. Carr, W. J. and A. A. Maradudin, Phys. Rev. 133 (1964).
  * Farmelo (2009) G. Farmelo, Physics Today 46 (2009).
  * Korn and Korn (1968) G. A. Korn and T. M. Korn, _Mathematical Handbook for Scientists and Engineers Second Ed._ (McGraw Hill, New York, 1968).
  * Porter (1965b) J. C. Porter, Ph.D. thesis, Purdue University (1965b).
  * (13) Eq.(2) should replace the corresponding expressions in references 3 and 5.
  * N (3) The final term of Eq.(18) was incorrectly omitted from reference 8\.
  * Onsager et al. (1966) L. Onsager, L. Mittag, and M. J. Stephen, Ann Physik 18, 71 (1966).
  * Rajagopal and Kimball (1977) A. K. Rajagopal and J. C. Kimball, Phys. Rev. B 15 (1977).
  * Montgomery and Tidman (1964) D. C. Montgomery and D. A. Tidman, _Plasma Kinetic Theory_ (McGraw Hill, New York, 1964), see Ch.6.