A simple correlation energy functional for the uniform electron gas is derived based on the second-order Moller-Plesset perturbation theory. It can reproduce the known correlation functional in the high-density limit, while in the mid-density range maintaining a good agreement with the near-exact correlation energy of the uniform electron gas to within 2 × 10^{−3} hartree. The correlation energy is a function of a density parameter r_{s} and is of the form $a*ln ( 1 + b r s + b r s 2 ) $. The constants “a” and “b” are derived from the known correlation functional in the high-density limit. Comparisons to the Ceperley-Alder’s near-exact Quantum Monte Carlo results and the Vosko-Wilk-Nusair correlation functional are also reported.

Electron correlation energy functional is an important part of the density functional theory.^{1} In the local-density approximation of a uniform electron gas, general expression of the correlation energy is not known,^{2} except for the high- and low-density limits. There are many parameterization schemes such as the Vosko-Wilk-Nusair^{3} and the Perdew-Wang^{4} forms. These are primarily based on fitting functions to the near-exact Quantum Monte Carlo results of Ceperley and Alder.^{5} However, some of these functional forms such as the Vosko-Wilk-Nusair form are usually analytically complicated and inherently contain numerical irregularities.^{4}

In this report, a simple analytical function describing electron correlation energy within the local-density approximation is derived based on the second-order Moller-Plesset perturbation theory (MP2).^{6–8} The functional is of the form

Not only it approaches the leading term^{9,10}

in the known high-density limit but is also accurate in describing the correlation energy in the mid-density range. Here, $ r s = 4 \pi \rho / 3 \u2212 1 / 3 $ is the Wigner-Seitz density parameter; *ρ* is the electron density. Atomic units are used throughout.

The constants “a” and “b” are

As shown in Fig. 1, when the expression is plotted against the near-exact Quantum Monte Carlo results,^{3,5} it is able to describe the correlation energy in the mid-density range (2 ≤ *r _{s}* ≤ 100) quite accurately, having the root-mean-squared error of only 1.856 × 10

^{−3}hartree.

The VWN functional also agrees with the Monte Carlo data as well. This is to be expected because the VWN parameters have been purposely fitted to these data points, whereas the correlation energy from Eq. (1) is derived independently. Interestingly, the root-mean-squared error of the VWN functional is 1.983 × 10^{−3} hartree, higher than that of Eq. (1). The detailed comparisons are provided in the supplementary material.^{11}

The form of the functional in Eq. (1) can be deduced as follows. From the expression of the MP2 correlation energy for a uniform electron gas^{6}

the integrand can be expanded as a polynomial in *ρ*^{1/3} to obtain the correlation energy at low density. For example,

where *η*_{1}, *η*_{2} depend on the integrating variables *u*, *v*. Then, the integral can be separated into multiple terms, representing each order of *ρ*^{1/3} in the polynomial:

where $ \eta 1 ( I ) , \eta 2 ( I ) $ are now constants.

However, Eq. (6) is not valid at high density because the expression does not approach infinity logarithmically as *ρ* → ∞. In other words, it is inconsistent with Eq. (2). The problem can be remedied by rewriting Eq. (6) inside a logarithmic function. Namely,

In the low-density range *ρ*^{1/3} ≪ 1, Eq. (7) falls back to Eq. (6) because $ln 1 + x \u2245x$ if *x* is appreciably small. Furthermore, in the high-density range *ρ*^{1/3} ≫ 1, the term $ \eta 2 ( I ) a \rho 2 / 3 $ dominates; and $ \epsilon c \u2192aln \eta 2 ( I ) a \rho 2 / 3 =2aln \rho 1 / 3 +Constant$, which is essentially the same expression as Eq. (2). Therefore, Eq. (7) is also in its correct form in the high-density limit.

Note that, the density parameter *r _{s}* is proportional to

*ρ*

^{−1/3}. Therefore, Eq. (7) can be written in terms of

*r*as follows:

_{s} where “b” and “c” are constants. Equation (8) is valid when the two extreme cases are considered. When *r _{s}* ≫ 1, the term $1+ b r s $ dominates, and

*ϵ*is correct up to the leading term $ a b r s $. On the other hand, when

_{c}*r*≪ 1, the term $ c r s 2 $ is “turned on” instead, causing

_{s}*ϵ*to approach the correct high-density original form. However, what happens in the mid-density range requires an assumption about the turning point at which electron gas “chooses” to behave like the low- or the high-density limit.

_{c}In the atomic unit, the Wigner-Seitz density parameter *r _{s}* is defined as $ 4 \pi \rho / 3 \u2212 1 / 3 $. It defines the radius of a sphere which contains exactly one electron. In other words, it tells how far apart the electrons are from their nearest neighbors.

In the SI unit, the parameter *r _{s}* is divided by the Bohr radius

*a*

_{0}, the radius of a hydrogen atom, which can be viewed as the “comfortable” radius of an electron cloud encapsulating a unit of background positive charge. Therefore, the situation

*r*< 1 can be interpreted as the case when the density of electron gas is “squeezed” to be denser than that of a hydrogen atom. This causes the electrons to overlap and the Pauli exclusion principle to be significant, giving rise to the exchange interaction along with other quantum mechanical properties of the system.

_{s}When *r _{s}* ≫ 1, the electrons are far apart. Each electron encapsulates one unit of background positive charge. Following Wigner’s view,

^{12}there is plenty of room for an electron, and there is no need to overlap with its neighbor. In this limit, the system behaves classically, and the Coulomb interaction becomes the dominating term. It is no surprise that the correlation energy (per electron) is proportional to 1/

*r*in the case of the low-density limit.

_{s}Inspecting Eq. (8), one realizes that the parameters “b” and “c” govern the relative strength between the terms $ b r s $ and $ c r s 2 $. For example, if $ b r s $ is the larger one, then the functional *ϵ _{c}* behaves like the low-density limit, and vice versa. At the tipping point between the low- and high-density limits, the two terms must be equal.

If it is assumed that such a tipping point occurs at *r _{s}* = 1, exactly at the Bohr radius, then it follows that

as indicated in Eq. (1). The constants “a” and “b” in Eq. (1) are calculated by requiring that Eq. (1) approaches Eq. (2) as *r _{s}* → 0. Namely, when

*r*is small, the term $ b r s 2 $ becomes the most dominant. In this case,

_{s} yielding $a=\u2212 A 2 ,b=exp \u2212 2 C A $. The values A and C have been analytically calculated^{10} using the perturbation theory. For example, in the paramagnetic case $A= 1 \u2212 ln 2 \pi 2 ,C=\u22120.046\u2009920\u20093$, giving the evaluated “a” and “b” in Eq. (3).

For a spin-polarized uniform electron gas, the correlation energy functional can be written as^{2}

where $\zeta \u2261 \rho \alpha \u2212 \rho \beta \rho , \u03f5 c 0 r s , \u03f5 c 1 r s $ are the spin-polarization parameter, the paramagnetic-, and the ferromagnetic-correlation energy functional, respectively. The function $f \zeta = 1 + \zeta 4 / 3 + 1 \u2212 \zeta 4 / 3 \u2212 2 2 2 1 / 3 \u2212 1 $ is a weighting factor between the two extreme cases: the paramagnetic and the ferromagnetic systems.

$ \u03f5 c 0 r s $ is essentially the same as Eq. (1). For the ferromagnetic $ \u03f5 c 1 r s $, the values A and C from the high-density limit are needed to compute “a” and “b.” Using $A= 1 \u2212 ln 2 2 \pi 2 ,C=\u22120.025\u2009737\u2009 5 10 $ yields

The correlation energy $ \u03f5 c 1 r s $ in the form given by Eq. (1) with the constants “a” and “b” listed above also agrees extremely well with the Monte Carlo results^{3,5} for the ferromagnetic case. The error in this case is well below 0.5 × 10^{−3} hartree as indicated in Fig. 3 of the supplementary material.^{11} Interestingly, this error is also comparable to the much more complicated formula obtained from an interpolation between high- and low-density limits.^{13}

In conclusion, a correlation energy functional is derived based on the second-order Moller-Plesset perturbation theory. The second-order energy correction is first expanded in powers of *ρ*^{1/3} to compute the correlation energy in the low-density range. Without losing significant accuracy, it is argued that when the expansion (up to the second order) is embedded inside a logarithmic function ln(1 + *x*), the analytical form of the high-density limit is recovered. Based on the assumption that the tipping point between the low- and high-density limits occurs when the Wigner–Seitz radius *r _{s}* is exactly one Bohr radius, a simple analytic function $ \u03f5 c r s =aln ( 1 + b r s + b r s 2 ) $ appears naturally. The constants “a” and “b” are then evaluated using the known high-density limit

^{9,10}from the perturbation theory. The correlation energy functional agrees very well with the Monte Carlo results

^{3,5}in the mid- and low-density range (2 ≤

*r*≤ 100), while approaching the known theoretical functional

_{s}*A*ln

*r*+

_{s}*C*in the high-density limit.

Its simple and accurate analytical form is useful as a starting point to further take into account the non-uniform density in molecules or bulk solids within the framework of the generalized gradient approximation (GGA).^{14–16}