In a recent paper, Lucco Castello *et al.* (arXiv:2107.03537) provided an accurate parameterization of classical one-component plasma bridge functions that was embedded in a novel dielectric scheme for strongly coupled electron liquids. Here, this approach is rigorously formulated, its set of equations is formally derived, and its numerical algorithm is scrutinized. A systematic comparison with available and new path integral Monte Carlo simulations reveals a rather unprecedented agreement especially in terms of the interaction energy and the long wavelength limit of the static local field correction.

## I. INTRODUCTION

The quantum one-component plasma or jellium (also referred to as uniform electron gas or homogeneous electron gas) is universally acclaimed as one of the most important idealized systems in condensed matter physics,^{1,2} statistical mechanics,^{3–6} and quantum chemistry.^{7,8} The jellium can serve as a model for many metals under the simplifying assumption that the charge density of the positive ionic background is uniformly smeared out.^{2–6} Hence, initial investigations focused on its ground state at metallic densities.^{9} The ground-state studies have led to remarkable insights such as the BCS superconductivity theory,^{10} Landau’s Fermi-liquid theory,^{11} the Bohm–Pines quasi-particle picture of collective excitations,^{12} Wigner’s crystallization paradigm,^{13} and even the Kohn–Sham formulation of density functional theory.^{14}

Increasing interest in warm dense matter,^{15,16} an exotic state of high-temperature highly compressed matter that is encountered in dense astrophysical objects (giant planet interiors, brown or white dwarfs, and neutron star crusts) and ultra-fast laser heating or shock compression of metals,^{17,18} has provided the impetus for the worldwide intense research activity targeted at the high-density finite temperature jellium.^{19} Only recently, an accurate description of the uniform electron gas has been achieved in this regime on the basis of a combination of novel path integral Monte Carlo (PIMC) methods.^{16,20–23}

Owing to the fact that low-density finite temperature inhomogeneous electron systems are inaccessible even to the contemporary state-of-the-art experiments, much less attention has been paid to the strongly coupled uniform electron liquid. In fact, there have been rather few computational and theoretical studies, despite the numerous intriguing physical phenomena that have been speculated to manifest themselves at strong coupling. These include the possibility of a charge-density wave instability,^{24–26} the emergence of a spin-density wave,^{26–28} the possibility of a continuous paramagnetic to ferromagnetic transition,^{29–31} and the emergence of an incipient excitonic mode.^{32,33} It should be pointed out that this constitutes a particularly challenging regime for theoretical approaches due to the complex interplay between quantum effects (exchange degeneracy and diffraction), thermal excitations, and strong Coulomb correlations.

The thermodynamic properties and static structure of strongly coupled electron liquids were recently probed by extensive PIMC simulations at finite temperatures.^{34} A dielectric formalism scheme that handles quantum mechanical effects at the random phase approximation level, assumes a frequency independent local field correction, and treats strong Coulomb correlations within the classical hypernetted-chain (HNC) approximation^{35} was revealed to provide consistently accurate predictions for the interaction energies, although there was ample room for improvement regarding structural properties.^{34} From the integral equation theory of classical liquids,^{36–38} it is known that as the coupling increases, the bridge diagrams, which are neglected in the HNC approximation, have an increasing importance for two-particle correlations. Hence, the HNC-based scheme is expected to gradually become inaccurate as crystallization is approached.

In a recent letter,^{39} we extracted the bridge functions of the classical one-component plasma at multiple thermodynamic state points, spanning the whole dense liquid region, from specially designed molecular dynamics simulations. With this input, we constructed a very accurate closed-form bridge function parameterization that covers the entire non-trivial range. This analytic description was incorporated into a novel dielectric scheme that naturally extends the HNC dielectric scheme by including the exact Coulomb bridge function. A comparison with available and new PIMC simulations of the strongly coupled electron liquid revealed that the novel scheme leads to significant improvements over the HNC-based scheme. It should be noted that this marked the first incorporation of bridge functions into a dielectric scheme, although bridge functions had been earlier embedded in variants of the so-called classical mapping method.^{40–42}

In the present work, the integral equation theory based dielectric scheme is formally introduced and its emerging self-consistent set of equations is derived. The basic assumptions and main drawbacks of this dielectric scheme are discussed together with possible further refinements. Mathematical techniques are presented that decrease the computational cost and facilitate convergence. Extensive PIMC simulations are analyzed that extend our current picture of the finite temperature strongly coupled electron liquid. A systematic comparison is carried out with other dielectric schemes and “exact” PIMC simulations in terms of interaction energies, static structure factors, static local field corrections, and static density responses.

## II. THEORETICAL

### A. The paramagnetic electron liquid

The uniform electron fluid (UEF) is a homogeneous quantum mechanical system that consists of electrons that are immersed in a rigid charge neutralizing background.^{43} In other words, the UEF constitutes the quantum mechanical analog of the classical one-component plasma (OCP). Since the electronic interactions involve not only Coulomb interactions but also exchange effects, the distinction between spin-up and spin-down electrons makes the UEF essentially a two-component system. As a consequence, the thermodynamic state points of the UEF are specified by three dimensionless parameters:^{16,19,43} the *quantum coupling or Brueckner parameter* *r*_{s} = *d*/*a*_{B}, where *d* is the Wigner–Seitz radius *d* = (4*πn*/3)^{−1/3} and *a*_{B} = *ℏ*^{2}/(*m*_{e}*e*^{2}) is the first Bohr radius; the *degeneracy parameter* *θ* = *k*_{B}*T*/*E*_{F}, where $EF=[(6\pi 2n\u2191)2/3/2](\u210f2/me)$ is the Fermi energy with respect to the Fermi wave-vector of spin-up electrons $kF\u2191=(6\pi 2n\u2191)1/3$; and the *spin polarization parameter* *ξ* = (*n*^{↑} − *n*^{↓})/*n* that is constrained within 0 ≤ *ξ* ≤ 1 under the spin-up choice *n*^{↑} ≥ *n*^{↓}. In the above description, *ℏ* is the reduced Planck constant, *k*_{B} is the Boltzmann constant, *e* is the electron charge, *m*_{e} is the electron mass, *n* = *n*^{↑} + *n*^{↓} is the electron density, and *T* is the temperature. In contrast to the UEF, the thermodynamic state points of the classical OCP are specified by one dimensionless parameter: the *classical coupling parameter* Γ = *e*^{2}/(*dk*_{B}*T*)^{19,43–45} for which Γ = 2*λ*^{2}*r*_{s}/*θ* with $\lambda 3=(kFd)\u22123=4/(9\pi )$.

In what follows, we shall restrict ourselves to the spin-unpolarized (paramagnetic) case, *n*^{↑} = *n*^{↓} or *ξ* = 0. The ideal (Lindhard) density response of free electrons, i.e., in the absence of Coulomb interactions, has the form^{19,43,46}

with *ϵ*(** q**) =

*ℏ*

^{2}

*q*

^{2}/(2

*m*

_{e}) being the electron kinetic energy,

*η*→ 0

^{+}owing to the adiabatic switching of the perturbation at

*t*→ −

*∞*that ensures causality, and

*f*

_{0}(

**) being the Fermi–Dirac distribution function that is given by**

*q*with $\mu \u0304=\beta \mu $ being the reduced chemical potential (*μ* is the chemical potential and 1/*β* = *k*_{B}*T*) that is determined by the normalization condition $\u222b[d3q/(2\pi )3]f0q=n/2$.

Finally, we shall mainly focus on the high-degeneracy moderate-density range that lies beyond the warm dense matter regime^{15,16,19} but simultaneously prior to the Wigner crystallization.^{13,47} Roughly demarcating the upper warm dense matter boundary with *r*_{s} ∼ 10,^{16,19} the above implies that we are interested in *r*_{s} ∼ 10–100 and *θ* ∼ 1. In this range, the UEF behaves as an electron liquid whose properties are determined by the interplay between strong Coulomb correlations and quantum effects.

### B. The dielectric formalism

The dielectric formalism is based on the following: (i) the quantum fluctuation dissipation theorem that connects the dynamic structure factor (DSF) *S*(** k**,

*ω*) with the imaginary part of the exact density–density response function

*χ*(

**,**

*k**ω*); (ii) the general relation that expresses

*χ*(

**,**

*k**ω*) in terms of the ideal (Lindhard) density response

*χ*

_{0}(

**,**

*k**ω*), which introduces the unknown dynamic local field correction (LFC)

*G*(

**,**

*k**ω*); and (iii) a functional relation connecting

*G*(

**,**

*k**ω*) with the static structure factor (SSF)

*S*(

**) that is obtained by perturbative quantum/classical BBGKY approaches or non-perturbative integral equation approaches.**

*k*^{19,43}The first two building blocks are exact and common to all the dielectric schemes, whereas the third building block is approximate and differs.

The quantum fluctuation dissipation theorem has the general form^{46,48}

Integration over the frequency domain together with the zero frequency moment rule *S*(** k**) =

*∫S*(

**,**

*k**ω*)

*dω*leads to $S(k)=\u2212[\u210f/(2\pi n)]\u222bcoth(\beta \u210f\omega )/2I{\chi (k,\omega )}d\omega $.

^{19,43}The extension of the

*χ*(

**,**

*k**ω*) domain with the aid of analytic continuation leads to the complex valued $\chi \u0303(k,z)$ and allows for integral computation with contour integration techniques.

^{49}Owing to the infinitely many poles of the integrand, the integral is converted to an infinite sum, i.e.,

where *ω*_{l} = 2*πl*/(*βℏ*) are the so-called bosonic Matsubara frequencies.^{19,49}

In addition, within the linear response theory, the polarization potential approach reveals that *χ*(** k**,

*ω*) can always be expressed in terms of

*χ*

_{0}(

**,**

*k**ω*) and the unknown LFC as

^{50–52}

with *U*(** k**) = 4

*πe*

^{2}/

*k*

^{2}being the regularized Fourier transform of the Coulomb pair interaction energy.

Furthermore, all dielectric formalism schemes approximate the LFC as a SSF functional, leading to closed self-consistent approaches.^{19,43} Most rigorous approaches treat quantum effects on the random phase approximation (RPA) level and correlation effects in a classical manner, such as the STLS scheme,^{53,54} VS scheme,^{55,56} CA scheme,^{57,58} and HNC scheme.^{34,35} These lead to a frequency independent LFC and closures of the form *G*(** k**) ≡

*F*[

*S*(

**k**)]. The same applies for the recently developed semi-empirical ESA scheme that utilizes established asymptotic limits and incorporates exact MC simulation results for the warm dense matter regime in order to construct an accurate frequency independent LFC.

^{59,60}The qSTLS scheme constitutes a notable exception since it captures beyond-RPA quantum effects by truncating the quantum BBGKY hierarchy within the Wigner representation at its first member with the introduction of the standard STLS closure condition.

^{61,62}This yields a dynamic LFC and a

*G*(

**,**

*k**ω*) ≡

*F*[

*S*(

**k**),

*ω*] closure.

Combining the above expressions, regardless of the approximation scheme, the dielectric formalism ultimately leads to a functional equation of the form

which can be solved in an iterative manner.

### C. The integral equation theory of liquids

For one-component pair-interacting classical systems, the integral equation theory (IET) of liquids consists of two formally exact equations, the Ornstein–Zernike (OZ) integral equation and the non-linear closure equation,^{36–38}

where *g*(*r*) is the radial distribution function or pair correlation function, *h*(*r*) = *g*(*r*) − 1 is the total correlation function, *c*(*r*) is the direct correlation function, and *b*(*r*) is the bridge function. It is noted that capital notations *H*(** k**),

*C*(

**),**

*k**B*(

**) are reserved for the respective Fourier transforms. Fourier transformed radial distribution functions do not emerge in what follows and, thus, should not be confused with the frequency independent LFC**

*k**G*(

**).**

*k*The bridge function is formally defined by a virial-type expansion whose unknown coefficients are given by multi-dimensional integrals with the integrands being products of Mayer functions *f*(*r*) = exp[ −*βu*(*r*)] − 1. The operation of topological reduction allows us to recast the f-bond expansion as an h-bond expansion,^{36,37} a convenient resummation for long range interactions. However, slow convergence implies that rigorous computation of the exact *b*[*h*] functional is impossible. Hence, most theoretical attempts have focused on the formulation of phenomenological *b*[*h*] closures, typically of the form *b*(*h* − *c*).^{63} The hypernetted-chain (HNC) approximation, which assumes that *b*[*h*] ≡ 0, serves as a prominent example. For the classical OCP, the HNC yields near-exact structural and thermodynamic results for Γ ≲ 1,^{64} but its accuracy strongly degrades as crystallization is approached.^{65–67}

A computationally costly but far more accurate alternative concerns the indirect bridge function extraction from computer simulations.^{68–70} Extraction efforts for the classical OCP date back to the 1980s,^{71–73} but their deficiencies have been documented.^{39,74,75} A very recent study managed to obtain very accurate OCP bridge functions for 17 state points Γ = 10–170 and construct an accurate analytic representation valid along the short range and intermediate range.^{39,75} In reduced units *x* = *r*/*d*, this parameterization reads as

where *b*_{S}(*x*, Γ) is the short range monotonic bridge function, *b*_{I}(*x*, Γ) is the intermediate range oscillatory decaying bridge function, *f*(*x*) is a sigmoid switching function, $si(\Gamma )=\u2211j=03sij\Gamma (ln\Gamma )j,li(\Gamma )=\u2211j=04lij\Gamma 1/6(ln\Gamma )j$ are monotonic functions of Γ, and the $sij,lij$ coefficients are tabulated.^{39,75} The closed system of Eqs. (6)–(8) leads to a near-exact description of strong correlations in the classical OCP.

### D. The IET dielectric scheme

The derivation of the functional closure of the IET-based scheme is based on the technique originally developed by Tanaka for the HNC-based scheme.^{35} The difference is that the IET non-linear closure equation for an arbitrary known *b*[*h*] ≡ *b*(*r*) function is used instead of the HNC non-linear closure for *b*(*r*) ≡ 0. In this section, the basic assumptions and mathematical steps will be outlined in sufficient detail for this article to remain self-contained.

The starting point is the classical (*θ* → *∞*) fluctuation dissipation theorem, which reads as^{37,48}

The integration over the frequency domain, together with the zero frequency moment rule *S*(** k**) =

*∫S*(

**,**

*k**ω*)

*dω*and the $\chi (k,0)=\pi \u22121\u222b[I{\chi (k,\omega )}/\omega ]d\omega $ Kramers–Kronig causality relation, leads to

^{76,77}

Substituting for the general form of the density–density response function *χ*(** k**,

*ω*), approximating the LFC with a frequency independent value

*G*(

**,**

*k**ω*) ≡

*G*(

**), employing the Maxwellian result for the static limit of the ideal classical density–density response**

*k**χ*

_{0}(

**, 0) = −**

*k**nβ*, using the connecting relation

*S*(

**) = 1 +**

*k**nH*(

**), and invoking the Fourier transformed OZ equation**

*k**H*(

**) =**

*k**C*(

**)/[1 −**

*k**nC*(

**)], one obtains**

*k*^{77–79}

Having established the above linear *G* ≡ *F*[*C*] functional, the objective is to exploit the two exact IET equations to derive an integral *G* ≡ *F*[*S*] functional. Application of the gradient operator to the non-linear closure equation and the OZ equation [see Eqs. (6) and (7)] yields

Equating the right-hand sides of the above equations, two of the convolution terms cancel out each other, resulting in

whose Fourier transform (after using the multiplication, convolution, and differentiation properties) reads as

We proceed with substituting the linear *G* ≡ *F*[*C*] functional relation on both sides, operating with (** k**·) on both sides, and solving for

*G*(

**), which leads to**

*k*Using again the connecting relation *S*(** k**) = 1 +

*nH*(

**) in order to dispose of the Fourier transformed total correlation function and substituting for**

*k**U*(

**) = 4**

*k**πe*

^{2}/

*k*

^{2}ultimately result in the sought-for IET functional

As expected, the IET functional collapses to the HNC functional for *b*(*r*) ≡ 0, which explicitly reads as^{34,35}

It is rather straightforward to express the IET functional in terms of the HNC functional,

Considering the short range of the OCP bridge function and the long range of Coulomb interactions, this relation suggests that the IET scheme results in small corrections to the LFC of the HNC scheme at the long and intermediate wavelength regions.

## III. NUMERICAL

### A. On the treatment of the IET closure functional

Considering the isotropy of homogeneous electron liquids, the IET closure functional *G*_{IET}(** k**) is formally expressed as a triple integral of the type

*∫d*

^{3}

*q*(

**·**

*k***)**

*q**I*(|

**−**

*k***|)**

*q**J*(|

**|); see Eq. (9). Such a triple integral type is also encountered in the course of the Wertheim–Thiele derivation of the exact solution of the Percus–Yevick approximation for hard spheres,**

*q*^{80,81}the Laplace transform based derivation of the exact solution of the soft mean spherical approximation for classical plasmas,

^{82,83}and the mathematical treatment of the asymptotic memory kernel in mode coupling theories of classical supercooled liquids.

^{84–86}This triple integral can be directly converted to a double integral with the introduction of azimuthally expanded two-center bipolar coordinates, which is equivalent to the sequential transformations

**=**

*q***+ (1/2)**

*p***, spherical coordinates for**

*k***assuming $k\Vert z\u0302$ without loss of generality, and**

*p**u*= |

**+ (1/2)**

*p***|,**

*k**w*= |

**− (1/2)**

*p***| with a surface element identity**

*k**p*

^{2}sin

*θdθdp*= (

*uw*/

*k*)

*dudw*. Ultimately, the IET closure functional becomes

Utilizing dimensionless variables for all the wave-vectors, i.e., *y* → *u*/*k*_{F}, *z* → *w*/*k*_{F}, and *x* → *k*/*k*_{F}, and rearranging the integrands, the double integral expression for the IET closure functional reads as

To our knowledge, bipolar coordinates have never been utilized in earlier applications of the CA scheme^{57,58} and the HNC scheme,^{34,35} whose LFC closure functional is also formally expressed as a triple integral of the type *∫d*^{3}*q*(** k** ·

**)**

*q**I*(|

**−**

*k***|)**

*q**J*(|

**|), although this would lead to a drastic reduction in the computational cost without invoking extra approximations. This comes in stark contrast to the so-called modified versions of the CA scheme and the HNC scheme,**

*q*^{35,87}where the triple integral is converted to a single integral after replacing

*S*(

**−**

*k***) by an ad hoc Yukawa screening function of a characteristic screening wave-number that is determined by an interaction energy constraint.**

*q*### B. On the convergence of the infinite Matsubara series

Independent of the dielectric scheme (STLS, CA, qSTLS, HNC, and IET), the Matsubara summation of Eq. (5) is slowly converging, especially for small values of the degeneracy parameter *θ*. The implementation of mathematical tricks that speeds up the convergence rate is especially important for computationally heavy schemes, such as the HNC, IET, and qSTLS.

A significant speed-up of the rate of convergence can be achieved by splitting the Hartree–Fock SSF, i.e., the SSF in the absence of Coulomb pair interactions *U*(*k*) = 0, since its respective Matsubara summation can be calculated exactly.^{49,54,62,88} Introducing the auxiliary complex function $\Phi (k,z)=\u2212(2EF)/(3n)\chi \u03030(k,z)$ and normalized wave-vectors *x* → *k*/*k*_{F}, Eq. (5) becomes

where, courtesy of logarithm properties and the product $\u220fn=\u2212\u221e+\u221ea2+\pi 2n2/b2+\pi 2n2=sinh2(a)/sinh2(b)$, the Hartree–Fock SSF $SHF(x)=(3/2)\theta \u2211l=\u2212\u221e\u221e\Phi (x,l)$ is given by the following integral:^{49,54,88}

Convergence can be further accelerated by splitting the square of the high frequency–short wavelength asymptotic form (*l*, *q* → *∞*) of the auxiliary complex function $\Phi \u221e(x,l)=(4/3)x2/[x4+(2\pi l\theta )2]+O(x\u22124,l\u22124)$ since its respective Matsubara summation can also be calculated exactly.^{49,54,62} This leads to

where, by differentiating both sides of the known formula $\u2211n=\u2212\u221e+\u221e(x4+a2n2)\u22121=[\pi /(ax2)]coth\pi x2/a$ with respect to *x*, the residual SSF correction that is defined by $S\u221e(x)=(6/\pi )\lambda rs\theta (1/x2)1\u2212G(x)\u2211l=\u2212\u221e\u221e\Phi \u221e2(x,l)$ is given by^{49,54}

where coth(·) is the hyperbolic cotangent and csch(·) is the hyperbolic cosecant.

On the practical side, with an accuracy goal of 0.001% in the interaction energies, the above tricks speed up the convergence rate by more than two orders of magnitude. To be more specific, for *r*_{s} = 100 and for *θ* = 0.5 or 4, this accuracy goal is well achieved within *l* = 512 when these two splitting procedures are implemented. On the other hand, for *r*_{s} = 100, the interaction energy accuracy at *l* = 51 200 is merely 0.1% for *θ* = 4 and 0.7% for *θ* = 0.5 when no splitting is implemented.

### C. On the asymptotic convergence of the IET local field correction

From the general SSF decomposition in terms of *S*_{HF}(*x*) and *S*_{∞}(*x*), it is rather straightforward to prove that the frequency independent LFC assumption *G*(** k**,

*ω*) ≡

*G*(

**) implies the (3**

*k**π*/8)

*x*

^{4}×[1 −

*S*(

*x*)] =

*λr*

_{s}[1 −

*G*(

*x*)] asymptotic limit.

^{49}When combined with Kimball’s expression

^{89}$\u2202g(r\u0303=0)/\u2202r\u0303=(3\pi /8)limx\u2192\u221ex4[1\u2212S(x)]$ with $r\u0303=rkF$ and the cusp relation

^{89,90}$\u2202g(r\u0303=0)/\u2202r\u0303=\lambda rsg(0)$, this yields the general asymptotic condition

*G*(

*x*→

*∞*) = 1 −

*g*(0),

^{49,91}where

*g*(0) is the contact or on-top value of the radial distribution function. This asymptotic condition is often considered as a self-consistency condition since it solely originates from the first two building blocks of the dielectric formalism and needs to be independently satisfied by the third building block (closure functional). This has been shown to be valid in the case of the STLS scheme

^{49}and can also be confirmed for other dielectric schemes (VS, CA, and HNC) and for the IET scheme, given the rapid

*B*(

*x*)/

*βU*(

*x*) decay to zero. Note that in the strongly coupled electron liquid regime,

*g*(0) ≃ 0 is expected,

^{59,92}leading to

*G*(

*x*→

*∞*) ≃ 1.

Nevertheless, the implicit form of the IET LFC can strongly inhibit its proper numerical convergence at short wavelengths. In particular, the converged LFC solution has been consistently observed to exhibit a short wavelength dependence, which abruptly drops from near-unity to near-zero. The unphysical asymptotic behavior always takes place close to the wavenumber cutoff considered in the integrations, regardless of its actual value. It translates to interaction energy errors of the order of 0.005%, which exceed our 0.001% accuracy goal.

In order to achieve proper convergence of the IET local field correction in the short wavelength limit, it is highly beneficial to split the STLS LFC from the IET LFC. The STLS LFC emerges by setting *B*(*k*) ≡ 0 and *G*(*k*) ≡ 0 in the right-hand side of the IET LFC. The standard single integral form of the STLS LFC is recovered from the double integral form of the IET LFC by employing the change in variables (*y*, *z*) → (*t*, *s*) of the form $y=x2+s2\u22122xst$, *z* = *s* and carrying out the *t*− integration. This leads to

In the above equation, *G*_{1}(*x*) denotes the STLS LFC with a numerical asymptotic behavior *G*_{1}(*x* → *∞*) ≃ 1 and *G*_{2}(*x*) denotes the residual IET LFC with a numerical asymptotic behavior *G*_{2}(*x* → *∞*) = 0. Thus, the IET LFC short wavelength limit now correctly converges to a value close to unity.

### D. Structure and details of the IET algorithm

The closed normalized set of equations for the IET-based dielectric scheme comprises of the following: the normalization condition of the Fermi–Dirac energy distribution function [see Eq. (11), Sec. II A], the ideal Lindhard density response expressed through the auxiliary complex function Φ(*x*, *l*) and evaluated at the imaginary Matsubara frequencies *ω*_{l} including the static limit [see Eqs. (12) and (13), Sec. II A], the Fourier transform of the classical OCP bridge function [see Eq. (14), Sec. II C], the infinite Matsubara summation expression for the SSF after separating the Hartree–Fock SSF and also the asymptotic component [see Eq. (15), Sec. III B], and the IET LFC double integral expression after utilizing bipolar coordinates and splitting the STLS LFC single integral expression [see Eq. (16), Secs. III A and III C],

Concerning the origin of Eq. (14), within the classical *x* = *r*/*d* or *q* = *kd* normalization, it is rather straightforward to show that $B(q)/\beta U(q)=(q/\Gamma )\u222b0\u221exb(x,\Gamma )sin(qx)dx$ for the ratio of spatial Fourier transforms, where *b*(*x*, Γ) is directly adopted from Eq. (8). Naturally, the above expression needs to be translated to the quantum *x* = *rk*_{F} or *q* = *k*/*k*_{F} normalization, which is formally equivalent to the substitution *q* → *q*/*λ*. Finally, it is apparent that the utilization of classical OCP bridge functions necessitates a mapping of the quantum states (*r*_{s}, *θ*) to classical states (Γ) via Γ = 2*λ*^{2}(*r*_{s}/*θ*). The *B*(*q*)/*βU*(*q*) contribution is illustrated in Fig. 1.

The accuracy goal was set to 0.001% with respect to the interaction energies. All improper integrals were numerically evaluated with the doubly adaptive Clenshaw–Curtis quadrature method, as implemented in the GNU Scientific Library, with a 0.1*k*_{F} grid resolution and with a 40*k*_{F} upper cutoff. The only exception was the complete Fermi–Dirac integral $I1/2(\mu \u0304)$, which is implemented as a special function in the GNU Scientific Library. It should be noted that non-adaptive quadrature rules were confirmed to require much denser grid spacing <0.005*k*_{F} in order to satisfy the accuracy goal. The infinite Matsubara summation was truncated at |*l*| = 512. Convergence studies were carried out, for representative quantum coupling parameters and degeneracy parameters, in order to ensure that the chosen resolution and cutoffs do not affect the thermodynamic and structural results.

The iteration cycle proceeds as follows: (i) The reduced chemical potential $\mu \u0304$ is calculated from Eq. (11) using a bisection root-finding algorithm. (ii) The Lindhard density responses (0 ≤ *l* ≤ 512) are computed from Eqs. (12) and (13). (iii) The Fourier transform of the classical OCP bridge function is computed from Eq. (14) and stored. (iv) With the RPA’s LFC as an initial guess, the SSF is evaluated from Eq. (15). (v) These LFC and SSF are substituted in the right-hand side of Eq. (16) for an initial evaluation of the IET LFC. (vi) The last two steps are repeated until the absolute relative difference between two successive IET LFC evaluations is smaller than 10^{−5} for all the grid points. Starting from the RPA solution, convergence is typically reached within 200 iterations. Especially for state points with *r*_{s} > 100 and *θ* < 1.0, Broyles’s technique of mixing iterates^{65,93} was necessary to speed up and sometimes even to achieve convergence.

## IV. COMPUTATIONAL

To compute accurate benchmark data for our new IET scheme, we have carried out direct PIMC simulations^{94–96} of the UEF without any nodal restrictions.^{97} As a consequence, our simulations are afflicted with the notorious fermion sign problem,^{98} which, in general, leads to an exponential increase in the compute time by increasing the system-size *N* or decreasing the temperature *T*; the reader is addressed to Refs. 99 and 100 for topical and accessible review articles. In practice, however, the sign problem is not severe as quantum exchange effects are effectively reduced by the strong Coulomb repulsion in the strongly coupled electron liquid regime.

The basic idea behind the PIMC method is to evaluate the (canonical, i.e., system-size *N*, volume *V*, and inverse temperature *β* = 1/*k*_{B}*T* are fixed) partition function in coordinate space, which gives

Specifically, $R=(r1,\u2026,rN)T$ contains the coordinates of all *N* = *N*^{↑} + *N*^{↓} electrons, and due to the antisymmetry of the fermionic density matrix under the exchange of particle coordinates, we have to explicitly take the sums over all elements *σ*_{i} of the respective permutation group $SNi$, *i* ∈ {↑, ↓}, where the sign function sgn(*σ*^{↑}, *σ*^{↓}) gives positive (negative) unity for an even (odd) number of pair permutations. Furthermore, the operators $\pi \u0302\sigma i$ realize the particular permutations for a corresponding element *σ*^{i}.

The problem with Eq. (17) is that the matrix elements of the density operator $e\u2212\beta H\u0302$ cannot be readily evaluated since the kinetic ($K\u0302$) and potential ($V\u0302$) contributions to the Hamiltonian $H\u0302=K\u0302+V\u0302$ do not commute, $e\u2212\beta H\u0302\u2260e\u2212\beta K\u0302e\u2212\beta V\u0302$. In order to overcome this obstacle, we utilize the exact semi-group property of $\rho \u0302$,

where the definition *ϵ* = *β*/*P* has been employed. Applying Eq. (18) to Eq. (17) and inserting *P* − 1 unity operators of the form $1\u0302=\u222bdR\alpha |R\alpha \u3009\u3008R\alpha |$ lead to the intermediate result

which is still exact. Evidently, Eq. (19) requires the evaluation of *P* density matrices, but at *P* times the temperature. For a sufficiently large *P*, each of these factors can be straightforwardly evaluated using a suitable high-temperature approximation, such as the primitive factorization $e\u2212\u03f5H\u0302\u2248e\u2212\u03f5K\u0302e\u2212\u03f5V\u0302$. In fact, the factorization error of the latter decays as ∼1/*P*^{2},^{101} and the convergence in the limit of *P* → *∞* is ensured by the celebrated Trotter formula,^{102}

In practice, we find *P* = 200 sufficient to reduce the factorization error substantially below the noise level of the Monte Carlo simulation. For completeness, we note that Eq. (20) only holds for the case of potentials $V\u0302$ that are bounded from below, which is, indeed, the case for the UEF. Attractive potentials such as the Coulomb interaction between positive and negative charges require a modified procedure such as the pair approximation, which is discussed in detail, for instance, in Ref. 103. In addition, we note that higher-order factorizations of the thermal density matrix that converge as ∼1/*P*^{4}^{104,105} and even ∼1/*P*^{6−8}^{106} have been presented, although we do not find them necessary for the present conditions.

The final result for the PIMC partition function can then be written in the abbreviated form as

where the integration over the $X=(R0,\u2026,RP\u22121)T$ meta-variable also contains the summation over all the possible permutations. Furthermore, the weight function *W*(**X**) can be readily evaluated for each individual configuration **X** and contains contributions from the potential energy and from the free thermal density matrix; see, for instance, Ref. 94 for an accessible review article. Evidently, Eq. (21) requires the evaluation of a *d* = 3*PN*-dimensional integral, which easily leads to *d* ∼ 10^{4} in the present study. While this renders the application of standard quadrature methods impractical due to the well-known curse of dimensionality, the Metropolis Monte Carlo method^{107} does not suffer from this drawback.

A graphical illustration of the PIMC approach is shown in Fig. 2. Specifically, we show a configuration **X** of *N* = 4 electrons with *P* = 6 high-temperature factors. It is evident that each particle is represented by a closed (i.e., *β*-periodic) path of particle coordinates in the imaginary time *τ* ∈ [0, *β*] (strictly speaking, it is *τ* ∈ *ℏ*/*i*[0, *β*], but it is conventional to drop the pre-factor for simplicity). This, in turn, corresponds to the famous *classical isomorphism*,^{108} where the complicated quantum many-body system of interest is mapped onto an effective classical system of interacting ring-polymers.

An additional difficulty arises due to the indistinguishable nature of electrons of the same spin-species, which requires us to sample all the possible permutations of particle coordinates. Within the PIMC picture, the latter manifest as the so-called permutation cycles,^{109} which are trajectories comprising more than a single particle. An example for such a permutation cycle can be identified at the center of Fig. 2. As a consequence of the depicted pair exchange, we have to move twice through the imaginary time to return to the point of origin. This, in turn, results in a negative sign of the corresponding configuration weight, *W*(**X**) < 0, and thereby contributes to the aforementioned fermion sign problem.

In practice, we construct a Markov chain of random configurations {**X**_{i}} using a canonical adaption^{110} of the worm algorithm by Boninsegni *et al*.^{111,112}

Let us next consider the PIMC estimation of the interaction energy per particle *V*_{N}/*N*, which is shown in Fig. 3 for *r*_{s} = 125 and *θ* = 2. In particular, PIMC simulations are, by default, only possible for a finite number of electrons *N*. The raw PIMC data for three different system sizes are shown as the red circles in Fig. 3 and exhibit a significant dependence on *N*. In practice, however, we are interested in the thermodynamic limit, i.e., in the limit of an infinite number of particles with the density being constant,

To eliminate the difference between the PIMC data for *V*_{N}/*N* and Eq. (22), we use the finite-size correction by Chiesa *et al.*,^{113} which has subsequently been adapted to finite temperatures in Ref. 114,

where the plasma frequency is given by $\omega p=3/rs3$ in Hartree atomic units. To be more specific, Eq. (23) is based on the insight that the system-size dependence of *V*_{N}/*N* is mainly the consequence of the approximation of a continuous integral [see Eq. (25)] by the sum over reciprocal lattice vectors due to the momentum quantization in a finite simulation cell. To the first order, this discretization error can be approximated by utilizing the exact long wave-length limit of the static structure factor^{115}

A detailed derivation is beyond the scope of the present work and was presented by Drummond *et al.*^{116} Adding the finite-size correction given in Eq. (23) to the PIMC data leads to the green crosses in Fig. 3. Evidently, the bulk of the finite-size errors have been removed, and the corrected data points fall into an interval of 0.01% (horizontal light gray lines) around their common average value (horizontal blue line). All the PIMC results for the interaction energy that are shown in this work have been obtained by following this procedure.

For completeness, it should be mentioned that finite-size effects are substantially more pronounced in the warm dense matter regime and that the simple first-order correction from Eq. (23) breaks down.^{19,117} A recent investigation of finite-size effects of a uniform electron gas at extreme densities and temperatures was presented by Dornheim and Vorberger.^{118}

## V. RESULTS

Here, we compare the paramagnetic electron liquid interaction energies and static properties as computed from different dielectric formalism schemes with their “exact” counterparts as extracted from our PIMC simulations. It has been demonstrated that beyond the warm dense matter regime and especially for *r*_{s} > 30, the STLS and the VS schemes yield increasingly inaccurate results.^{34} Only the IET, HNC, and qSTLS schemes will be numerically solved herein since the inclusion of the STLS and VS schemes has been judged to be meaningless. The HNC and IET comparison will lead to direct conclusions regarding the impact of the classical OCP bridge function inclusion, while the qSTLS and IET comparison will lead to indirect conclusions regarding the significance of the beyond-RPA quantum effects. Our HNC algorithm can be essentially obtained from the IET algorithm described in Sec. III D by setting *B*(*q*)/[*βU*(*q*)] ≡ 0 and has been systematically validated against the published results of Tanaka.^{35} On the other hand, our qSTLS algorithm has marked differences from the IET algorithm described in Sec. III D owing to the dynamic nature of the LFC and has been benchmarked against the published results of Schweng and Böhm.^{62}

### A. Interaction energy

The interaction energy per particle is generally obtained by $U=(1/2)\u222b[d3\u2061k/(2\pi )3]U(k)S(k)\u22121$. Substituting for the Coulomb potential energy, introducing the (*r*_{s}, *θ*) variables, and employing normalized wavenumbers lead to the standard UEF expression^{19,43}

where $u\u0303(rs,\theta )$ denotes the interaction energy per particle normalized by the Hartree energy *E*_{h} = *e*^{2}/*a*_{B}.

The normalized interaction energies extracted from the PIMC simulations and computed with the qSTLS, HNC, and IET schemes are listed in Table I. The absolute relative deviations between the theoretical and “exact” interaction energies are also tabulated therein. (i) The qSTLS interaction energies are relatively inaccurate with 5.0%–8.7% relative deviations from PIMC results. The errors systematically decrease as *θ* increases and increase as *r*_{s} increases. (ii) The HNC interaction energies are very accurate with 0.47%–1.37% relative deviations from PIMC results. The errors systematically decrease as *θ* increases (diminishing quantum effects) and increase as *r*_{s} increases (stronger beyond-HNC classical pair correlations^{65–67}). (iii) The IET interaction energies are revealed to be the most accurate with merely 0.05%–0.68% relative deviations from the PIMC results. No systematic tendencies have been observed in the errors with respect to either *θ* or *r*_{s}. (iv) The IET scheme provides the most accurate interaction energy predictions for all 20 thermodynamic state points investigated, with appreciable improvements over the HNC interaction energy predictions. (v) Despite their very high accuracy, the IET (as well as the HNC) interaction energies cannot reproduce the exact *θ*-dependence of $u\u0303(rs,\theta )$ within the highly degenerate *θ* ≲ 1 range. To be more specific, the PIMC results reveal a monotonic $|u\u0303|$ decrease as *θ* increases, while the IET (and HNC) results exhibit a monotonic $|u\u0303|$ increase as *θ* increases within *θ* ≲ 1 and the correct monotonic $|u\u0303|$ decrease only as *θ* increases within *θ* ≳ 1. For this reason, we did not construct an exchange–correlation free energy parameterization via the adiabatic connection formula.

r_{s}
. | θ
. | Γ . | $u\u0303PIMC$ . | $u\u0303qSTLS$ . | e_{qSTLS} (%)
. | $u\u0303HNC$ . | e_{HNC} (%)
. | $u\u0303IET$ . | e_{IET} (%)
. |
---|---|---|---|---|---|---|---|---|---|

100 | 0.50 | 108.6 | −0.008 255 00 | −0.007 623 82 | 7.646 | −0.008 158 66 | 1.167 | −0.008 221 81 | 0.402 |

100 | 0.75 | 72.40 | −0.008 245 70 | −0.007 647 55 | 7.254 | −0.008 164 90 | 0.980 | −0.008 225 44 | 0.246 |

100 | 1.00 | 54.30 | −0.008 234 90 | −0.007 667 12 | 6.895 | −0.008 166 18 | 0.834 | −0.008 225 59 | 0.113 |

100 | 2.00 | 27.15 | −0.008 176 50 | −0.007 688 92 | 5.963 | −0.008 129 05 | 0.580 | −0.008 190 66 | 0.173 |

100 | 4.00 | 13.58 | −0.008 006 23 | −0.007 603 08 | 5.035 | −0.007 968 33 | 0.473 | −0.008 031 43 | 0.315 |

50 | 0.50 | 54.30 | −0.016 007 00 | −0.014 998 07 | 6.303 | −0.015 898 41 | 0.678 | −0.016 035 10 | 0.176 |

60 | 0.50 | 65.16 | −0.013 453 10 | −0.012 560 41 | 6.636 | −0.013 348 04 | 0.781 | −0.013 460 14 | 0.052 |

70 | 0.50 | 76.02 | −0.011 611 75 | −0.010 807 32 | 6.928 | −0.011 509 38 | 0.882 | −0.011 603 90 | 0.068 |

80 | 0.50 | 86.88 | −0.010 219 37 | −0.009 485 45 | 7.182 | −0.010 120 12 | 0.971 | −0.010 201 49 | 0.175 |

90 | 0.50 | 97.74 | −0.009 128 62 | −0.008 452 87 | 7.403 | −0.009 032 93 | 1.048 | −0.009 104 15 | 0.268 |

110 | 0.50 | 119.5 | −0.007 526 42 | −0.006 943 43 | 7.746 | −0.007 440 12 | 1.147 | −0.007 496 75 | 0.394 |

125 | 0.50 | 135.8 | −0.006 654 21 | −0.006 124 35 | 7.963 | −0.006 573 77 | 1.209 | −0.006 622 68 | 0.474 |

125 | 0.75 | 90.50 | −0.006 650 53 | −0.006 143 35 | 7.626 | −0.006 578 38 | 1.085 | −0.006 625 56 | 0.442 |

125 | 1.00 | 67.88 | −0.006 643 36 | −0.006 159 66 | 7.281 | −0.006 579 99 | 0.954 | −0.006 626 47 | 0.254 |

125 | 1.50 | 45.25 | −0.006 625 35 | −0.006 179 40 | 6.731 | −0.006 574 32 | 0.770 | −0.006 621 12 | 0.064 |

125 | 2.00 | 33.94 | −0.006 602 98 | −0.006 184 76 | 6.334 | −0.006 559 00 | 0.666 | −0.006 607 12 | 0.063 |

150 | 0.50 | 162.9 | −0.005 581 77 | −0.005 119 15 | 8.288 | −0.005 508 21 | 1.318 | −0.005 547 97 | 0.606 |

150 | 1.00 | 81.45 | −0.005 571 34 | −0.005 148 88 | 7.583 | −0.005 513 37 | 1.040 | −0.005 551 32 | 0.359 |

200 | 0.50 | 217.2 | −0.004 222 44 | −0.003 855 65 | 8.687 | −0.004 164 45 | 1.373 | −0.004 193 73 | 0.680 |

200 | 1.00 | 108.6 | −0.004 217 10 | −0.003 878 21 | 8.036 | −0.004 168 13 | 1.161 | −0.004 195 59 | 0.510 |

r_{s}
. | θ
. | Γ . | $u\u0303PIMC$ . | $u\u0303qSTLS$ . | e_{qSTLS} (%)
. | $u\u0303HNC$ . | e_{HNC} (%)
. | $u\u0303IET$ . | e_{IET} (%)
. |
---|---|---|---|---|---|---|---|---|---|

100 | 0.50 | 108.6 | −0.008 255 00 | −0.007 623 82 | 7.646 | −0.008 158 66 | 1.167 | −0.008 221 81 | 0.402 |

100 | 0.75 | 72.40 | −0.008 245 70 | −0.007 647 55 | 7.254 | −0.008 164 90 | 0.980 | −0.008 225 44 | 0.246 |

100 | 1.00 | 54.30 | −0.008 234 90 | −0.007 667 12 | 6.895 | −0.008 166 18 | 0.834 | −0.008 225 59 | 0.113 |

100 | 2.00 | 27.15 | −0.008 176 50 | −0.007 688 92 | 5.963 | −0.008 129 05 | 0.580 | −0.008 190 66 | 0.173 |

100 | 4.00 | 13.58 | −0.008 006 23 | −0.007 603 08 | 5.035 | −0.007 968 33 | 0.473 | −0.008 031 43 | 0.315 |

50 | 0.50 | 54.30 | −0.016 007 00 | −0.014 998 07 | 6.303 | −0.015 898 41 | 0.678 | −0.016 035 10 | 0.176 |

60 | 0.50 | 65.16 | −0.013 453 10 | −0.012 560 41 | 6.636 | −0.013 348 04 | 0.781 | −0.013 460 14 | 0.052 |

70 | 0.50 | 76.02 | −0.011 611 75 | −0.010 807 32 | 6.928 | −0.011 509 38 | 0.882 | −0.011 603 90 | 0.068 |

80 | 0.50 | 86.88 | −0.010 219 37 | −0.009 485 45 | 7.182 | −0.010 120 12 | 0.971 | −0.010 201 49 | 0.175 |

90 | 0.50 | 97.74 | −0.009 128 62 | −0.008 452 87 | 7.403 | −0.009 032 93 | 1.048 | −0.009 104 15 | 0.268 |

110 | 0.50 | 119.5 | −0.007 526 42 | −0.006 943 43 | 7.746 | −0.007 440 12 | 1.147 | −0.007 496 75 | 0.394 |

125 | 0.50 | 135.8 | −0.006 654 21 | −0.006 124 35 | 7.963 | −0.006 573 77 | 1.209 | −0.006 622 68 | 0.474 |

125 | 0.75 | 90.50 | −0.006 650 53 | −0.006 143 35 | 7.626 | −0.006 578 38 | 1.085 | −0.006 625 56 | 0.442 |

125 | 1.00 | 67.88 | −0.006 643 36 | −0.006 159 66 | 7.281 | −0.006 579 99 | 0.954 | −0.006 626 47 | 0.254 |

125 | 1.50 | 45.25 | −0.006 625 35 | −0.006 179 40 | 6.731 | −0.006 574 32 | 0.770 | −0.006 621 12 | 0.064 |

125 | 2.00 | 33.94 | −0.006 602 98 | −0.006 184 76 | 6.334 | −0.006 559 00 | 0.666 | −0.006 607 12 | 0.063 |

150 | 0.50 | 162.9 | −0.005 581 77 | −0.005 119 15 | 8.288 | −0.005 508 21 | 1.318 | −0.005 547 97 | 0.606 |

150 | 1.00 | 81.45 | −0.005 571 34 | −0.005 148 88 | 7.583 | −0.005 513 37 | 1.040 | −0.005 551 32 | 0.359 |

200 | 0.50 | 217.2 | −0.004 222 44 | −0.003 855 65 | 8.687 | −0.004 164 45 | 1.373 | −0.004 193 73 | 0.680 |

200 | 1.00 | 108.6 | −0.004 217 10 | −0.003 878 21 | 8.036 | −0.004 168 13 | 1.161 | −0.004 195 59 | 0.510 |

It is important to point out that the employed classical OCP bridge function parameterization is strictly valid for 10 ≤ Γ ≤ 170.^{39,75} This upper threshold is surpassed by the state point (*r*_{s}, *θ*) = (200, 0.50), which corresponds to Γ = 217.2. Since all the *s*_{i}(Γ) and *l*_{i}(Γ) coefficients involved in Eq. (8) are monotonic functions of Γ, it can be expected from the well-known continuity between the stable and metastable liquid states that the analytic classical OCP bridge function expression can be extrapolated without large errors toward the supercooled liquid regime Γ > 171.8. The same applies for extrapolations in the weakly interacting regime Γ < 10, but these are less significant due to the diminished bridge function impact.

Furthermore, we compare with the interaction energies as computed from three accurate parameterizations of the exchange–correlation free energy through the expression

which is acquired after the differentiation of the thermodynamic formula $f\u0303xc(rs,\Theta )=rs\u22122\u222b0rsrs\u2032u\u0303int(rs\u2032,\Theta )drs\u2032$.^{19} In particular, we consider the parameterization by Groth *et al.* (GDSMFB) that is based on simulation results obtained by various novel PIMC methods within 0.1 ≤ *r*_{s} ≤ 20 and 0.5 ≤ *θ* ≤ 8,^{119} the parameterization by Karasiev *et al.* (KSDT) that is based on simulation results obtained by the restricted PIMC method within 1.0 ≤ *r*_{s} ≤ 40 and 0.062 5 ≤ *θ* ≤ 8,^{120} and the corrected version of the parameterization by Karasiev *et al.* (corrKSDT) that is also based on the above restricted PIMC data.^{121} Despite some documented deficiencies of the restricted PIMC data input^{114} (concerning the uncontrolled fixed node approximation^{122} and the unsatisfactory treatment of finite-size effects^{117}) and despite a procedural mistake in the original fitting procedure (concerning the utilization of an analytic ground-state fit instead of the actual ground-state MC data^{123}), the KSDT interaction energies exhibit 0.39%–1.63% relative deviations from PIMC results, whereas the GDSMFB interaction energies exhibit 1.08%–4.94% relative deviations and the corrKSDT interaction energies exhibit 1.58%–6.90% relative deviations from PIMC results. Therefore, extrapolated GDSMFB, KSDT, and corrKSDT interaction energies are more accurate than the qSTLS results but less accurate than the IET results. More specifically, the KSDT interaction energies are even more accurate than the HNC results at some investigated states but never more accurate than the IET results. Hence, even though the GDSMFB and corrKSDT exchange–correlation free energy parameterization is much more accurate than the KSDT parameterization within the warm dense matter regime,^{119} it can be concluded that the latter extrapolates better within the strongly coupled electron regime, at least as far as interaction energies are concerned. For completeness, we point out that the GDSMFB, KSDT, and corrKSDT discrepancies in warm dense matter ranges do not impact density functional theory calculations.^{124}

Finally, we also compare with the interaction energies as computed from the direct parameterization of the interaction energy by Ichimaru *et al.* (IIT), which is based on the STLS interaction energies after their correction, to comply with variable coupling QMC simulations for the ground-state limit (*θ* → 0) and with variable coupling MC simulations for the classical limit (*θ* → *∞*), via the implementation of an ad hoc *θ*− interpolation function.^{43,125,126} Although the interpolation function accuracy is unclear for intermediate degeneracies,^{19} the incorporation of exact strong coupling ground-state results^{127} and classical results^{128,129} hints that the IIT parameterization might be very accurate. This is verified by the comparison that reveals that the IIT interaction energies exhibit 0.06%–0.63% relative deviations from PIMC results. Remarkably, IIT interaction energies are much more accurate than the qSTLS, more accurate than the HNC, and even as accurate as the IET results. In particular, the IET interaction energies are more accurate for 11 and the IIT interaction energies more accurate for nine of the studied states, while the mean absolute relative errors are 0.29% for the IET and 0.35% for the IIT.

### B. Static structure factor

Characteristic static structure factors extracted from our PIMC simulations and computed with the three dielectric formalism schemes are illustrated in Figs. 4 and 5. For all the 20 investigated state points, the magnitudes and positions of the first SSF peak resulting from the PIMC simulations and from the qSTLS, HNC, and IET schemes are listed in Table II. The absolute relative deviations between the theoretical and the “exact” SSF peak values are also tabulated therein.

r_{s}
. | θ
. | $SPIMCmax$ . | $SqSTLSmax$ . | e_{qSTLS} (%)
. | $SHNCmax$ . | e_{HNC} (%)
. | $SIETmax$ . | e_{IET} (%)
. | $argqSPIMCmax$ . | $argqSqSTLSmax$ . | e_{qSTLS} (%)
. | $argqSHNCmax$ . | e_{HNC} (%)
. | $argqSIETmax$ . | e_{IET} (%)
. |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

100 | 0.50 | 1.494 | 1.224 | 18.074 | 1.202 | 19.546 | 1.306 | 12.598 | 2.248 | 1.770 | 21.247 | 2.220 | 1.225 | 2.200 | 2.115 |

100 | 0.75 | 1.449 | 1.191 | 17.779 | 1.203 | 16.932 | 1.301 | 10.159 | 2.248 | 1.770 | 21.247 | 2.220 | 1.225 | 2.210 | 1.670 |

100 | 1.00 | 1.406 | 1.166 | 17.127 | 1.199 | 14.724 | 1.290 | 8.285 | 2.191 | 1.790 | 18.288 | 2.230 | 1.797 | 2.220 | 1.341 |

100 | 2.00 | 1.285 | 1.102 | 14.204 | 1.164 | 9.366 | 1.227 | 4.461 | 2.248 | 1.870 | 16.798 | 2.270 | 0.999 | 2.250 | 0.110 |

100 | 4.00 | 1.166 | 1.051 | 9.912 | 1.105 | 5.263 | 1.138 | 2.381 | 2.303 | 2.040 | 11.422 | 2.340 | 1.605 | 2.320 | 0.736 |

50 | 0.50 | 1.252 | 1.107 | 11.569 | 1.093 | 12.677 | 1.136 | 9.211 | 2.248 | 1.910 | 15.018 | 2.290 | 1.889 | 2.250 | 0.110 |

60 | 0.50 | 1.327 | 1.131 | 14.767 | 1.115 | 15.956 | 1.171 | 11.802 | 2.231 | 1.860 | 16.635 | 2.260 | 1.294 | 2.230 | 0.051 |

70 | 0.50 | 1.351 | 1.156 | 14.450 | 1.137 | 15.822 | 1.205 | 10.861 | 2.188 | 1.830 | 16.355 | 2.250 | 2.843 | 2.220 | 1.471 |

80 | 0.50 | 1.391 | 1.180 | 15.174 | 1.159 | 16.674 | 1.238 | 10.990 | 2.144 | 1.800 | 16.029 | 2.230 | 4.030 | 2.210 | 3.097 |

90 | 0.50 | 1.438 | 1.202 | 16.418 | 1.181 | 17.905 | 1.272 | 11.551 | 2.231 | 1.790 | 19.772 | 2.220 | 0.499 | 2.210 | 0.947 |

110 | 0.50 | 1.481 | 1.248 | 15.766 | 1.223 | 17.456 | 1.339 | 9.591 | 2.144 | 1.750 | 18.362 | 2.210 | 3.097 | 2.200 | 2.631 |

125 | 0.50 | 1.538 | 1.288 | 16.310 | 1.254 | 18.502 | 1.390 | 9.661 | 2.188 | 1.720 | 21.383 | 2.200 | 0.557 | 2.200 | 0.557 |

125 | 0.75 | 1.557 | 1.243 | 20.160 | 1.256 | 19.369 | 1.384 | 11.132 | 2.231 | 1.730 | 22.461 | 2.210 | 0.947 | 2.200 | 1.395 |

125 | 1.00 | 1.501 | 1.209 | 19.455 | 1.251 | 16.671 | 1.369 | 8.765 | 2.188 | 1.750 | 20.011 | 2.220 | 1.471 | 2.210 | 1.014 |

125 | 1.50 | 1.436 | 1.163 | 19.000 | 1.231 | 14.258 | 1.331 | 7.338 | 2.231 | 1.790 | 19.772 | 2.230 | 0.051 | 2.220 | 0.499 |

125 | 2.00 | 1.365 | 1.133 | 17.039 | 1.209 | 11.464 | 1.292 | 5.360 | 2.274 | 1.820 | 19.952 | 2.250 | 1.040 | 2.230 | 1.919 |

150 | 0.50 | 1.695 | 1.355 | 20.022 | 1.305 | 23.018 | 1.474 | 13.022 | 2.144 | 1.700 | 20.694 | 2.200 | 2.631 | 2.200 | 2.631 |

150 | 1.00 | 1.586 | 1.255 | 20.916 | 1.300 | 18.038 | 1.447 | 8.768 | 2.231 | 1.720 | 22.909 | 2.210 | 0.947 | 2.210 | 0.947 |

200 | 0.50 | 1.817 | 1.436 | 20.962 | 1.403 | 22.807 | 1.638 | 9.833 | 2.231 | 1.700 | 23.806 | 2.190 | 1.844 | 2.200 | 1.395 |

200 | 1.00 | 1.738 | 1.337 | 23.112 | 1.395 | 19.757 | 1.600 | 7.943 | 2.231 | 1.700 | 23.806 | 2.200 | 1.395 | 2.200 | 1.395 |

r_{s}
. | θ
. | $SPIMCmax$ . | $SqSTLSmax$ . | e_{qSTLS} (%)
. | $SHNCmax$ . | e_{HNC} (%)
. | $SIETmax$ . | e_{IET} (%)
. | $argqSPIMCmax$ . | $argqSqSTLSmax$ . | e_{qSTLS} (%)
. | $argqSHNCmax$ . | e_{HNC} (%)
. | $argqSIETmax$ . | e_{IET} (%)
. |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

100 | 0.50 | 1.494 | 1.224 | 18.074 | 1.202 | 19.546 | 1.306 | 12.598 | 2.248 | 1.770 | 21.247 | 2.220 | 1.225 | 2.200 | 2.115 |

100 | 0.75 | 1.449 | 1.191 | 17.779 | 1.203 | 16.932 | 1.301 | 10.159 | 2.248 | 1.770 | 21.247 | 2.220 | 1.225 | 2.210 | 1.670 |

100 | 1.00 | 1.406 | 1.166 | 17.127 | 1.199 | 14.724 | 1.290 | 8.285 | 2.191 | 1.790 | 18.288 | 2.230 | 1.797 | 2.220 | 1.341 |

100 | 2.00 | 1.285 | 1.102 | 14.204 | 1.164 | 9.366 | 1.227 | 4.461 | 2.248 | 1.870 | 16.798 | 2.270 | 0.999 | 2.250 | 0.110 |

100 | 4.00 | 1.166 | 1.051 | 9.912 | 1.105 | 5.263 | 1.138 | 2.381 | 2.303 | 2.040 | 11.422 | 2.340 | 1.605 | 2.320 | 0.736 |

50 | 0.50 | 1.252 | 1.107 | 11.569 | 1.093 | 12.677 | 1.136 | 9.211 | 2.248 | 1.910 | 15.018 | 2.290 | 1.889 | 2.250 | 0.110 |

60 | 0.50 | 1.327 | 1.131 | 14.767 | 1.115 | 15.956 | 1.171 | 11.802 | 2.231 | 1.860 | 16.635 | 2.260 | 1.294 | 2.230 | 0.051 |

70 | 0.50 | 1.351 | 1.156 | 14.450 | 1.137 | 15.822 | 1.205 | 10.861 | 2.188 | 1.830 | 16.355 | 2.250 | 2.843 | 2.220 | 1.471 |

80 | 0.50 | 1.391 | 1.180 | 15.174 | 1.159 | 16.674 | 1.238 | 10.990 | 2.144 | 1.800 | 16.029 | 2.230 | 4.030 | 2.210 | 3.097 |

90 | 0.50 | 1.438 | 1.202 | 16.418 | 1.181 | 17.905 | 1.272 | 11.551 | 2.231 | 1.790 | 19.772 | 2.220 | 0.499 | 2.210 | 0.947 |

110 | 0.50 | 1.481 | 1.248 | 15.766 | 1.223 | 17.456 | 1.339 | 9.591 | 2.144 | 1.750 | 18.362 | 2.210 | 3.097 | 2.200 | 2.631 |

125 | 0.50 | 1.538 | 1.288 | 16.310 | 1.254 | 18.502 | 1.390 | 9.661 | 2.188 | 1.720 | 21.383 | 2.200 | 0.557 | 2.200 | 0.557 |

125 | 0.75 | 1.557 | 1.243 | 20.160 | 1.256 | 19.369 | 1.384 | 11.132 | 2.231 | 1.730 | 22.461 | 2.210 | 0.947 | 2.200 | 1.395 |

125 | 1.00 | 1.501 | 1.209 | 19.455 | 1.251 | 16.671 | 1.369 | 8.765 | 2.188 | 1.750 | 20.011 | 2.220 | 1.471 | 2.210 | 1.014 |

125 | 1.50 | 1.436 | 1.163 | 19.000 | 1.231 | 14.258 | 1.331 | 7.338 | 2.231 | 1.790 | 19.772 | 2.230 | 0.051 | 2.220 | 0.499 |

125 | 2.00 | 1.365 | 1.133 | 17.039 | 1.209 | 11.464 | 1.292 | 5.360 | 2.274 | 1.820 | 19.952 | 2.250 | 1.040 | 2.230 | 1.919 |

150 | 0.50 | 1.695 | 1.355 | 20.022 | 1.305 | 23.018 | 1.474 | 13.022 | 2.144 | 1.700 | 20.694 | 2.200 | 2.631 | 2.200 | 2.631 |

150 | 1.00 | 1.586 | 1.255 | 20.916 | 1.300 | 18.038 | 1.447 | 8.768 | 2.231 | 1.720 | 22.909 | 2.210 | 0.947 | 2.210 | 0.947 |

200 | 0.50 | 1.817 | 1.436 | 20.962 | 1.403 | 22.807 | 1.638 | 9.833 | 2.231 | 1.700 | 23.806 | 2.190 | 1.844 | 2.200 | 1.395 |

200 | 1.00 | 1.738 | 1.337 | 23.112 | 1.395 | 19.757 | 1.600 | 7.943 | 2.231 | 1.700 | 23.806 | 2.200 | 1.395 | 2.200 | 1.395 |

For the state points of interest, the liquid character of the UEF becomes apparent from the relatively large magnitude of the first SSF peak (which is generally around ∼1.5 and even reaches 1.82), the well-resolved first SSF trough around 3*k*_{F} (especially for *r*_{s} ≳ 100), and the well-resolved second SSF peak above 4*k*_{F} (only for *r*_{s} ≳ 100). (i) The IET scheme always generates the most accurate SSF across the entire interval: within the long wavelength range of *k* ≲ 2*k*_{F}, in the Lorentzian shaped region that surrounds the first maximum *k* ∼ 2 − 2.5*k*_{F}, and within the short wavelength range of *k* ≳ 2.5*k*_{F}. (ii) The qSTLS scheme has the worst performance. It strongly underestimates the position of the first SSF peak by around 20%. The same behavior was earlier observed for the STLS and the VS scheme.^{34} These results can be anticipated from classical OCP liquids, where it has been established that stronger Coulomb correlations (neglected in STLS, VS, and qSTLS) not only increase the SSF peak magnitude but also slightly displace it toward longer wavenumbers. (iii) On the other hand, the IET and HNC schemes provide very accurate predictions for the SSF peak positions, 0.05%–3.10% and 0.05%–4.03% relative deviations from the PIMC results, respectively, with the IET scheme having the slight edge. Thus, it can be concluded that the first SSF peak position is mainly controlled by strong correlations. (iv) Regardless of the state, the IET scheme greatly improves the HNC prediction for the SSF peak magnitude, with the relative deviations from PIMC results being 2.38%–13.02% and 5.26%–23.02%, respectively. (v) For all states, the IET SSF is characterized by a very accurate long wavelength behavior, especially for *k* ≲ *k*_{F}.

Despite the large improvements over the HNC scheme, the IET scheme seems to generate SSFs that are not accurate enough to justify the very accurate predictions of interaction energies. Detailed inspection of Figs. 4 and 5 reveals that the very accurate interaction energies are the result of favorable error cancellation in Eq. (25) since the IET SSF tends to be slightly too large within *k*_{F} ≤ *k* ≤ 2*k*_{F}, becomes clearly too small within 2.0*k*_{F} ≤ *k* ≤ 2.5*k*_{F}, and tends to be somewhat too large within 2.5*k*_{F} ≤ *k* ≤ 4.0*k*_{F}. The same reasoning applies for the HNC scheme. Note that a similar favorable error cancellation is responsible for the success of STLS generated interaction energies within the warm dense matter regime.^{19,126}

Let us briefly discuss the radial distribution functions (RDFs) *g*(*r*) extracted from PIMC simulations and computed with the three dielectric formalism schemes. Since the RDFs encode the same correlation information as SSFs, it is unsurprising that the IET scheme greatly improves the predictions of the HNC scheme for the first maximum and first non-zero minimum magnitude, that the IET and HNC schemes provide very accurate similar predictions for the first maximum and first non-zero minimum positions, and that the qSTLS scheme has the worst performance. Finally, it is important to discuss the un-physical negative RDF region near the origin, a common pathology of all dielectric formalism schemes stemming from the approximate treatment of quantum effects.^{19} The IET and HNC schemes have negative regions of similar extent and signal power, as expected from the RPA-level treatment of quantum effects. On the other hand, the qSTLS scheme has a negative region of less extent and power, as expected from its beyond-RPA quantum features.

### C. Static local field correction and static density response

Characteristic static local field corrections and static density responses extracted from our PIMC simulations and computed with the three dielectric formalism schemes are illustrated in Figs. 6 and 7.

Let us first discuss the *static local field correction*: (i) The PIMC LFC exhibits a pronounced first local maximum of magnitude ∼1.10–1.25 located at *k* ∼ 2.35–2.60*k*_{F}, which is followed by a shallow local minimum located at *k* ∼ 3.05–3.45*k*_{F} and is ultimately succeeded by a rather sharp short wavelength increase. The asymptotic behavior does not contradict the short wavelength limit *G*(*k* → *∞*) = 1 − *g*(0)^{49,91} that was discussed earlier since this is valid for frequency independent LFCs (as confirmed by the HNC and the IET asymptotics), but not for the exact static LFC. In the ground state *θ* → 0, the asymptotic behavior of the exact static LFC is described by the Holas expansion *G*(*k* → *∞*, 0) = *B* + *Ck*^{2}, which predicts a parabolic divergence.^{130} Such a parabolic asymptote has been empirically observed to persist also when *θ* ∼ 1,^{34,131} but not in the classical limit *θ* ≫ 1, where one gets *G*(*k* → *∞*, 0) = 1.^{43} The above observation is confirmed by the present PIMC LFC since as *θ* increases, the asymptotic limit gradually switches from parabolically diverging to nearly unity in the depicted wavenumber range; see Figs. 7(a)–7(d). (ii) Regardless of the state point of interest, the IET scheme generates the most accurate static LFC and the qSTLS scheme generates the least accurate. The IET scheme improves the HNC prediction for the LFC first peak magnitude, with the relative deviations from PIMC results being 1.4%–10.1% and 2.6%–12.0%, respectively. However, the HNC prediction for the LFC first peak position is always more accurate. (iii) Perhaps, the most remarkable feature of the IET scheme’s LFC concerns its very accurate long wavelength behavior that extends up to nearly *k* ≲ 2*k*_{F}. This does not necessarily guarantee that the IET scheme satisfies to a high degree the compressibility sum rule (CSR), an exact relation connecting the long wavelength static LFC to the second density derivative of the exchange–correlation free energy, i.e., in (*n*, *T*) thermodynamic variables and cgs units,^{19,43,51,52}

or in (*r*_{s}, *θ*) thermodynamic variables and Hartree units,

This is simply because each dielectric scheme satisfies its individual CSR that involves its individual exchange–correlation free energy and not the exact exchange–correlation free energy. However, when combined with the very accurate IET interaction energies, the very accurate IET LFC long wavelength limit makes it very likely that the IET scheme satisfies the CSR to a large degree. We postpone such an investigation to a future work owing to its high computational cost since as Eq. (28) suggests, it requires an accurate parameterization of the IET exchange–correlation free energy. Nevertheless, this IET feature can be exploited to provide an alternative route for the accurate calculation of the paramagnetic electron liquid’s isothermal compressibility without involving neither PIMC simulations nor thermodynamic integration.

Let us now discuss the *static density response function* defined by *χ*(** k**) ≡

*χ*(

**,**

*k**ω*= 0). (i) The qualitative behavior of the “exact”

*χ*(

**) of the paramagnetic electron liquid has been discussed earlier**

*k*^{34}and is confirmed by our PIMC results. It is evident that the magnitude of the pronounced

*χ*(

**) extremum increases as**

*k**r*

_{s}decreases (see Fig. 6) and as

*θ*decreases (see Fig. 7) and that the width of the

*χ*(

**) extremum becomes more sharp as**

*k**r*

_{s}increases (see Fig. 6) and especially as

*θ*decreases (see Fig. 7). (ii) Since

*χ*(

**) is obtained by setting**

*k**ω*= 0 to Eq. (4) and only involves the common among schemes

*χ*

_{0}(

**) and the varying among schemes**

*k**G*(

**), a comparison with respect to**

*k**χ*(

**) should reflect the comparison with respect to**

*k**G*(

**). In fact, the IET scheme generates the most accurate**

*k**χ*(

**) and the qSTLS scheme, the least accurate**

*k**χ*(

**). The IET scheme strongly improves the HNC prediction for the extremum magnitude, but the HNC prediction for the extremum position is somewhat more accurate. (iii) Finally, the IET scheme’s static density response is again nearly exact for long wavelengths up to nearly**

*k**k*≲ 2

*k*

_{F}.

### D. Static approximation and effective static approximation

The *static approximation* utilizes the exact static LFC PIMC results to close the set of equations of the dielectric formalism.^{33,132} In other words, this approach replaces the exact dynamic LFC *G*(** k**,

*ω*) with its exact static limit

*G*(

**) =**

*k**G*(

**,**

*k**ω*= 0) that is accessible from PIMC simulations. The approach has been demonstrated to yield basically exact results for the dynamic structure factor and other related dynamic quantities in their entire non-trivial range, provided that

*r*

_{s}≲ 4, and to even reproduce the most prominent characteristics of the same quantities even at stronger coupling.

^{132}The approach has also been revealed to lead to accurate results for the SSF only for

*k*≲ 2.5

*k*

_{F}but to overestimate short-range correlations, a shortcoming that stems from the divergence of the exact static LFC short wavelength limit.

^{59}The SSF-related accuracy of the static approximation is expected to deteriorate further at strong coupling, where the dynamic nature of the LFC has been documented to be more important even in the classical limit.

^{52,78,79}In order to draw more quantitative conclusions, the exact static LFC PIMC results were employed as the closure of the exact two building blocks of the dielectric formalism; see Eqs. (3) and (4). It is evident that the static approximation has a low computational cost similar to that of the RPA.

Systematic comparisons revealed the following: (i) Regardless of the state point, the static approximation begins to overestimate the SSF prior to the main peak. For some state points, the overestimation takes place already at *k* ≲ 1.5*k*_{F}, i.e., right after the long wavelength range. (ii) As expected, the static approximation is much less accurate in the strongly coupled regime than in the warm dense matter regime. However, as *θ* increases, the accuracy level becomes significantly higher and the deviations become restricted at progressively shorter wavelengths, as revealed in Fig. 8. (iii) Since the static approximation either reproduces or overestimates the SSF in the entire wavenumber range, there is no possibility of favorable error cancellation during the computation of the interaction energy. As a result, the absolute relative deviations in the interaction energies always exceed 10% and can even reach 40%.

Moreover, given the success of the IIT interaction energy parameterization in the strongly coupled regime (see Sec. V A), it is tempting to check whether an extrapolation of the *effective static approximation* (ESA)^{59,60} scheme would also prove to be successful. The ESA utilizes an analytical representation for the frequency independent LFC,^{60} which is based on a neural net representation of PIMC static LFC results,^{131} but also incorporates the static LFC exact long wavelength limit as expressed by the CSR (employing the GDSMFB exchange–correlation free energy parametrization^{119}) and the frequency independent LFC exact short wavelength limit as expressed by the asymptotic self-consistency condition [employing a *g*(0) parameterization based on extrapolated restricted PIMC data^{59}]. Owing to the latter feature, the ESA scheme does not suffer from the aforementioned drawback of the static approximation. In fact, the ESA was recently demonstrated to be a fast and reliable tool for the computation of numerous thermodynamic, static, and dynamic quantities within the warm dense matter relevant range of 0.7 ≤ *r*_{s} ≤ 20 and 0 ≤ *θ* ≤ 4.^{59,60} Taking into account the satisfactory performance of the extrapolation of the GDSMFB parameterization toward low densities and the very small (albeit unphysical negative) values of the extrapolated *g*(0) parameterization toward low densities, the ESA performance for the paramagnetic electron liquid mainly depends on whether the analytical representation of the intermediate wavelength PIMC static LFC results can be successfully extrapolated toward high coupling parameters. In order to draw conclusions, the analytical ESA LFC was used as the closure of the exact two building blocks of the dielectric formalism; see Eqs. (3) and (4). It is evident that the ESA scheme also has a low computational cost similar to that of the RPA.

Systematic comparisons revealed the following: (i) As expected, the ESA LFC always extrapolates very well toward the long and the short wavelength ranges. (ii) At least for 0.75 ≤ *θ* ≤ 2.00 and regardless of *r*_{s}, the ESA well describes the position and the magnitude of the LFC peak. On the other hand, the ESA strongly overestimates the magnitude of the LFC peak for *θ* = 0.5 especially at strong coupling and underestimates the magnitude of the LFC peak for *θ* = 4.0. (iii) Regardless of the state point of interest, the main drawback of the ESA LFC lies at its transition to the short wavelength regime. This transition is always too abrupt and is not accompanied by small oscillations, as one would expect from a frequency independent LFC at strong coupling (see the IET LFC and the HNC LFC). These characteristics are not present in the warm dense matter regime or in dielectric schemes that are only equipped to describe weak correlations (see the qSTLS LFC). This drawback can be traced back to the implementation of a sigmoid activation function during the construction of ESA.^{60} (iv) Within 1 ≤ *θ* ≤ 2 and regardless of *r*_{s}, the ESA accurately estimates the position and magnitude of the SSF peak. On the other hand, the ESA strongly overestimates the magnitude of the SSF peak for *θ* = 0.5, 0.75 especially at strong coupling and underestimates the magnitude of the SSF peak for *θ* = 4.0. (v) Regardless of the state point, the ESA SSF directly transitions to its unity short wavelength limit after the Lorentzian peak and does not feature a shallow minimum followed by a second maximum. (vi) For most of the state points, the ESA SSF overestimates the PIMC SSF not only prior to the peak but also after. As a consequence, it leads to very inaccurate interaction energies. To sum up, the ESA does not extrapolate well to strong coupling, which can mostly be attributed to the abrupt transition from the LFC peak to the short wavelength range. Thus, it can be concluded that the future utilization of an ESA-like scheme at low densities would require a more involved activation function. Some characteristic examples are illustrated in Fig. 9.

## VI. SUMMARY AND DISCUSSION

In this work, we investigated the thermodynamic and the structural properties of the paramagnetic electron liquid with different dielectric theories and with *ab initio* path integral Monte Carlo simulations. First and foremost, we formulated a novel IET scheme that combines the dielectric formalism of many fermion systems with the integral equation theory of classical liquids. Essentially, the IET scheme incorporates a near-exact recent parameterization of the classical OCP bridge function into the existing HNC scheme. In addition, we carried out extensive PIMC simulations of the paramagnetic electron liquid for 16 state points. Combined with the recent availability of PIMC simulations for 20 state points, these new results substantially extend our current picture of the finite temperature UEF in the strongly coupled regime.

After an extensive comparison with “exact” PIMC results, various dielectric formalism schemes, and numerous extrapolated high-density parameterizations, we have demonstrated the following: (i) The IET scheme yields more accurate results for the interaction energy, static structure factor, static local field correction, and static density response of the paramagnetic electron liquid than all other known schemes (HNC, qSTLS, VS, and STLS). (ii) The IET interaction energies exhibit a remarkable agreement with the PIMC interaction energies having an accuracy within 0.68% and an average accuracy of 0.29%. This has been attributed to a favorable error cancellation in the course of the static structure factor integration. (iii) The IET long wavelength static local field correction, which is connected with the isothermal compressibility through the eponymous sum rule, is nearly indistinguishable from its “exact” PIMC counterpart up to nearly *k* ≲ 2*k*_{F}. (iv) The utilization of the IIT interaction energy parameterization leads to very accurate results for the paramagnetic electron liquid (similar accuracy level to the IET), in contrast to extrapolations of the GDSMFB, KSDT, and corrKSDT exchange–correlation free energy parameterizations. This has been attributed to the incorporation of exact ground-state and classical limit results and to the introduction of an accurate ad hoc *θ*− interpolation function. (v) The static approximation is much less accurate in the strongly coupled regime than in the warm dense matter regime. This confirms that the dynamic nature of the local field correction is an essential ingredient for the paramagnetic electron liquid. On the other hand, the analytic effective static approximation, which is very reliable in the warm dense matter regime, cannot be extrapolated to the strongly coupled regime primarily due to its direct transition to the short wavelength limit after the peak of the frequency independent local field correction.

At this point, it is important to bring forth the two main drawbacks of the IET scheme. *First*, the classical OCP bridge function is introduced as an analytic *b*(*r*, Γ) parameterization and not as a *b*[*h*] functional. As a consequence, the bridge function only reacts to the classical Coulomb interactions and not to quantum mechanical interactions (exchange degeneracy and diffraction effects). This shortcoming manifests itself in the lack of an appropriate ground-state limit, given the Γ = 2*λ*^{2}(*r*_{s}/*θ*) mapping of quantum to classical states that stems from the classical coupling parameter definition Γ = *e*^{2}/(*dk*_{B}*T*). A remedy can be found by recalling the ground-state quantum coupling parameter definition Γ_{q} = *e*^{2}/(*dE*_{F}) that leads to Γ_{q} = 2*λ*^{2}*r*_{s}. This suggests an effective coupling parameter definition $\Gamma eff=e2/{d[(kBT)2+EF2]1/2}$^{16} that leads to the new $\Gamma =2\lambda 2rs/1+\theta 2$ mapping. The utilization of the effective mapping *in lieu* of the classical mapping in the IET scheme led to nearly identical results for the state points investigated here. Nevertheless, this should be further tested for low-density state points near the ground-state (*θ* → 0). It should be mentioned that different empirical mappings could also be utilized based on enforcing consistency with the compressibility sum rule, based on enforcing consistency with some PIMC results (e.g., the static structure factor peak magnitude) or based on adopting empirical relations for a quantum temperature from the classical mapping method.^{133–135} Such options would probably lead to an improved IET scheme that is even more accurate, but they are not desirable because they add a strong phenomenological element to an otherwise rigorous theoretical approach. *Second*, quantum effects are treated on the random phase approximation level. A more sophisticated quantum treatment can be achieved by essentially combining the qSTLS and IET schemes. Given the studied static structure factors of these schemes, this seems to be a very promising strategy that could improve the peak magnitude prediction of the IET scheme. Unfortunately, it also adds an inconsistent element to the approach because the bridge function incorporation technique is based on the classical fluctuation dissipation theorem for a frequency independent local field correction, while the qSTLS scheme necessarily leads to a dynamic local field correction.

Future work will focus on the exploration of the aforementioned possibilities to remedy the drawbacks of the current version of the IET scheme. Moreover, the characterization of the IET scheme’s self-consistency concerning the compressibility sum rule and the application of the IET scheme for different values of the spin polarization parameter will also be pursued in the future.

Finally, benefitting from the recent availability of extensive path integral Monte Carlo simulation results for the finite temperature *θ* ∼ 1 paramagnetic electron fluid (non-ideal gas and liquid), the following practical guidelines can be formulated concerning the optimal dielectric formalism scheme. Within the warm dense matter regime of *r*_{s} ≲ 20, the semi-empirical ESA constitutes the most accurate scheme. At the moderate coupling 20 ≲ *r*_{s} ≲ 50 regime, the first-principle HNC scheme is the most accurate. At the strong coupling regime of *r*_{s} ≳ 50, the novel first-principle IET scheme is the most accurate. Naturally, the exact quantum coupling parameter boundaries between these three regimes depend on the specific degenerate parameter value.

## ACKNOWLEDGMENTS

The present work was partly funded by the Swedish National Space Agency under Grant No. 143/16. The present work was also partly funded by the Center for Advanced Systems Understanding (CASUS), which is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Culture and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament. All the PIMC simulations were carried out at the Norddeutscher Verbund für Hoch-und Höchstleistungsrechnen (HLRN) under Grant No. shp00026 and on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at the Technische Universität Dresden. The IET, HNC, and qSTLS schemes were numerically solved on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the NSC (Linköping University) that is partially funded by the Swedish Research Council under Grant Agreement No. 2018-05973.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request and within the article.