Integral equation theory based dielectric scheme for strongly coupled electron liquids

In a recent paper, Lucco Castello et al. [arXiv:2107.03537] provided an accurate parametrization 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. 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 (referred to also 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][4][5][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][3][4][5][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, neutron star crusts) and ultra-fast laser heating or shock compression of metals [17,18], has provided the impetus for 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][21][22][23].
Owing to the fact that low density finite temperature inhomogeneous electron systems are inaccessible even to 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, in spite of 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][25][26], the emergence of a spin-density wave [26][27][28], the possibility of a continuous paramagnetic to ferromagnetic transition [29][30][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][37][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 parametrization 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.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][41][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 analogue 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 ) the first Bohr radius, the degeneracy parameter θ = k B T /E F where E F = [(6π 2 n ↑ ) 2/3 /2]( 2 /m e ) is the Fermi energy w.r.t the Fermi wave-vector of spin-up electrons k ↑ F = (6π 2 n ↑ ) 1/3 , the spin polarization parameter ξ = (n ↑ − n ↓ )/n that is constrained within 0 ≤ ξ ≤ 1 under the spin-up choice n ↑ ≥ n ↓ .In the above, is the reduced Planck constant, k B the Boltzmann constant, e the electron charge, m e the electron mass, n = n ↑ + n ↓ the electron density, T 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][44][45] for which Γ = 2λ 2 r s /θ with λ 3 = (k F d) −3 = 4/(9π).
In what follows, we shall restrict ourselves to the spinunpolarized (paramagnetic) case, n ↑ = n ↓ or ξ = 0.The ideal (Lindhard) density response of free electrons, i.e. in absence of Coulomb interactions, has the form [19,43,46] χ 0 (k, ω) = 2 d 3 q (2π) 3  f 0 (q) − f 0 (q + k) ω + ǫ(q) − ǫ(q + k) + ıη , with ǫ(q) = 2 q 2 /(2m e ) the electron kinetic energy, η → 0 + owing to the adiabatic switching of the perturbation at t → −∞ that ensures causality and f 0 (q) the Fermi-Dirac distribution function that is given by with μ = βµ the reduced chemical potential (µ is the chemical potential and 1/β = k B T ) that is determined by the normalization condition [d 3 q/(2π) 3 ]f 0 (q) = 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 & quantum effects.

B. The dielectric formalism
The dielectric formalism is based on (i) the quantum fluctuation dissipation theorem which connects the dynamic structure factor S(k, ω) (DSF) 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 G(k, ω) (LFC), (iii) a functional relation connecting G(k, ω) with the static structure factor S(k) (SSF) that is obtained by perturbative quantum / classical BBGKY approaches or non-perturbative integral equation approaches [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 [19,43].The extension of the χ(k, ω) domain with the aid of analytic continuation leads to the complex valued χ(k, z) and allows 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.
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][51][52] with U (k) = 4πe 2 /k 2 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 selfconsistent 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 in 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 that 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][37][38]  The bridge function is formally defined by a virial-type expansion whose unknown coefficients are given by multidimensional integrals with the integrands being products of Mayer functions f (r) = exp [−βu(r)] − 1.The operation of topological reduction allows 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 & thermodynamic results for Γ 1 [64], but its accuracy strongly degrades as crystallization is approached [65][66][67].
A computationally costly but far more accurate alternative concerns the indirect bridge function extraction from computer simulations [68][69][70].Extraction efforts for the classical OCP date back to the 80's [71][72][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 as well as construct an accurate analytic representation valid along the short range and intermediate range [39,75].In reduced units x = r/d, this parametrization 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, (ln Γ) j are monotonic functions of Γ and the s j i , l j i coefficients have been tabulated [39,75].The closed system of Eqs.(6,7,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 χ(k, 0) = π −1 [ℑ{χ(k, ω)}/ω]dω 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(k), employing the Maxwellian result for the static limit of the ideal classical density-density response χ 0 (k, 0) = −nβ, using the connecting relation S(k) = 1 + nH(k) and invoking the Fourier transformed OZ equation 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,7), yields Equating these r.h.s., 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(k), which leads to Using again the connecting relation S(k) = 1 + nH(k) in order to dispose of the Fourier transformed total correlation function and substituting for U (k) = 4πe 2 /k 2 , ultimately results 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(|q|), 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 [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][85][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 = p + (1/2)k, spherical coordinates for p assuming k||ẑ without loss of generality, and u = |p + (1/2)k| , w = |p − (1/2)k| with a surface element identity 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 , 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(|q|), although this would lead to a drastic reduction of 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 [35,87], where the triple integral is converted to a single integral after replacing S(k − q) by an ad hoc Yukawa screening function of a characteristic screening wave-number that is determined by an interaction energy constraint.

B. On the convergence of the infinite Matsubara series
Independent of the dielectric scheme (STLS, CA, qSTLS, HNC, IET), the Matsubara summation of Eq.( 5) is slowly converging, especially for small values of the degeneracy parameter θ.The implementation of mathematical tricks that speed-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 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 Φ(k, z) = −(2E F )/(3n) χ 0 (k, z) and normalized wave-vectors x → k/k F , Eq.(5) becomes , where, courtesy of logarithm properties and the product given by the 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 Φ ∞ (x, l) = (4/3)x 2 /[x 4 + (2πlθ) 2 ] + O(x −4 , l −4 ), 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 where coth (•) is the hyperbolic cotangent and csch(•) 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 = 51200 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(k) implies the (3π/8)x 4 [1−S(x)] = λr s [1−G(x)] asymptotic limit [49].When combined with Kimball's expression [89] ∂g( r = 0)/∂ r = (3π/8) lim x→∞ x 4 [1−S(x)] with r = rk F and the cusp relation [89,90] ∂g( r = 0)/∂ r = λr s g(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, HNC) as well as 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] 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 cut-off 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 r.h.s. 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 of variables (y, z) → (t, s) of the form y = √ x 2 + s 2 − 2xst, z = s and carrying out the t− integration.This leads to In the above, 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 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,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], the IET LFC double integral expression after utilizing bipolar coordinates and splitting the STLS LFC single integral expression [see Eq.( 16), Secs.III A,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)/βU (q) = (q/Γ) ∞ 0 xb(x, Γ) 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 For the computation, Eq.( 14) was combined with Eq.( 8).The numerical integration was performed with the doubly adaptive Clenshaw-Curtis quadrature method.
states (Γ) via Γ = 2λ 2 (r s /θ).The B(q)/βU (q) contribution has been 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.1k F grid resolution and with a 40k F upper cut-off.The only exception was the complete Fermi-Dirac integral I 1/2 (μ), 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.005k 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.
(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 r.h.s. 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' 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][95][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 with increasing the system size N or decreasing the temperature T ; the reader is addressed to Refs.[99,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 = (r 1 , . . ., r N ) 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 S N i , i ∈ {↑, ↓}, where the sign function sgn(σ ↑ , σ ↓ ) gives positive (negative) unity for an even (odd) number of pair permutations.Further, the operators πσ 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 −β Ĥ cannot be readily evaluated, since the kinetic ( K) and potential ( V ) contributions to the Hamiltonian Ĥ = K + V do not commute, e −β Ĥ = e −β K e −β V .In order to overcome this obstacle, we utilize the exact semi-group property of ρ where the definition ǫ = β/P has been employed.Applying Eq.( 18) to Eq.( 17) and inserting P −1 unity operators of the form 1 = dR α |R α R α | leads 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 hightemperature approximation, like the primitive factorization e −ǫ Ĥ ≈ e −ǫ K e −ǫ V .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] lim 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 that are bounded from below, which is indeed the case for the UEF.Attractive potentials like 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 abbreviated form as where the integration over the X = (R 0 , . . ., R P −1 ) 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 = 3P Ndimensional 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 wellknown 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 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 ω p = 3/r 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) below, 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 has been 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 grey 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 finitesize 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 has been 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 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 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 published results of Schweng & Böhm [62].

A. Interaction energy
The interaction energy per particle is generally obtained by . Substituting for the Coulomb potential energy, introducing the (r s , θ) variables and employing normalized wavenumbers, this leads to the standard UEF expression [19,43] where u(r s , θ) 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 & IET schemes have been listed in Table I.The absolute relative deviations between the theoretical and "exact" interaction energies have also been 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][66][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) In spite of their very high accuracy, the IET (as well as the HNC) interaction energies cannot reproduce the exact θ-dependence of u(r s , θ) within the highly degenerate θ range.To be more specific, the PIMC results reveal a monotonic | u| decrease as θ increases, while the IET (and HNC) results exhibit a monotonic | u| increase as θ increases within θ 1 and the correct monotonic | u| decrease only as θ increases within θ 1.For this reason, we did not construct an exchange-correlation free energy parametrization via the adiabatic connection formula.
It is important to point out that the employed classical OCP bridge function parametrization 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 towards 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 parametrizations of the exchange-correlation free energy through the expression that is acquired after the differentiation of the thermodynamic formula f xc (r s , Θ) = r −2 s rs 0 r ′ s u int (r ′ s , Θ)dr ′ s [19].In particular, we consider the parametrization by Groth et al. (GDSMFB) that is based on simulation results obtained by various novel PIMC methods within 0.1 ≤ r s ≤ 20 & 0.5 ≤ θ ≤ 8 [119], the parametrization by Karasiev et al. (KSDT) that is based on simulation results obtained by the restricted PIMC method within 1.0 ≤ r s ≤ 40 & 0.0625 ≤ θ ≤ 8 [120] and the corrected version of the parametrization 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 in spite of 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 TABLE I: The interaction energy u(rs, θ) per particle (expressed in Hartree units) of the paramagnetic electron liquid: comparison of the finite-size corrected "exact" PIMC results with the predictions of the qSTLS, HNC and IET dielectric formalism schemes.The PIMC data for the first 6 state points are adopted from Ref. [34], while the PIMC data for the remaining 14 state points are new.The absolute relative deviations between the dielectric scheme results and the PIMC results are also reported.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 parametrization is much more accurate than the KSDT parametrization 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, 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 parametrization 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 parametrization might be very accurate.This is verified by the comparison which 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 9 of the studied states, while the mean absolute relative errors are 0.29% for the IET, 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 have been illustrated in Fig. 4 and Fig. 5.For all the 20 investigated state points, the magnitudes and positions of the first SSF peak resulting from the PIMC simulations as well as from the qSTLS, HNC & IET schemes have been listed in Table II.The absolute relative deviations between the theoretical and the "exact" SSF peak values have also been tabulated therein.
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 3k F (especially for r s 100) and the wellresolved second SSF peak above 4k 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 2k F , in the Lorentzian shaped region that surrounds the first maximum k ∼ 2 − 2.5k F and within the short wavelength range of k 2.5k 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 towards 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 concluded that the first SSF peak position is mainly controlled by strong correlations.(iv) Regardless of 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,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 ≤ 2k F , becomes clearly too small within 2.0k F ≤ k ≤ 2.5k F and tends to be somewhat too large within 2.5k F ≤ k ≤ 4.0k 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 g(r) (RDFs) extracted from PIMC simulations and computed with the three dielectric formalism schemes.Since 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 & 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 & 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.

S(q)
FIG. 4: Dependence of the paramagnetic electron liquid static structure factor on the quantum coupling parameter rs.Results from the IET scheme (red solid line), the HNC scheme (blue dashed line), the qSTLS scheme (green dot-dashed line) and PIMC simulations (black crosses) for θ = 0.50 and varying rs = 50, 100, 150, 200 in the normalized wavenumber range k ≤ 4.5kF.The superiority of the new IET scheme within the long wavelength range, in the maximum vicinity and within the short wavelength range is obvious.

S(q)
FIG. 5: Dependence of the paramagnetic electron liquid static structure factor on the quantum degeneracy parameter θ. Results from the IET scheme (red solid line), the HNC scheme (blue dashed line), the qSTLS scheme (green dot-dashed line) and PIMC simulations (black crosses) for rs = 125 and varying θ = 0.5, 1.0, 1.5, 2.0 in the normalized wavenumber range k ≤ 4.5kF.The superiority of the new IET scheme within the long wavelength range, in the maximum vicinity and within the short wavelength range is obvious.
FIG. 6: Dependence of the paramagnetic electron liquid static local field correction (main figure) and static density response (inset figure) on the quantum coupling parameter rs.Results from the IET scheme (red solid line), the HNC scheme (blue dashed line), the qSTLS scheme (green dot-dashed line) and our PIMC simulations (black crosses) for θ = 1.00 and varying rs = 50, 100, 150, 200 in the normalized wavenumber range k ≤ 4.5kF.For both quantities, the superiority of the new IET scheme within the long wavelength range and in the vicinity of the maximum is obvious.
FIG. 7: Dependence of the paramagnetic electron liquid static local field correction (main figure) and static density response (inset figure) on the quantum degeneracy parameter θ. Results from the IET scheme (red solid line), the HNC scheme (blue dashed line), the qSTLS scheme (green dot-dashed line) and our PIMC simulations (black crosses) for rs = 100 and varying θ = 0.5, 1, 2, 4 in the normalized wavenumber range k ≤ 4.5kF.For both quantities, the superiority of the new IET scheme within the long wavelength range and in the vicinity of the maximum is obvious.

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 have been illustrated in Fig. 6 and Fig. 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.60k F , which is followed by a shallow local minimum located at k ∼ 3.05 − 3.45k 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 Fig. 7(a-d).(ii) Regardless of the state point of interest, the IET scheme generates the most accurate and the qSTLS scheme the least accurate static LFC.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 2k 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 exchangecorrelation free energy and not the exact exchangecorrelation 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 parametrization 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" χ(k) of the paramagnetic electron liquid has been discussed earlier [34] and is confirmed by our PIMC results.It is evident that the magnitude of the pronounced χ(k) extremum increases as r s decreases (see Fig. 6) and as θ decreases (see Fig. 7) and that the width of the χ(k) extremum becomes more sharp as r s increases (see Fig. 6) and especially as θ decreases (see Fig. 7).(ii) Since χ(k) is obtained by setting ω = 0 to Eq.( 4) and only involves the common among schemes χ 0 (k) and the varying among schemes G(k), comparison with respect to χ(k) should reflect the comparison with respect to G(k).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 2k 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.5k F but to overestimate short-range correlations; a shortcoming which 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,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.5k 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 parametrization in the strongly coupled regime (see section 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 exchangecorrelation 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) parametrization 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 parametrization towards low densities as well as the very small (albeit unphysical negative) values of the extrapolated g(0) parametrization towards 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 towards 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,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 towards the long and the short wavelength ranges.(ii) At least for 0.75 ≤ θ ≤ 2.00 and regardless of r s , the ESA 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 coupling and underestimates the magnitude of the LFC peak for θ = 4.0.(iii) Regardless of the state point FIG. 9: The paramagnetic electron liquid's static local field correction (main) and static structure factor (inset) for rs = 125 and θ = 0.75 (a), θ = 1.00 (b), θ = 1.50 (c).Results from the IET scheme (red solid line), the ESA scheme (cyan dashed line) and PIMC simulations (black crosses) in the normalized wavenumber range k ≤ 4.5kF.The ESA shortcoming that concerns the abrupt LFC (and thus SSF) transition from the peak vicinity to the short wavelength limit is apparent.
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 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
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 parametrization of the classical OCP bridge function to 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 parametrizations, we have demonstrated that; (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, 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, that is connected with the isothermal compressibility through the eponymous sum rule, is nearly indistinguishable from its "exact" PIMC counterpart up to nearly k 2k F .(iv) The utilization of the IIT interaction energy parametrization 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 parametrizations.This has been attributed to the incorporation of exact ground-state and classical limit results as well as to the introduction of an accurate ad hoc θ−interpolation function.(v) The static approximation is much less accurate in the strongly coupled regime than 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, that 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, Γ) parametrization 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 Γ eff = e 2 /{d[(k B T ) 2 + E 2 F ] 1/2 } [16] that leads to the new Γ = 2λ 2 r s / √ 1 + θ 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][134][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 lo-cal 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.
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(k) 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 G(k).

FIG. 2 :
FIG.2: Illustration of the PIMC method.A configuration of N = 4 particles is shown in the τ -x-plane for P = 6 hightemperature factors.Notice the single pair exchange of the two particles in the center leading to a negative sign for this particular configuration, W (X) < 0.
FIG.3:The finite-size extrapolation of the PIMC data.The UEF interaction energy per particle is shown at rs = 125 and θ = 2 as a function of the inverse electron number.The red circles depict the raw PIMC results for different numbers of electrons N and the green crosses have been obtained by adding onto the former the finite-size correction from Eq.(23).The horizontal blotted blue line shows the average value of the corrected data points, and the light dotted grey lines indicate a deviation of 0.01%.

TABLE II :
[34]peak magnitude S max and peak position arg q S max of the static structure factor S(k/kF; rs, θ) of the paramagnetic electron liquid: comparison of the "exact" PIMC results with the predictions of the qSTLS, HNC and IET dielectric formalism schemes.The PIMC data for the first 6 state points are adopted from Ref.[34], while the PIMC data for the remaining 14 state points are new.The absolute relative deviations between the dielectric scheme results and the PIMC results are also reported.