Performance of one-body reduced density matrix functionals for the homogeneous electron gas
N.\
Lathiotakis
The subject of this study is the exchange-correlation-energy functional of reduced density matrix functional theory. Approximations of this functional are tested by applying them to the homogeneous electron gas. We find that two approximations recently proposed by Gritsenko, Pernal, and Baerends, J. Chem. Phys., {122}, 204102 (2005), yield considerably better correlation energies and momentum distributions than previously known functionals. %MODIFICATION: %We propose new functionals We introduce modifications to these functionals which, by construction, reproduce the exact correlation energy of the homogeneous electron gas.
jelium.tex
%\documentclass[twocolumn,showpacs,superscriptaddress,prb,aps,floatfix]{revtex4}
%\documentclass[superscriptaddress,prb,preprint,aps,showpacs,floatfix]{revtex4}
\documentclass[preprint,showpacs,superscriptaddress,aps,floatfix]{revtex4}

\usepackage{graphicx}
\usepackage{rotating}
\usepackage{amsmath}

% Useful Macros
\tolerance 1000

% new commands
\newcommand{\br}{{\bf r}}
\newcommand{\bR}{{\bf R}}
\newcommand{\bx}{{\bf x}}
\newcommand{\bk}{{\bf k}}
\newcommand{\brp}{{\bf r'}}
\newcommand{\bRp}{{\bf R'}}

\begin{document}

\title{Performance of one-body reduced density matrix functionals for the homogeneous electron gas}

\author{N.\,N.\,Lathiotakis} 
\affiliation{Institut f{\"u}r Theoretische Physik, Freie Universit{\"a}t Berlin, Arnimallee 14, 
D-14195 Berlin, Germany}
\affiliation{European Theoretical Spectroscopy Facility}

\author{N.\,Helbig}
\affiliation{Institut f{\"u}r Theoretische Physik, Freie Universit{\"a}t Berlin, Arnimallee 14, 
D-14195 Berlin, Germany}
\affiliation{European Theoretical Spectroscopy Facility}
\affiliation{Fritz-Haber-Institut der Max-Planck-Gesellschaft, Faradayweg 4-6, 14195 Berlin,
Germany}
\affiliation{Unit\'e de Physico-Chimie et de Physique des Materi\'eaux, Universit\'e Catholique de Louvain, B-1348 Louvain-la-Neuve, Belgium}

\author{E.\,K.\,U.~Gross}
\affiliation{Institut f{\"u}r Theoretische Physik, Freie Universit{\"a}t Berlin, Arnimallee 14, 
D-14195 Berlin, Germany}
\affiliation{European Theoretical Spectroscopy Facility}

\begin{abstract}
The subject of this study is the exchange-correlation-energy functional of reduced density 
matrix functional theory. Approximations of this functional are tested by applying them to 
the homogeneous electron gas. We find that two approximations recently proposed by
Gritsenko,  Pernal, and Baerends, J. Chem. Phys., {\bf 122}, 204102 (2005),
yield considerably better correlation energies 
and momentum distributions than previously known functionals. 
%MODIFICATION:
%We propose new functionals 
We introduce modifications to these functionals
which, by construction, reproduce the exact correlation energy of the homogeneous electron gas. 
\end{abstract}
\pacs{71.10.-w,71.10.Ca,05.30.Fk,31.25.-v}
\maketitle

\section{Introduction}
Reduced-density-matrix-functional theory (RDMFT) is one possible way
to tackle the problem of electronic correlation. 
It is based on Gilbert's theorem~\cite{gilbert}, which is an extension of the
Hohenberg-Kohn theorem to non-local external potentials. 
Gilbert's theorem guarantees that the ground-state expectation value 
of any observable of a quantum mechanical system is a unique
functional of the ground-state one-body reduced-density matrix (1-RDM).
Thus, the fundamental quantity in RDMFT is the 1-RDM instead of the 
electronic density on which DFT is built upon.  
The properties of the exact 1-RDM have been the subject 
of theoretical studies for a long time~\cite{gilbert,valone,zumbach,mueller}. 
Nevertheless, only relatively recently approximate total-energy functionals of the 1-RDM 
were used in practical 
applications~\cite{bb0,UG,csanyi,csgoe,yasuda,kollmar,kios1,kios2,gritsenko,openshell,pernal,pirin,kollmar2}.
All these approximate functionals give a satisfactory account of electronic correlations in small 
atoms and molecules at the equilibrium distance. The latest generation of 
functionals performs equally well at the molecular dissociation limit. Especially this
latest success makes RDMFT particularly appealing.

A great advantage of RDMFT, 
compared to DFT, is that the exact many-body kinetic energy is easily 
expressed in terms of the 1-RDM. More specifically, the total energy of a
many-electron system in its ground state,
characterized by
an external potential $V(\br)$, can be expressed in terms of the ground-state
1-RDM, $\gamma$ as
\begin{multline}
\label{eq:E_of_g}
E_{\rm tot}\left[  \gamma \right] = E_{\rm ion} +
\int d^3\br \int  d^3\brp\: \delta(\br - \brp ) \left[ -\frac{1}{2} \nabla_{\!\!\br}^2 \right] \gamma(\br, \brp) \\
+\int d^3\br \int  d^3\brp\: \delta(\br - \brp ) V(\br) \gamma(\br, \brp) \\
+\frac{1}{2} \int d^3\br \int  d^3\brp \frac{\gamma(\br,\br)\: \gamma(\brp, \brp)}{|\br-\brp|}
+ E_{\rm xc} \left[  \gamma \right] \,.
\end{multline}
The first term in the above expression is a constant not related to the electronic degrees 
of freedom, for example the ion-ion repulsion energy. The next three terms are,
respectively, the kinetic, the external potential and the electrostatic
energy and they are known, explicit functionals of $\gamma$. 
Finally, the last term is the exchange and correlation (xc) energy 
which contains all the remaining electronic contributions to the total energy. 
The exact form of this functional
is unknown and for practical applications needs to be approximated. 
Most of the approximate functionals for the xc energy that have been introduced 
so far %~\cite{mueller,bb0,UG,csanyi,csgoe,yasuda,kollmar,kios1,kios2,gritsenko,pernal,pirin,kollmar2}
are implicit functionals of $\gamma$. They depend explicitly on
the natural orbitals $\varphi_j$, and the corresponding occupation numbers $n_j$ 
which are defined as the eigenfunctions and eigenvalues of the 1-RDM:
\begin{equation}
\int d^3\br'\: \gamma(\br, \brp) \: \varphi_j (\brp) = n_j\: \varphi_j (\br)\,.
\end{equation}

Viewed as a functional of arbitrary $\gamma$ in an appropriately defined domain, 
the functional given in  Eq.~(\ref{eq:E_of_g}) has a minimum value at the ground-state $\gamma$.
This appropriate domain is defined through subsidiary conditions for $\gamma$ known
as $N$-representability conditions. It was shown by Coleman~\cite{coleman} that there
 are two such conditions for $\gamma$, and they concern the occupation numbers
\begin{equation}
\label{eq:conditions}
\sum_{j=1}^\infty n_j = N\,, \quad 0 \leq n_j \leq 1\,,
\end{equation}
where $N$ is the total number of electrons. These conditions guarantee that $\gamma$ corresponds 
to either a pure many-body state or an ensemble of pure states. The first 
condition can be enforced in the minimization with respect to $\gamma$ through
the Lagrange-multipliers method. In that way, the quantity to be minimized becomes
\begin{equation}
F = E_{\rm tot} - \mu \left(\sum_{j=1}^\infty n_j - N\right)\,,
\label{eq:F}
\end{equation}
where $\mu$ is the corresponding Lagrange multiplier. $\mu$ was shown to be equal to the
chemical potential, i.e. the derivative of the total energy with respect to the total number
of electrons~\cite{our_gap}. Interestingly, this allows one to exploit the discontinuity of
$\mu$ as a function of the particle number for the calculation of the fundamental 
gap of materials and is a motivation for the development of RDMFT schemes for periodic systems in
order to address questions like the semiconductor and insulator gaps. 

The second of the $N$-representability conditions~(\ref{eq:conditions}) has a dramatic consequence: 
it allows for border minima in the occupation number optimization. In other words, it
allows for the the existence of a subset of the optimal occupation numbers 
which are exactly equal to either one or zero and do not satisfy the condition 
$\partial F/ \partial n_j = 0$.  
We refer to the corresponding states as pinned states.
It is rather unlikely for the exact theory to produce pinned states for most 
systems of interest. For an occupation number being exactly equal to one, the corresponding
natural orbital would be present in all determinants of the full CI expansion
with non-zero coefficients. A situation like that has not been found for small atoms
and molecules or for the HEG, where the exact solution can be calculated.
While this is true for the exact xc functional, 
all the approximate functionals which are considered 
in this work, yield pinned states with $n_j=1$ for all systems we applied them to, except
for the two-electron systems. These pinned states are core states and,
in the exact theory, they correspond to occupation numbers which are marginally smaller 
than one. Hence, as far as the optimal $\gamma$ is concerned, the approximate result, $n_j=1$, 
is perfectly satisfactory.
The important implication, however, is that $\partial F /  \partial n_j \neq 0$ for the pinned states,
and consequently $\delta F / \delta \gamma(\br, \brp)  \neq 0$ at the optimal $\gamma$.

A number of approximate functionals for $E_{\rm xc}$, including those of interest in the present 
work, can be cast into the form
\begin{widetext}
\begin{equation}
\label{eq:Exc}
E_{\rm xc} \left[  \gamma \right] = E_{\rm xc} \left[ \{ n_j \}, \{ \varphi_j \} \right]=
-\frac{1}{2} \sum_{j,l=1}^{\infty}  \int d^3\br \int  d^3\brp \: f(n_j,n_l) \:
\frac{\varphi_j^*(\br)\: \varphi_l^*(\brp)\: \varphi_l(\br)\: \varphi_j(\brp)}{|\br-\brp|} \,,
\end{equation}
\end{widetext}
i.e. they have the form of the usual Hartree-Fock exchange modified by the function
$f(n_j,n_l)$ of the occupation numbers. The first such approximation was introduced
by M\"uller~\cite{mueller} and corresponds to the function $f(n_j,n_l) = \sqrt {n_j n_l}$.
M\"uller considered a more general exponent for the occupation number product in the 
exchange-like term and found an optimal exponent of 1/2.  By modelling 
the exchange and correlation hole, Buijse and Baerends~\cite{bb0} arrived 
at the same functional. Goedecker and Umrigar~\cite{UG} (GU) considered a 
modification by explicitly removing the self-interaction (SI) terms. They also
presented~\cite{UG} a direct minimization with respect to the natural orbitals and the
occupation numbers and found correlation energies for small atomic 
systems which are in very good agreement with the exact results.  Later however, it
was realized that the GU
functional fails to reproduce the correct dissociation limit for
small molecules~\cite{staroverov,herbert}. On the other hand, the M\"uller 
functional yields the correct dissociation limit but, in all cases, overestimates
substantially the correlation energy~\cite{staroverov,herbert}. 

In the last decade, several other functionals of the 1-RDM have been 
introduced~\cite{csgoe,yasuda,kollmar,kios1,kios2,gritsenko,openshell,pernal,pirin,kollmar2} and 
applied to atomic and molecular systems. 
Recently, Gritsenko et al.~\cite{gritsenko} proposed 
improved 1-RDM functionals based on a hierarchy of repulsive corrections 
to the M\"uller functional. In that way, they attempted to correct the 
overcorrelation of this functional. The functionals corresponding to these hierarchical corrections, 
are called BBC1, BBC2, and BBC3. For all these functionals, it is essential to
divide the natural orbitals
into strongly and weakly occupied ones. This distinction appears naturally
for finite systems since usually a subset of the orbitals corresponds to occupation numbers close to one, and the rest to occupation numbers close to zero. For the BBC1 and BBC2  functionals, the function $f(n_j,n_l)$ is:
\begin{widetext}
\begin{equation}
\mbox{BBC1:\ \ \ \  }f(n_j,n_l) = \left\{  \begin{array}{rl}
                                 -\sqrt{n_j\: n_l}\,, & j\neq l, \mbox{and\ } j,l\mbox{\ weakly occupied,} \\
                                 \sqrt{n_j\: n_l}\,, & \mbox{\ otherwise,} 
                             \end{array}
\right.
\end{equation}
\begin{equation}
\mbox{BBC2:\ \ \ \  }f(n_j,n_l) = \left\{  \begin{array}{rl}
                                 -\sqrt{n_j\: n_l}\,, & j\neq l, \mbox{and\ } j,l\mbox{\ weakly occupied,} \\
                                 n_j\: n_l\,, & j\neq l, \mbox{and\ } j,l\mbox{\ strongly occupied,} \\
                                 \sqrt{n_j\: n_l}\,, & \mbox{\ otherwise.} 
                             \end{array}
\right.
\end{equation}
\end{widetext}
Finally, in the BBC3 functional the anti-bonding orbital is treated as strongly occupied
orbital. Additionally, the self-interaction terms are removed 
as in the GU functional, except for the bonding and anti-bonding orbitals. 
Gritsenko et al.~\cite{gritsenko} applied the BBC functionals to 
diatomic molecules and showed that they give an accurate description of
these molecules at both the equilibrium distance and the dissociation limit.

There is a strong motivation for the extension and application of 1-RDM functionals
to solid-state systems. This motivation stems from the success of these functionals 
in the description of electron correlation for finite systems, as well as the difficulties
of DFT methods in describing certain materials and properties such as the band gap of
semiconductors and insulators~\cite{ortiz,filippi,terakura,sharma} or the band
width of the conduction band in Na~\cite{sodium}.

A very important prototype system, which serves as a benchmark
for the performance of approximate 1-RDM functionals when applied to
periodic systems, is the homogeneous electron gas (HEG). Furthermore, this 
system can serve as a laboratory for the development of approximate functionals, 
in a fashion similar to DFT. As a consequence of translational invariance, the natural orbitals
can be chosen as plane waves and the 
search for the ground state 1-RDM is restricted to the optimization of the momentum distribution 
$n(\bk)$, which is the 
occupation number that corresponds to the plane-wave natural orbital with wavevector $\bk$. 
An important point to note is that the self-interaction terms for the plane-wave 
natural orbitals vanish. Consequently, the GU and the M\"uller functionals are identical.
Finally, approximations that involve a special treatment of single orbitals have zero effect 
in the continuous wavevector case. Thus, the BBC3 functional is identical to BBC2.

As a consequence of the rotational invariance of the HEG the occupation numbers have the 
property $n(\bk)=n(k)$, i.e. all the plane-wave natural orbitals corresponding to the same 
absolute value $k$ are degenerate (with respect to the occupation number). This allows one
to perform unitary transformations among the degenerate plane waves leading, e.g., to 
angular-momentum eigenfunctions
\begin{equation}
\label{eq:sphwaves}
\varphi_{klm} (\br) = j_l(k r)\: Y_{lm}(\Omega) 
\end{equation}
for the natural orbitals, where $j_l(k r)$ are spherical Bessel functions and $Y_{lm}$
spherical harmonics.
 Since these functions are localized in real space it is conceivable to include 
self-interaction corrections in terms of the natural orbitals (\ref{eq:sphwaves}) for the HEG. 
To our knowledge, this possibility has not been explored so far and is also beyond the aim of the
present work.

Choosing the natural orbitals as plane waves, the 
1-RDM of the HEG can be written as
\begin{equation}
\gamma(\br, \brp) = \frac{2}{V}\sum_{\bk} n(k)\: e^{i\bk(\br-\brp)}\,, 
\end{equation}
where $V$ is the volume of the system ($V\rightarrow \infty$).
Substituting this expression in Eqs.~(\ref{eq:E_of_g}) and~(\ref{eq:Exc}) we obtain for the total energy
\begin{widetext}
\begin{equation}
E_{\rm tot} = 2 \sum_{\bk_1} \frac{\bk_1^2}{2} \: n({k_1}) - \frac{1}{V}
\sum_{\bk_1, \bk_2} f\left(n({k_1}\right), n({k_2}))\: 
\frac{4\pi}{\left| \bk_1 - \bk_2 \right|^2 +\alpha^2} \,,
\label{eq:etotheg}
\end{equation}
where $\bk_1$, $\bk_2$ are wavevector indices, and $\alpha$ is a small quantity (usually 
$\sim 10^{-8}$) included 
for numerical stability.  As in 
Hartree-Fock, the external potential energy and the electronic Coulomb repulsion energy,
i.e. the third and fourth terms in Eq.~(\ref{eq:E_of_g})
cancel exactly with the first term i.e. the ion-ion interaction. 
The quantity $F$ of Eq.~(\ref{eq:F}) per particle then becomes
\begin{multline}
\label{eq:energy}
\frac{F}{N} = \frac{3}{2k_{\rm F}^3} \int_0^\infty dk_1\: 
k_1^2\,(k_1^2 -2\mu)\: n(k_1) \\ + 
\frac{3}{4\pi k_{\rm F}^3} \int_0^\infty dk_1 \int_0^\infty dk_2 \: k_1\: k_2\:
\log \left[ \frac{\left( k_1 -k_2  \right)^2+\alpha^2}
{\left( k_1 + k_2\right)^2+\alpha^2} \right]  f(n(k_1), n(k_2)) + \mu \,,
\end{multline}
\end{widetext}
where $k_{\rm F}=(9\pi /4)^{1/3}r_s^{-1}$ is the 
Fermi-wavevector of the non-interacting HEG and $r_s$ is the radius (in atomic units)
of the sphere with volume equal to the volume per electron.

Cioslowski and Pernal~\cite{ciospernal} applied the M\"uller functional
to the HEG and calculated analytically the resulting momentum distribution 
\begin{equation}
\label{eq:kios_per}
n(k) = 512 \pi \rho \; (1+4k^2 )^{-4}\,,
\end{equation}
where  $\rho $ is the electron-density-per-spin,
$\rho= 3 (8\pi r_s^{3})^{-1}$. The corresponding total energy 
per particle is independent of the density and equal to $-1/8$ Hartree. 
It is obvious that the solution of the Eq.~(\ref{eq:kios_per}) is consistent 
with the second $N$-representability constraint of Eq.~(\ref{eq:conditions})
only for $\rho\leq (512\pi)^{-1}$, i.e. $r_s \geq 5.77$. 
In other words, the M\"uller functional gives a solution without pinned states only
for $r_s \geq 5.77$.  In addition, Cioslowski and Pernal demonstrated
that the Oxford-Lieb~\cite{oxlieb} bound is violated for $\rho \geq 1.65\times10^{-3}$,
i.e. $r_s \leq 4.167$. The solution with pinned states for $r_s < 5.77$, was 
calculated by Cs\'anyi and Arias~\cite{csanyi}. More specifically, for
$r_s < 5.77$, one gets an optimal momentum distribution $n(k)$ with 
$n(k)=1$ for $k$ below a certain value $k_p$ and fractional $n(k) <1$ for $k>k_p$. 
This behavior is in complete analogy to the case of finite systems for the M\"uller
functional. Unfortunately, it is in conflict with the fact that the 
exact momentum distribution~\cite{exactmom2,gori} is a monotonically decreasing function of $k$ and
is strictly smaller than 1, i.e. there are no pinned states. Additionally,
the exact momentum distribution is concave for $k<k_{\rm F}$, it shows a
discontinuity at $k_{\rm F}$, and for $k>k_{\rm F}$ it goes to zero asymptotically. The size of the 
discontinuity is decreasing with $r_s$.

In addition to the M\"uller functional, Cs\'anyi and Arias~\cite{csanyi} considered a similar 
functional derived from a tensor product expansion of the two-body density matrix, which they
called Corrected Hartree-Fock (CHF). Unfortunately,
CHF gives zero correlation for the HEG in the high-density limit ($r_s\rightarrow 0$), coinciding
with Hartree-Fock. In the opposite limit, it strongly overcorrelates 
giving the same results as the M\"uller functional. In the intermediate region,
including the metallic densities,
the result for the correlation energy is close to the exact but its dependence on $r_s$ is
monotonically decreasing instead of increasing. In an attempt to improve over the M\"uller functional
and CHF, Cs\'anyi, Goedecker and Arias considered an improved tensor product expansion of
the two-particle density matrix~\cite{csgoe}. The resulting functional, which is called CGA,
performs very well in the high-density regime and significantly better than the previous two
functionals in the region of metallic densities. The deviation from the exact correlation 
energy increases
with $r_s$ and at higher densities CGA coincides with CHF and the M\"uller functional.

%As already mentioned, BBC3 coincides with BBC2 for the HEG when one assumes that
%the natural orbitals are plane waves.
%For the application of BBC1,2 one has to appropriately define 
%the strongly and weakly occupied orbitals, since the corrections
%are based on that distinction.
%We choose a critical wavevector $k_c$ which is equal to the Fermi wavevector,
%$k_{\rm F}$ of the non-interacting HEG, 
%below which we consider the states to be 
%strongly occupied. 
In the present work, we apply the BBC1 and BBC2 functionals
of Gritsenko et al.~\cite{gritsenko} to the HEG and compare with previous
functionals as far as the resulting correlation energies are concerned.
We also investigate other features like the resemblance of the resulting momentum 
distribution to the exact and the state-pinning.  In order to apply the BBC1,2 
functionals to the case of the HEG through Eq.~(\ref{eq:energy}), we need to 
distinguish between strongly and weakly occupied 
orbitals. For finite systems,  Gritsenko et al.~\cite{gritsenko} chose the first $N/2$ 
natural orbitals to be strongly occupied. In complete analogy, we can use a critical 
wavevector $k_c=k_{\rm F}$ below which all states are assumed to be strongly occupied 
while above they are weakly occupied.

An additional goal of the present work is to demonstrate that the HEG can be used to develop 
functionals suitable for metallic systems.  The idea is to modify the BBC1 functional 
in such a way that it yields the exact correlation energy for the HEG. This 
is achieved in two different ways: (i) For each given density we choose $k_c$ such that 
BBC1 reproduces the exact correlation energy of the HEG at that density. We call this functional 
$k_c$-functional. (ii) We introduce a function $s(r_s)$ multiplying the xc terms of Eq.~(\ref{eq:Exc})
for two weakly occupied orbitals, keeping $k_c=k_{\rm F}$. In this way, we replace the sign 
change of the BBC1 functional with the parameter $s$. Accordingly, we call this functional $s$-functional. 

In the following section we present details of the numerical implementation as well as the results 
of applying the BBCs and the 
%MODIFICATION
%two new functionals 
$k_c$- and $s$-functional
%END MODIFICATION
to the HEG.  

%We also 
%are presented in Section~\ref{sec:results}. We close this paper with  concluding remarks in 
%Section~\ref{sec:conclusion}.

%\section{RDMFT for the homogeneous electron gas\label{sec:theory}}

%Most of the functionals of the 1-RDM that have been introduced so far 
%can be written in the general form 
%\begin{multline}
%E_{\rm tot} = 
%\sum_i n_i \: \phi_i^*(\br) \left[  -\frac{\nabla^2}{2} + V_{\rm ext}(\br)  \right]
%\phi_i(\br)  \\
%+ \sum_{ij} n_i n_j \int\int 
%\frac{\phi_i^*(\br) \: \phi_j^*(\brp) \; \phi_i(\br) \: \phi_j(\brp)}{|\br -\brp|} d^3\br d^3\brp  \\
%- \sum_{ij} f(n_i, n_j) \int\int 
%\frac{\phi_i^*(\br) \: \phi_j^*(\brp) \; \phi_j(\br) \: \phi_i(\brp)}{|\br -\brp|} d^3\br d^3\brp
%\,,
%\label{eq:mueller}
%\end{multline}
%where $V_{\rm ext}(\br)$ is an external local potential like the 
%attractive ion potential. 
%The last term is exchange-like and the functional of Eq.~(\ref{eq:mueller})
%reduces to  Hartree-Fock for   $f(n_i, n_j)=n_i n_j$.
%For the M\"uller~\cite{mueller,bb0} and Goedecker-Umrigar~\cite{UG}  functionals 
%$f(n_i, n_j)=\sqrt{n_i n_j}$, but for the 
%Goedecker-Umrigar functional the self-interaction terms are excluded from
%the direct Coulomb and the exchange-like terms.
%\end{widetext}
%The minimization procedure involves the functional derivative 
%\begin{widetext}
%\begin{multline}
%\frac{\delta (F/N)}{\delta n(k)}
%= \frac{3}{2k_{\rm F}^3}\: k^2 \,(k^2 -2\mu) \\
%+\frac{3}{4\pi k_{\rm F}^3} \int_0^\infty dk^\prime \: k\: k^\prime \:
%\log \left[ \frac{\left( k -k^\prime  \right)^2+\alpha^2}
%{\left( k +k^\prime  \right)^2+\alpha^2} \right]
%2\, f^\prime(n(k), n(k^\prime)) 
%\,,
%\end{multline}
%%\end{widetext}
%where $f^\prime(n(k), n(k^\prime))$ is the derivative of $f$ with respect 
%to its first argument.
%For example, for the M\"uller functional we have 
%\begin{equation}
%f^\prime(n(\epsilon), n(\epsilon^\prime)) = \frac{1}{2} \sqrt{ \frac{n(\epsilon^\prime)}{n(\epsilon)}}\,,
%\end{equation}
%For example, for the  BBC1 functional
%\begin{equation}
% f^\prime(n(k), n(k^\prime)) =  \left\{
% \begin{array}{rl}-\frac{1}{2}\sqrt{ \frac{n(k)}{ n(k^\prime)}}\,, & 
%                                 \mbox{$\bk$, $\bk^\prime$ weakly occupied} \\
%                   \frac{1}{2}\sqrt{ \frac{n(k)}{ n(k^\prime)}}\,, & \mbox{otherwise.} \\
% \end{array} \right.
%\end{equation}
%and for the BBC2
%\begin{equation}
% f^\prime(n(\epsilon), n(\epsilon^\prime)) =  \left\{
% \begin{array}{rl}-\frac{1}{2}\sqrt{\frac{n(\epsilon)}{ n(\epsilon^\prime)}}\,, & 
%                                  \mbox{$\epsilon$, $\epsilon^\prime$ weakly occupied} \\
%                   n(\epsilon^\prime ) \,, & 
%                                  \mbox{$\epsilon$, $\epsilon^\prime$ strongly occupied} \\
%    \frac{1}{2}\sqrt{\frac{n(\epsilon)}{ n(\epsilon^\prime)}}\,, & \mbox{otherwise} \\
% \end{array} \right.
%\end{equation}


\section{Numerical Implementation, results\label{sec:results}}

The minimization of the energy expression~(\ref{eq:energy}) with respect to
$n(k)$ is performed using the steepest-descent method. We choose to work in 
energy-space instead
of $k$-space, i.e. we perform the variable substitution $\epsilon=k^2/2$ and solve numerically 
the minimization problem for  $n(\epsilon)$. Working in energy-space 
rather than $k$-space improves the 
stability of the numerical treatment. The energy $\epsilon$ is discretized using a
double-logarithmic mesh centered at $\epsilon_{\rm F} = k_{\rm F}^2/2$, where
the occupation varies the most. 
The upper limit of integration is chosen such that  
the momentum-distribution function has dropped to values smaller than $10^{-6}$. 
The double integration with respect to the energy is carried out using 
an adaptive grid technique capable of treating integrable singularities
like the logarithmic singularity of the present problem. 
The values of the energy-distribution function $n(\epsilon)$ in between the mesh points,
necessary for the adaptive grid method, are obtained from an interpolation scheme. 
The $N$-representability constraint $0 \leq n(\epsilon) \leq 1$ is implemented through 
the substitution $n(\epsilon) = \sin^2 [ \pi \: \theta(\epsilon) /2 ]$ and
variation with respect to $\theta(\epsilon)$. Extra care is required to avoid $n(\epsilon)$ being
falsely pinned to 0 or 1. Indeed, if for a particular 
$\epsilon$, $n(\epsilon)$ gets very close to 1 or 0 during the variation it would stay pinned at that
point. The variation with
respect to the Lagrange-multiplier $\mu$ is implemented as an external iterative procedure, 
thus achieving convergence for each value of $\mu$. 
The correct value of $\mu$ is selected by requiring the 
momentum-distribution function to integrate to the correct number of electrons for a given 
value of $r_s$. Finally, we found that a reasonable value for the parameter $\alpha$ 
in Eqs.~(\ref{eq:etotheg}) and (\ref{eq:energy}) is $10^{-8}$, which avoids both numerical 
problems as well as the dependence of the results on $\alpha$.

\subsection{Application of the BBC functionals to the HEG}

\begin{figure}[t]
\includegraphics[width=0.9\columnwidth,clip]{Jelium_energy_vs_rs.eps}
\caption{\label{fig:correl} (Color online) The correlation energy of the HEG as a function of $r_s$
calculated with the BBC1 and BBC2 functionals compared to various other calculations.
The Monte Carlo result represents the Perdew-Wang fit~\cite{PW} of the DMC 
data of Ortiz and Ballone~\cite{exactmom2}. 
The dotted line corresponds to the numerical results by Cs\'anyi and Arias~\cite{csanyi} 
employing the M\"uller functional,
for $r_s < 5.77$. Its continuation, the dash-dotted line, for $r_s > 5.77$, stands for the
analytical results of Cioslowski and Pernal~\cite{ciospernal} employing the same functional.
The results for the CHF functional~\cite{csanyi} as well as the CGA~\cite{csgoe} are
also shown.}
\end{figure}

In Fig.~\ref{fig:correl},  we show the correlation energy of the HEG as a function of $r_s$. The
correlation energy calculated with the BBC1 and BBC2 functionals is significantly closer to
the exact than any other functional over the whole range of $r_s$.  Both functionals
also seem to reproduce the correct asymptotic limit of zero correlation for
the dilute HEG. For small densities up to metallic densities, the
BBC functionals under-correlate, i.e. the absolute value of the correlation energy is too small. 
In the dense limit, they over-correlate and the crossover is
at around $r_s$=0.5 and 0.3 for the BBC1 and BBC2, respectively. Unfortunately, in the area of low metallic
densities both functionals yield correlation energies which deviate from the 
exact values by 50\%. 
In absolute numbers, the error of BBC1 and BBC2 is of the same order as the RPA result~\cite{PW}. 
Nevertheless, in the range $0.1<r_s<1$, the BBC functionals perform remarkably well. 
Compared to all previous 1-RDM functionals, BBC1 and BBC2 offer a much better account of the 
electron correlation for the HEG. Although less accurate than the CGA in the high density 
region, they perform better for metallic densities and they reproduce the limit of zero 
correlation at the dilute HEG limit where the M\"uller functional, CHF and CGA 
fail. 

A feature of the exact momentum distribution, namely the discontinuity at the 
Fermi wavevector $k_{\rm F}$,  is reproduced by the BBC functionals. 
The discontinuity is more pronounced for the BBC1 functional, as can be seen in 
Fig.~\ref{fig:momdistr},
where we plotted the momentum distribution of the HEG with $r_s=1$ and $r_s=5$  
using the M\"uller functional as well as BBC1 and BBC2. Contrary to BBC1 and BBC2 
the M\"uller functional does not yield a discontinuity. To our knowledge,
there is no report of any other 1-RDM functional reproducing this feature of the HEG.
To extract the size of the discontinuity quantitatively, 
we used two energy mesh points very close to $\epsilon_{\rm F}$ 
(at a distance of $\pm 10^{-8}\epsilon_{\rm F}$).
In Fig.~\ref{fig:zeta1}, we plot the size of the 
discontinuity $\Delta n$ as a function of $r_s$. As we see, it increases monotonically 
with $r_s$ and has the tendency to saturate for large $r_s$ for both BBC1 and BBC2. 
For BBC2 the discontinuity is substantially smaller than for BBC1 and it goes to zero 
at $0.6 < r_s < 0.7$. 
Both the size and the dependence on $r_s$ are in complete disagreement
with the exact theory, where $\Delta n$ is substantially bigger and decreases with $r_s$
as one can see from the two fits to the DMC data~\cite{exactmom2,gori}.

\begin{figure*}
\begin{tabular}{cc}
\includegraphics[angle=0,width=0.45\textwidth,clip]{n_rs_1.0_BB0.eps}& \ \ \ \ \
\includegraphics[angle=0,width=0.45\textwidth,clip]{n_rs_5.0_BB0.eps}\\
\includegraphics[angle=0,width=0.45\textwidth,clip]{n_rs_1.0_BB1.eps}& \ \ \ \ \ 
\includegraphics[angle=0,width=0.45\textwidth,clip]{n_rs_5.0_BB1.eps}\\
\includegraphics[angle=0,width=0.45\textwidth,clip]{n_rs_1.0_BB2.eps}&  \ \ \ \ \ 
\includegraphics[angle=0,width=0.45\textwidth,clip]{n_rs_5.0_BB2.eps}\\
\end{tabular}
\caption{\label{fig:momdistr} (Color online) The momentum distribution $n(k)$ of the HEG for $r_s$=1 (left)
and $r_s=5$ (right) 
calculated with the M\"uller functional, BBC1, and BBC2. 
BBC1 and BBC2 show a discontinuity of the momentum distribution at $k_{\rm F}$. For
comparison we include the fit to the DMC data
of Ortiz-Ballone~\cite{exactmom2}.  The derivative 
$\delta (F/N) / \delta n(k) $ is also plotted. 
Note that for small $k$ the derivative is not zero and  $n(k)$ is pinned at one.
}
\end{figure*}

\begin{figure}
\includegraphics[angle=0,width=0.9\columnwidth,clip]{deltan_rs.eps}
\caption{\label{fig:zeta1} (Color online) The discontinuity $\Delta n$ of the momentum distribution at $k_{\rm F}$ for BBC1 and
BBC2 as a function of $r_s$, compared to the fits of Ortiz and Ballone (OB)~\cite{exactmom2}, and
Gori-Giorgi and Ziesche (GZ)~\cite{gori} to the DMC data.}
\end{figure}

%One of the features of the M\"uller functional is that the momentum distribution 
%is equal to 1 for small wavevectors $k$. For these wavevectors, $n(k)$ is determined
%by the $N$-representability constraint and therefore the solution is not variational
%($\delta(F/N)/\delta n(k) \neq 0$). We found that this feature is preserved in the
%BBC1 and BBC2 corrections. This can be seen in Fig.~\ref{fig:momdistr} where we also show
%the derivative $\delta(F/N)/\delta n(k)$ as a function of $k$.
%For small $k$ the derivative is negative indicating that, without the 
%$N$-representability constraint, $n(k)$ would increase beyond 1. Then, after a critical
%wavevector $k_p$ the derivative is zero and the momentum distribution becomes smaller than one.

We now turn to the question of state pinning. As we see in Fig.~\ref{fig:momdistr}, state
pinning is a common feature of all the functionals we employed. 
To verify that the states are truly pinned, we plot the functional derivative 
$\delta F /\delta n(k)$ which is non-zero for pinned states. As one expects,
the number of pinned states decreases with increasing $r_s$.
We define a wavevector $k_p$ below 
which the corresponding states are pinned, i.e. $n(k)=1$ for $k < k_p$.
In Fig.~\ref{fig:ep_rs}, we plot $k_p$
as a function of $r_s$ for the M\"uller and the two BBC functionals.
For all three functionals, $k_p$ decreases monotonically with $r_s$. 
For the M\"uller functional our numerical calculation confirms the analytic 
result~\cite{ciospernal} that above a critical value of $r_s=5.77$ there are
no pinned states, and therefore $k_p$ goes to zero at this value. 
This can already be seen in Fig.~\ref{fig:momdistr}, where, for the M\"uller functional 
at $r_s=5$, the derivative is very close to zero even for small wavevectors.
Interestingly, for BBC1 and BBC2 we found no such critical value up to
$r_s=20$. Indeed, for both BBC1
and BBC2 the decrease of $k_p$ is much smaller than for the M\"uller functional,
$k_p$ being almost constant for BBC2.

As we have seen, the performance of the BBC1,2 functionals for the HEG is 
improved significantly compared to previous functionals. This is especially 
remarkable, given that they were originally constructed to describe the 
dissociation of small molecules. 

\begin{figure}
\includegraphics[angle=0,width=0.9\columnwidth,clip]{ep_rs.eps}
\caption{\label{fig:ep_rs} (Color online) The wavevector $k_p$, where the optimal
momentum distribution becomes fractional, as a function of $r_s$.
 }
\end{figure}

\subsection{Improved functionals for the HEG}

%MODIFICATION:
%In this section we are going beyond the BBC functionals. We 
%make use of the application of the BBCs to the HEG in order to develop functionals
%that reproduce the correlation energy of the HEG exactly.
%The new functionals are based on the BBC1 form where an additional function of
%$r_s$ is introduced such that the exact correlation energy of the HEG is
%reproduced for each $r_s$. 
In this section, we attempt to improve over the BBC functionals
for the HEG. We present two functionals which are simple modifications of the 
BBC1 functional. In both cases, the BBC1 functional is generalized 
by introducing an additional function of $r_s$
such that the correlation energy of the HEG is reproduced exactly for each $r_s$.
%END MODIFICATION
As exact results we regard the Perdew-Wang fit~\cite{PW} of the 
correlation energy obtained from diffusion Monte Carlo (DMC) calculations by 
Ceperley and Alder~\cite{exactmom1} and
Ortiz and Ballone~\cite{exactmom2}. The two Monte Carlo calculations yield almost identical 
correlation energies.

For the first functional, we adjust the critical wavevector $k_c$, which 
is used to distinguish between strongly and weakly occupies states, instead
of using $k_c=k_{\rm F}$, as in BBC1. We call this functional the $k_c$-functional.
The corresponding function $f$ in Eq.~(\ref{eq:Exc}) then reads
\begin{equation}
 f(n({k_1}), n({k_2})) =  \left\{
 \begin{array}{rl}-\sqrt{n({k_1})\: n({k_2})}\,, & \quad \mbox{$k_1,k_2>k_c(r_s)$,} \\
                   \sqrt{n({k_1})\: n({k_2})}\,, & \quad \mbox{otherwise.} \\
 \end{array} \right.
\label{eq:ourfunct_0}
\end{equation}
We perform  the fitting of $k_c$ 
over the range of metallic densities, $0.5 \leq r_s \leq 5$. 
The results are compiled in Table~\ref{tab:ec}. For $r_s = 0.5$, $k_c\approx k_{\rm F}$, since this
point is almost exactly reproduced by BBC1 (see Fig.~\ref{fig:correl}). For $r_s > 0.5$, $k_c$ is a monotonically increasing
function of $r_s$.  Fitting $k_c$ has a strong impact on the momentum distribution which is 
displayed in Fig.~\ref{fig:ec_rs5}, for $r_s=1$ and $r_s=5$. 
It is not surprising that the discontinuity is 
displaced from $k_{\rm F}$ to $k_c$. Additionally, its size is reduced significantly compared
to BBC1. Both the displacement of the discontinuity and the decrease in the
step size are in disagreement with the exact result.

\begin{table}[t]
\begin{tabular}{c||c|c|c|c|c|c}
\hline\hline
$r_s$ & 0.5 & 1.0 & 2.0 & 3.0 & 4.0 & 5.0 \\ \hline
$k_c$/$k_{\rm F}$ & 0.994 & 1.032 & 1.085 & 1.122 & 1.155 & 1.172 \\
\hline\hline
\end{tabular}
\caption{\label{tab:ec} The the critical wavevector $k_c$ for some metallic 
densities. }
\end{table}

An alternative idea is to keep $k_c=k_{\rm F}$ fixed and consider a fitting parameter 
$s$ multiplying the exchange-like terms when both states, $\bk_1$ and $\bk_2$, are weakly occupied, i.e.
\begin{equation}
 f(n({k_1}), n({k_2})) =  \left\{
 \begin{array}{rl}-s(r_s)\,\sqrt{n({k_1})\: n({k_2})}\,, & \quad \mbox{$k_1,k_2>k_{\rm F}$,} \\
                   \sqrt{n({k_1})\: n({k_2})}\,, & \quad \mbox{otherwise.} \\
 \end{array} \right.
\label{eq:ourfunct}
\end{equation}
In this way, the parameter $s$ which is a kind of strength of the 
xc terms becomes a function of $r_s$, as indicated 
in Eq.~(\ref{eq:ourfunct}). We call this functional the $s$-functional.
 The values of $s$ for the fitting to the  Ortiz and 
Ballone DMC results are included in Table~\ref{table:s_functional}.
As we see, $s$ varies between $4$ and $-0.26$ 
for the range of densities, $0.1 \leq r_s \leq 10$, we considered.
In Fig.~\ref{fig:s_rs}, we show the resulting, monotonically decreasing
function $s(r_s)$.

\begin{figure}
\begin{tabular}{cc}
\includegraphics[angle=0,width=0.45\columnwidth,clip]{n_rs_1.0_BB1_ec.eps}  \ \ \ & \ \ \ 
\includegraphics[angle=0,width=0.45\columnwidth,clip]{n_rs_5.0_BB1_ec.eps} \\
\end{tabular}
\caption{\label{fig:ec_rs5} (Color online) The momentum distribution for the $k_c$-functional 
compared with the BBC1 and the fit to the DMC data of Ortiz-Ballone~\cite{exactmom2}, 
for  $r_s=1$ and $r_s=5$. For the $k_c$-functional the 
discontinuity is moved from $k_{\rm F}$ to $k_c$, i.e to $1.032k_{\rm F}$ and $1.172k_{\rm F}$
respectively.}
\end{figure}

\begin{table}
\setlength{\tabcolsep}{0.6truecm}
\begin{tabular}{cc}
\hline\hline
\begin{tabular}{cc}
$r_s$ &  s \\ \hline
0.1 & 4.913 \\
0.2 & 2.751 \\
0.3 & 1.867 \\
0.4 & 1.390 \\
0.5 & 1.087 \\
0.6 & 0.877 \\ 
0.7 & 0.727 \\
0.8 & 0.602 \\ 
\end{tabular} & 
\begin{tabular}{cc}
$r_s$ & s \\ \hline
1.0 & 0.435 \\
1.5 & 0.190 \\
2.0 & 0.059 \\
3.0 & -0.074 \\
4.0 & -0.146 \\
5.0 & -0.189 \\
7.0 & -0.234 \\
10.0 & -0.263 \\ 
\end{tabular} \\ \hline\hline
\end{tabular}
\caption{\label{table:s_functional}The fitted values of $s$ for various values of $r_s$ for the 
$s$-functional. $s$ was fitted to reproduce the correlation energies of the DMC calculation 
of Ortiz and Ballone~\cite{exactmom2}.}
\end{table}

\begin{figure}
\includegraphics[angle=0,width=0.7\columnwidth,clip]{s_rs.eps}
\caption{\label{fig:s_rs} (Color online) The dependence of the fitting parameter 
$s$ on $r_s$ resulting
from fits to the Perdew-Wang parameterization~\cite{PW} of the correlation energy from
two different sets of DMC results:
the Ceperley and Alder~\cite{exactmom1} and the Ortiz and Ballone~\cite{exactmom2}.} 
\end{figure}

From Fig.~\ref{fig:n_k_fitts} one can see that one of the advantages of fitting $s$ instead of $k_c$
is that the discontinuity of the momentum distribution remains fixed at $k_{\rm F}$. 
In addition, its size is almost constant ($\approx 0.2$) as a function of $r_s$.
As we see in Fig.~\ref{fig:zeta1}, the exact discontinuity is significantly higher and it is a decreasing
function of $r_s$. Therefore, concerning the size of the discontinuity, the $s$-functional 
does not improve over the BBC1, which,
in the dilute gas limit, is close to the exact result.
However, the increasing behavior of BBC1  is improved by the $s$-functional which yields 
discontinuities that are almost constant as a function of $r_s$.

In conclusion, the $s$-functional results in momentum distributions that resemble 
the exact ones over the whole range of $r_s$ more closely than any of the other functionals
considered here.

\begin{figure}
\begin{tabular}{cc}
\includegraphics[angle=0,width=0.45\columnwidth,clip]{n_rs_1.0_BB1_s.eps} \ \ \ & \ \ \ 
\includegraphics[angle=0,width=0.45\columnwidth,clip]{n_rs_5.0_BB1_s.eps} \\
\end{tabular}
\caption{\label{fig:n_k_fitts} (Color online) The momentum distribution given by the $s$-functional,
compared with the BBC1 and the fit to the DMC data of Ortiz-Ballone~\cite{exactmom2},
for $r_s=1$ ($s=0.435$), and $r_s=5$ ($s=-0.189$).  The discontinuity remains at 
$k_{\rm F}$ and is approximately equal to 0.2.}
\end{figure}

%It is possible to extend the functional we just 
%introduced for the HEG so it is applicable 
%to a larger variety of system by circumventing 
%One way is the following:
%Lets define the quantity $\alpha$ as the ratio $\alpha=E_c/E_{\rm kin}^{\rm (HF)}$, where 
%$E_c$ is the correlation energy and $E_{\rm kin}^{\rm (HF)}$ is the kinetic Hartree Fock energy.
%Obviously, $\alpha$ is also a function of $r_s$. We will demonstrate that a monotonic 
%dependence $s=s(\alpha)$ can be established. The gain is that instead of 
%$r_s$,  $alpha$ is well defined
%quantity for all kinds of systems. The function $s=s(\alpha)$ can be exploited in 
%a self-consistent
%manner: a starting value for  $\alpha$ is chosen and in every iteration it 
%is recalculated 
%by the resulting correlation energy. We determined the function $s=s(\alpha)$ using 
%both the CA and OB Monte Carlo results. Our numerical results are in 
%Table~\ref{table:s_a} and are also shown in Fig.~\ref{fig:s_a}. This functional was tested on the 
%HEG and for the whole ranges of densities considered ($0.1<rs<10$). The implicit function 
%$s^{(n+1)}=f(s^{(n)})$ where by $f$ we mean the procedure of determining the next value of 
%$s$ in the self consistent scheme has always a fixed attractive point in the whole range 
%of densities considered. The functional we just described is exact for the HEG and is applicable 
%to any kind of system in the same fashion that LDA functionals originate from the 
%correlation energy of the HEG. 
%
%\begin{table}
%\begin{tabular}{cc}
%\hline
%\multicolumn{2}{c}{Ceperley-Alder} \\
%\hline\hline
%$\alpha$ & $s$ \\
%-1.0939794$\times 10^{-3}$  & 4.80     \\
%-3.66072678$\times 10^{-3}$ &  2.6366     \\
%-7.331623097$\times 10^{-3}$ &  1.7585     \\
%-1.73353975$\times 10^{-2}$ &  0.981     \\
%-3.02585$\times 10^{-2}$    &  0.6253    \\
%-5.40964154$\times 10^{-2}$ &  0.34203     \\
%-0.1620329d0   & -0.0064     \\
%-0.3008926129 & -0.1176     \\
%-0.638405507 & -0.202     \\
%-1.680826122 & -0.24     \\
%\hline
%\end{tabular}
%
%\begin{tabular}{cc}
%\hline
%\multicolumn{2}{c}{Ortiz Ballone} \\
%\hline\hline
%$\alpha$ & $s$ \\
%-1.086$\times 10^{-3}$ & 4.913     \\
%-3.616$\times 10^{-3}$ & 2.751    \\
%%-7.215$\times 10^{-3}$ & 1.867    \\
%-1.170$\times 10^{-2}$ & 1.390    \\
%-1.696$\times 10^{-2}$ & 1.087    \\
%-2.291$\times 10^{-2}$ & 0.877    \\
%-2.949$\times 10^{-2}$ & 0.727   \\
%-3.664$\times 10^{-2}$ & 0.602    \\
%-5.252$\times 10^{-2}$ & 0.435   \\
%-0.100 & 0.190    \\
%-0.157 & 0.059    \\
%-0.293 & -0.074    \\
%-0.453 & -0.146    \\
%-0.632 & -0.189    \\
%-1.036 & -0.234    \\
%-1.729 & -0.263    \\
%\hline
%\end{tabular}
%\caption{\label{table:s_a}  Tabulated values of the dependence of the fitting strength $s$ 
%as a function of a, i.e. the ratio 
%of the correlation energy over the Hartree-Fock kinetic energy, for two sets of correlation energies:
%the Ceperley and Alder and Ortiz and Ballone Dynamic Monte-Carlo results.}
%\end{table}
%
%\begin{figure}
%\includegraphics[angle=0,width=0.9\columnwidth,clip]{s_a.eps}
%\caption{\label{fig:s_a} The dependence of the fitting strength $s$ as a function of a, i.e. the ratio 
%of the correlation energy over the Hartree-Fock kinetic energy, for two sets of correlation energies:
%the Ceperley and Alder and Ortiz and Ballone Dynamic Monte-Carlo results.}
%\end{figure}

\section{conclusion\label{sec:conclusion}}
We have applied a variety of 1-RDM functionals to the  
HEG. We show that the BBC functionals~\cite{gritsenko} yield a significant 
improvement over previous functionals as far as the correlation energy is concerned. In
addition, they yield a discontinuity of the momentum distribution at
the Fermi wavevector in resemblance of the exact HEG theory. However,
the size and the dependence on the density of this discontinuity are not 
in agreement with the quantum Monte Carlo results. 

%MODIFICATION
By introducing an appropriately fitted function of $r_s$ in the BBC1 functional, we
demonstrate that the exact correlation energy of the HEG can be reproduced with a
smooth and monotonic fitting function.
%END MODIFICATION
For this function, we either use the critical wavevector $k_c(r_s)$ which
distinguishes between the strongly and weakly occupied states, or 
a strength $s(r_s)$ multiplying the exchange-like terms for two weakly occupied states. 
Both of these procedures were applied to the BBC1 functional. The two functionals,
yielding by construction the exact correlation energy of the HEG, are assessed by the quality
of the resulting momentum distributions. We show that choosing the second procedure, i.e. the 
$s$-functional, is superior to fitting $k_c$. The discontinuity $\Delta n$ of the momentum distribution 
resulting from the $s$-functional is nearly constant as a function of $r_s$ and hence
represents a significant improvement over BBC1 and BBC2.
%MODIFICATION
However, the momentum distribution obtained by the $s$-functional still deviates significantly
from the exact one. To remedy this, more complicated strategies have to be considered, possibly 
with the introduction of more fitting parameters.
%END MODIFICATION

Our functional of choice, being derived from the HEG, is expected to yield 
good results for metallic systems. The application to finite as well as 
non-metallic periodic systems is not straightforward because $r_s$ is not well-defined
in these cases. Hence, the necessity arises to map $r_s$ onto other quantities characterizing 
these systems, or to involve a LDA-type prescription relating $s(r_s)$ to the local density.

\begin{acknowledgements}
This work was supported in part by the  Deutsche Forschungsgemeischaft within the program SPP 1145,
by the EXCITING Research and Training Network, and by the 
EU's Sixth Framework Program through the Nanoquanta Network of Excellence (NMP4-CT-2004-500198).
\end{acknowledgements}

\begin{thebibliography}{10}
%\addcontentsline{toc}{chapter}{Bibliography}
\bibitem{gilbert}
T.\,L. Gilbert, Phys. Rev. B{\bf 12}, 2111 (1975).

\bibitem{valone} 
S.\,M. Valone, J. Chem. Phys. {\bf 73}, 1344 (1980); {\bf 73}, 4653 (1980).

\bibitem{zumbach}
G. Zumbach and K. Maschke, J. Chem. Phys. {\bf 82}, 5604 (1985).

\bibitem{mueller}
A.\,M.\,K. M\"uller, Phys. Lett. A {\bf 105}, 446 (1984). 

\bibitem{bb0}
M.\,A. Buijse, PhD Thesis, Vrije Universiteit Amsterdam (1991);
M.\,A. Buijse, E.\,J. Baerends, Mol. Phys. {\bf 100}, 401 (2002).

\bibitem{UG}
S. Goedecker, C.\,J. Umrigar, Phys. Rev.~Lett.~{\bf 81}, 866 (1998).

\bibitem{csanyi}
G. Cs\'anyi, T.\,A. Arias, Phys. Rev. B{\bf 61} (2000) 7348.

\bibitem{csgoe} 
G. Cs\'anyi, S. Goedecker,  T.\,A. Arias, Phys. Rev. A{\bf 65} 032510 (2002).

\bibitem{yasuda}
K. Yasuda, Phys. Rev. Lett. {\bf 88}, 053001 (2002).

\bibitem{kollmar} 
C. Kollmar, J. Chem. Phys. {\bf 121}, 11581 (2004).

\bibitem{kios1} 
J. Cioslowski, K. Pernal, M. Buchowiecki, J. Chem. Phys. {\bf 119}, 6443 (2003).

\bibitem{kios2} 
K. Pernal, J. Cioslowski, J. Chem. Phys. {\bf 120}, 5987 (2004).

\bibitem{gritsenko}
O. Gritsenko, K. Pernal, E.\,J.\, Baerends, J. Chem. Phys., {\bf 122}, 204102 (2005).

\bibitem{openshell}  
N.\,N. Lathiotakis, N. Helbig, and E.\,K.\,U. Gross, Phys. Rev. A {\bf 72}, 030501 (2005).

\bibitem{pernal} J. Cioslowski and K. Pernal, Phys. Rev. B {\bf 71}, 113103 (2005).

\bibitem{pirin} P. Leiva, and M. Piris, J. Chem. Phys., {\bf 123}, 214102 (2005).

\bibitem{kollmar2} C. Kollmar, J. Chem. Phys., {\bf 125}, 084108 (2006).

\bibitem{coleman} A. Coleman, Rev. Mod. Phys. {\bf 35},  668  (1963).

\bibitem{our_gap} 
N. Helbig, N.\,N. Lathiotakis, M. Albrecht, E.\,K.\,U. Gross, 
%MIDIFICATION
%(to be published), cond-mat/0504436.
Europhys. Lett. {\bf 77}, 67003 (2007).

\bibitem{staroverov} 
V.\,N. Staroverov and G.\,E. Scuseria, J. Chem. Phys. {\bf 117}, 2489 (2002).

\bibitem{herbert} 
J.\,M. Herbert, J.\,E. Harriman, Chem. Phys. Lett. {\bf 382}, 142 (2003).

\bibitem{ortiz}
G. Ortiz, Phys. Rev. B {\bf 45}, 11328 (1992); 

\bibitem{filippi}
C. Filippi, D.\,J. Singh, C.\,J. Umrigar,
Phys. Rev. B {\bf 50} 14947 (1994); 

\bibitem{terakura}
K. Terakura, T. Oguchi, A.\,R. Williams, J. K\"ubler,
Phys. Rev. B {\bf 30} 4734 (1984); 

\bibitem{sharma}
Sharma S., Dewhurst J. K., and Ambrosch-Draxl C., Phys. 
Rev. Lett. {\bf 95}, 136402 (2005).

\bibitem{sodium}
H. Yasuhara, S. Yoshinaga, and M. Higuchi, Phys. Rev. Lett. {\bf 83}, 3250 (1999).

\bibitem{ciospernal} 
J. Cioslowski and K. Pernal, J. Chem. Phys. {\bf 111}, 3396 (1999).

\bibitem{oxlieb} 
E.\,H. Lieb and S. Oxford, Int. J. Quantum Chem. {\bf 19}, 427 (1981).

%\bibitem{exactmom1} 
%B. Farid, V. Heine, G. E. Engel, and I. J. Robertson, Phys. Rev. B{ \bf 48}, 11~602 (1993).

\bibitem{exactmom2} 
G. Ortiz, P. Ballone, Phys. Rev. B{\bf 50}, 1391 (1994);
{\bf 56}, 9970 (1997).

\bibitem{gori}
P. Gori-Giorgi and P. Ziesche, Phys. Rev. B{\bf 66}, 235116 (2002).

\bibitem{PW} J.\,P. Perdew, Y. Wang, Phys. Rev. B {\bf 45}, 13244 (1992).

\bibitem{exactmom1} 
D.\,M. Ceperley and B.\,J. Alder, Phys. Rev.~Lett.~{\bf 45}, 566
(1980).

\end{thebibliography}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

%\bibitem{ciosyas} J. Cioslowski, K. Pernal, J. Chem. Phys. {\bf 117}, 67 (2002).
%
%\bibitem{extens} J. Cioslowski, K. Pernal, J. Chem. Phys. {\bf 120}, 10364 (2004).
%
%\bibitem{svane} 
%A. Svane, Phys. Rev B{\bf 53}, 4275 (1996).
%
%\bibitem{goda}
%S. Goedecker, C.\,J. Umrigar, Phys. Rev. A{\bf 55}, 1765 (1997).
%
%\bibitem{cohen}
%A.\,J.\, Cohen, E.\,J. Baerends, Chem. Phys. Lett. {\bf 364}, 409 (2002).
%
%
\end{document}

cond-mat0605531.extracted_bib
\begin{thebibliography}{10}
%\addcontentsline{toc}{chapter}{Bibliography}
\bibitem{gilbert}
T.\,L. Gilbert, Phys. Rev. B{\bf 12}, 2111 (1975).

\bibitem{valone} 
S.\,M. Valone, J. Chem. Phys. {\bf 73}, 1344 (1980); {\bf 73}, 4653 (1980).

\bibitem{zumbach}
G. Zumbach and K. Maschke, J. Chem. Phys. {\bf 82}, 5604 (1985).

\bibitem{mueller}
A.\,M.\,K. M\"uller, Phys. Lett. A {\bf 105}, 446 (1984). 

\bibitem{bb0}
M.\,A. Buijse, PhD Thesis, Vrije Universiteit Amsterdam (1991);
M.\,A. Buijse, E.\,J. Baerends, Mol. Phys. {\bf 100}, 401 (2002).

\bibitem{UG}
S. Goedecker, C.\,J. Umrigar, Phys. Rev.~Lett.~{\bf 81}, 866 (1998).

\bibitem{csanyi}
G. Cs\'anyi, T.\,A. Arias, Phys. Rev. B{\bf 61} (2000) 7348.

\bibitem{csgoe} 
G. Cs\'anyi, S. Goedecker,  T.\,A. Arias, Phys. Rev. A{\bf 65} 032510 (2002).

\bibitem{yasuda}
K. Yasuda, Phys. Rev. Lett. {\bf 88}, 053001 (2002).

\bibitem{kollmar} 
C. Kollmar, J. Chem. Phys. {\bf 121}, 11581 (2004).

\bibitem{kios1} 
J. Cioslowski, K. Pernal, M. Buchowiecki, J. Chem. Phys. {\bf 119}, 6443 (2003).

\bibitem{kios2} 
K. Pernal, J. Cioslowski, J. Chem. Phys. {\bf 120}, 5987 (2004).

\bibitem{gritsenko}
O. Gritsenko, K. Pernal, E.\,J.\, Baerends, J. Chem. Phys., {\bf 122}, 204102 (2005).

\bibitem{openshell}  
N.\,N. Lathiotakis, N. Helbig, and E.\,K.\,U. Gross, Phys. Rev. A {\bf 72}, 030501 (2005).

\bibitem{pernal} J. Cioslowski and K. Pernal, Phys. Rev. B {\bf 71}, 113103 (2005).

\bibitem{pirin} P. Leiva, and M. Piris, J. Chem. Phys., {\bf 123}, 214102 (2005).

\bibitem{kollmar2} C. Kollmar, J. Chem. Phys., {\bf 125}, 084108 (2006).

\bibitem{coleman} A. Coleman, Rev. Mod. Phys. {\bf 35},  668  (1963).

\bibitem{our_gap} 
N. Helbig, N.\,N. Lathiotakis, M. Albrecht, E.\,K.\,U. Gross, 
%MIDIFICATION
%(to be published), cond-mat/0504436.
Europhys. Lett. {\bf 77}, 67003 (2007).

\bibitem{staroverov} 
V.\,N. Staroverov and G.\,E. Scuseria, J. Chem. Phys. {\bf 117}, 2489 (2002).

\bibitem{herbert} 
J.\,M. Herbert, J.\,E. Harriman, Chem. Phys. Lett. {\bf 382}, 142 (2003).

\bibitem{ortiz}
G. Ortiz, Phys. Rev. B {\bf 45}, 11328 (1992); 

\bibitem{filippi}
C. Filippi, D.\,J. Singh, C.\,J. Umrigar,
Phys. Rev. B {\bf 50} 14947 (1994); 

\bibitem{terakura}
K. Terakura, T. Oguchi, A.\,R. Williams, J. K\"ubler,
Phys. Rev. B {\bf 30} 4734 (1984); 

\bibitem{sharma}
Sharma S., Dewhurst J. K., and Ambrosch-Draxl C., Phys. 
Rev. Lett. {\bf 95}, 136402 (2005).

\bibitem{sodium}
H. Yasuhara, S. Yoshinaga, and M. Higuchi, Phys. Rev. Lett. {\bf 83}, 3250 (1999).

\bibitem{ciospernal} 
J. Cioslowski and K. Pernal, J. Chem. Phys. {\bf 111}, 3396 (1999).

\bibitem{oxlieb} 
E.\,H. Lieb and S. Oxford, Int. J. Quantum Chem. {\bf 19}, 427 (1981).

%\bibitem{exactmom1} 
%B. Farid, V. Heine, G. E. Engel, and I. J. Robertson, Phys. Rev. B{ \bf 48}, 11~602 (1993).

\bibitem{exactmom2} 
G. Ortiz, P. Ballone, Phys. Rev. B{\bf 50}, 1391 (1994);
{\bf 56}, 9970 (1997).

\bibitem{gori}
P. Gori-Giorgi and P. Ziesche, Phys. Rev. B{\bf 66}, 235116 (2002).

\bibitem{PW} J.\,P. Perdew, Y. Wang, Phys. Rev. B {\bf 45}, 13244 (1992).

\bibitem{exactmom1} 
D.\,M. Ceperley and B.\,J. Alder, Phys. Rev.~Lett.~{\bf 45}, 566
(1980).

\end{thebibliography}