Revisiting the Vashishta-Singwi dielectric scheme for the warm dense uniform electron fluid

The finite temperature version of the Vashishta-Singwi (VS) dielectric scheme for the paramagnetic warm dense uniform electron fluid is revisited correcting for an earlier thermodynamic derivative error. The VS scheme handles quantum mechanical effects at the level of the random phase approximation and treats correlations via the density expansion of a generalized Singwi-Tosi-Land-Sj\"olander (STLS) closure that inserts a parameter determined by enforcing the compressibility sum rule. Systematic comparison with quasi-exact results, based on quantum Monte Carlo simulations, reveals a structural superiority of the VS scheme towards strong coupling and a thermodynamic superiority of the STLS scheme courtesy of a favorable cancellation of errors. Guidelines are provided for the construction of dielectric schemes that are expected to be more accurate but computationally costly.

The main pre-requisite for the understanding of WDM is the accurate description of the properties of the warm dense uniform electron fluid (UEF), a fundamental homogeneous model system that does not require explicit treatment of the ionic component [30][31][32].For instance, the parametrization of the UEF exchange correlation free energy constitutes essential input for thermal density functional theory (DFT) [33][34][35][36][37], while the dynamic local field correction of the UEF is formally equivalent to the exchange correlation kernel of linear-response timedependent DFT [38][39][40] and can be utilized to include correlations into quantum hydrodynamic models [41,42].
With this investigation, we aim to improve the theoretical landscape by revisiting the finite temperature version of the Vashishta-Singwi dielectric scheme [70] which supplements the archetypal Singwi-Tosi-Land-Sjölander closure [69] with a density expansion that contains an adjustable parameter determined by the exact enforcement of the compressibility sum rule.Our self-consistent formulation corrects for a thermodynamic derivative error that was present in the earlier formulation of Sjostrom and Dufty [71].An efficient algorithm is devised that is implemented in a hybrid C++/python code.The scheme is numerically solved for thousands of paramagnetic UEF state points that span the WDM regime.The predictions for thermodynamic (interaction energy, exchange correlation free energy) and structural properties (static structure factor, static density response function, pair correlation function, static local field correction) are systematically compared with available quasi-exact results based on path integral Monte Carlo (PIMC) simulations.Finally, we hint at future possibilities for constructing more accurate, albeit more complex, dielectric schemes.

II. THEORETICAL
The UEF is a spatially homogeneous quantum model system consisting of electrons immersed in a rigid uniform ionic background that is solely characterized by the neutralizing charge density −ne, where n is the electron density [30][31][32].It constitutes the quantum mechanical analogue of the classical one-component plasma (OCP) [90][91][92], thus it is often also referred to as the quantum OCP.The UEF thermodynamic state points are specified by three dimensionless parameters [30][31][32]: (1) the quantum coupling parameter r s = d/a B with d the Wigner-Seitz radius d = (4πn/3) −1/3 and a B = 2 /(m e e 2 ) the Bohr radius, (2) the quantum degeneracy parameter Θ = T /E F with E F = [(6π 2 n ↑ ) 2/3 /2]( 2 /m e ) the Fermi energy with respect to the Fermi wave-vector of the spin-up electrons k ↑ F = (6π 2 n ↑ ) 1/3 and T the temperature in energy units, (3) the spin polarization parameter ξ = (n ↑ − n ↓ )/n for which 0 ≤ ξ ≤ 1 within the standard convention n ↑ ≥ n ↓ .In this work, we exclusively study the paramagnetic (or unpolarized) case of equal spin-up and -down electrons, ξ = 0.It is noted that, in contrast to the UEF, OCP state points are specified by a unique dimensionless parameter [90][91][92]: the classical coupling parameter Γ = e 2 /(dT ) for which Γ = 2λ 2 r s /Θ with λ 3 = (k F d) −3 = 4/(9π).
Within the UEF phase diagram, the WDM regime can be roughly demarcated by 0.1 r s , Θ 10 [31,32].The lack of small parameters, r s ∼ Θ ∼ 1, is indicative of the complex interplay between quantum effects (exchange, diffraction), bare Coulomb interactions and thermal excitations which make the theoretical treatment of the warm dense UEF a formidable task.

A. Dielectric formalism
The self-consistent dielectric formalism constitutes one of the most accurate and versatile microscopic frameworks for the description of the thermodynamic and static properties of interacting uniform quantum systems such as the UEF [30][31][32][86][87][88][89].It combines fundamental results of the density version of linear response theory [30] with an approximate closure stemming from perturbative quantum/classical kinetic theories of non-ideal gases [95,96] or from non-perturbative integral equation theories / memory function approaches of classical liquids [96,97].
In the polarization potential approach [96], the density response function χ(k, ω) is expressed very generally in terms of the ideal (Lindhard) density response χ 0 (k, ω) and the dynamic local field correction G(k, ω) (LFC) as where U (k) is the Fourier transformed pair interaction energy, with U (k) = 4πe 2 /k 2 for Coulomb interactions.
In addition, for finite temperature systems, the combination of the zero frequency moment sum rule, the quantum fluctuation-dissipation theorem and the analytic continuation of the causal χ(k, ω) to the complex frequency plane yield a static structure factor S(k) (SSF) relation that involves the infinite Matsubara summation [68] S with χ(k, z) the analytically-continued density response function and ω l = 2πl/(β ) the bosonic Matsubara frequencies.Finally, the dynamic LFC, which incorporates Pauli exchange, quantum diffraction and Coulomb correlation effects beyond the mean field description, is given by a complicated SSF functional of the general form The combination of Eqs.(4,5,6) leads to a non-linear functional equation of the type to be solved for the SSF [31,68].Multiple dielectric schemes have been developed over the last six decades that differ in the approximate treatment of the exact (unknown) LFC functional.They can be broadly categorized in the following manner; (i) semiclassical schemes based on classical kinetic equations [68][69][70][71], (ii) pure quantum schemes based on quantum kinetic equations [72][73][74][75], (iii) semi-classical schemes based on various integral equation theories [76][77][78][79][80], (iv) quantized schemes also based on integral equation theories [81], (v) memory-function or viscoelastic schemes incorporating the elusive third frequency moment sum rule [82][83][84][85], (vi) semi-empirical schemes based on exact QMC results and incorporating known asymptotic limits [62,63].
B. Vashishta-Singwi scheme: General formulation Vashishta-Singwi (VS) theory is a semi-classical scheme of the self-consistent dielectric formalism [70].The treatment of correlations is classical, since the s = 1 member of the classical Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy of s-reduced distribution functions is truncated with an equilibrium factorization Ansatz.The treatment of quantum mechanical effects is restricted to the random phase approximation level, since the resulting ideal Vlasov density response is substituted with the ideal Lindhard density response.It is noted that equilibrium closures to the classical BBGKY hierarchy imply a static LFC (SLFC), thus the VS closure leads to a SLFC.
After substitution of the VS closure to the first member of the classical BBGKY hierarchy and use of g eq (|r − r ′ |) = 1 + h(|r − r ′ |) with h(•) the total correlation function, the VS kinetic equation reads as with U ext (r, t) the external potential energy and U (r−r ′ ) the pair interaction potential energy.At the RHS, the first bracketed term represents the mean field contribution, the second bracketed term represents the STLS polarization field contribution and the last two bracketed terms correspond to the VS polarization field corrections.
Application of an external potential energy perturbation to the system (initially in equilibrium), linearization with respect to the small perturbation strength, use of spatiotemporal Fourier transforms, consideration of the adiabatic switching of the perturbation and substitution for the SSF via Integration over the momenta to obtain the density perturbation δn k,ω = δf p k,ω d 3 p, solution of the resulting linear equation with respect to the density perturbation δn k,ω , identification of the ideal Vlasov density response .
After term-by-term comparison with the general form of the density response function, Eq.( 4), and substitution for the Fourier transformed Coulomb interaction U (k) = 4πe 2 /k 2 , it is evident that the VS SLFC is given by The second bracketed factor can be identified to be equal to the STLS SLFC.Thus, the VS SLFC can be compactly rewritten as Finally, the CSR is most generally expressed through the long-wavelength limit of the static density response (SDR) as [30,87,89] where f (n, T ) denotes the total free energy per particle.Application of the above formula for the interacting and non-interacting case, subtraction by parts, substitution for the general density response function from Eq.( 4) and introduction of the exchange-correlation free energy (per particle) f xc , allows to alternatively express the CSR via the long-wavelength limit of the SLFC as [31,32] Taking into account that G VS (k) explicitly depends on α and that f xc (n, T ) implicity depends on α, Eq.( 9) essentially constitutes a non-linear equation that determines the self-consistency parameter α for every state point.It is also important to emphasize that the 1/2 ≤ α 1 UEF ground state bounds [70] and the α ≃ 2/3 UEF ground state approximation [70] are not applicable to the finite temperature UEF.
C. Vashishta-Singwi scheme: Correct finite temperature formulation In the standard normalized units, for the correct calculation of the density derivative present in the VS SLFC, one should take into account that the STLS SLFC depends on the density through the quantum coupling parameter r s (see the Wigner-Seitz radius), the quantum degeneracy parameter Θ (see the Fermi energy) and the normalized wavenumber x = k/k F (see the Fermi wavenumber).Thus, the chain differentiation rule leads to the VS SLFC with the STLS SLFC given by [31,68] In the standard normalized units where the energies are normalized with the Hartree energy, the CSR becomes The RHS is first computed by taking into account that the exchange-correlation free energy depends on the density through the quantum coupling parameter r s (see the Wigner-Seitz radius) and the quantum degeneracy parameter Θ (see the Fermi energy).Thus, repeated applications of the chain differentiation rule yield The LHS is then computed with the aid of the straightforward [see Eq.( 11)] long-wavelength limits where the fundamental expression for the interaction energy of the UEF [30,31,89] has been employed Substituting for the VS SLFC from Eq.( 10) and utilizing the above two limits, one expresses the LHS as Finally, combining the RHS, LHS expressions with the CSR expression, solving the linear explicit dependence with respect to the self-consistency parameter α(r s , Θ), utilizing the differential version of the adiabatic connection formula in order to symmetrize the numerators and denominators [30,31,89] and carrying out some simple algebraic manipulations, the CSR expression for the finite temperature VS scheme is ultimately transformed to a compact equation for the self-consistency parameter α.It is emphasized that this equation is highly non-linear owing to the implicit depen-dence of the interaction energy and exchange-correlation energy on α, as discerned by inspecting Eqs.(12,13) in view of the α−dependent VS SLFC, see Eqs. (10,11).The compact equation reads as To sum up, in normalized units, the correct finite temperature formulation of the VS closure consists of Eqs.(10,11) for the SLFC and Eq.( 15) for the CSR.As aforementioned, the VS scheme was originally formulated for the UEF ground state [70].It is noted that the ground state assumption has strong implications for the self-consistent dielectric formalism itself; when T = 0 there are no hyperbolic cotangent poles in the frequencyintegrated fluctuation-dissipation theorem, the bosonic Matsubara frequencies collapse and the Matsubara infinite sum is substituted with a positive frequency integral.The first rigorous attempt to extend the VS scheme to the finite temperature UEF was due to Sjostrom and Dufty [71].In what follows, we shall refer to it as VS-SD.Unfortunately, when progressing from the general formulation to the finite temperature formulation, the authors did not take into account the density dependence of the degeneracy parameter Θ in both the VS SLFC and the CSR expression, see their Eq.(15).As a consequence, the VS-SD SLFC reads as and the VS-SD self-consistency parameter equation that stems from the CSR expression reads as It is straightforward that the VS-SD formulation can be directly obtained from the correct VS formulation by setting Θ = 0 whenever Θ appears explicitly, leaving the implicit Θ−dependence of thermodynamic and structural properties intact.It is also apparent that the VS-SD finite temperature formulation results in the same SLFC closure and the same self-consistency parameter equation as the original VS ground state formulation.Nevertheless, the dielectric formalism machinery and thus the numerical solution of the two formulations are different.In spite of the thermodynamic error of the Sjostrom and Dufty analysis, the VS-SD finite temperature scheme is still valuable, since it is numerically much simpler as well as computationally much less costly than the correct VS finite temperature scheme but also since the two schemes should lead to identical predictions at the limit of high degeneracy and / or strong coupling.To sum up, in normalized units, the SD finite temperature formulation of the VS closure consists of Eqs.(11,16) for the SLFC and Eq.( 17) for the CSR.

III. COMPUTATIONAL DETAILS A. Set of equations
The closed normalized set of equations emerges by combining the general ingredients of the dielectric formalism with the specific VS closure.The VS equations feature: (1) The normalization condition of the Fermi-Dirac energy distribution function that allows the computation of the reduced chemical potential μ = βµ [see Eq.( 18)].
(3) The VS SLFC as extensively discussed in Sec.II C [see Eq.( 21)].(4) The Matsubara summation expression for the SSF [see Eq.( 22)].In practice, in the numerical implementation, the auxiliary complex function Φ(x, l) and the square of the combined high frequency large wavenum-ber asymptotic expansion of the square of the auxiliary complex function Φ 2 (x, l → ∞) are split up from the infinite series coefficients.Then, the resulting infinite sums can be evaluated exactly leading to the non-interacting (Hartree-Fock) SSF S HF (x) and a residual SSF correction S ∞ (x) [68,74,99].As a direct consequence, there is a sig-nificant acceleration of the convergence rate of the Matsubara summation even at strong degeneracy [79][80][81]. (5) The non-linear equation for the theoretical value of the self-consistency parameter a th , i.e. the value that exactly satisfies the CSR, as extensively discussed in Sec.II C [see Eq.( 23)].

B. Numerical algorithm
The closed normalized set of equations is solved in a home-made hybrid code, where C++ is used for the backend (where all the numerics are handled) and python is used for the frontend (where all the user inputs are introduced, all scheme outputs are extracted and the primary data post-processing is done).The first step concerns the computation of the reduced chemical potential μ(Θ) for each value of the degeneracy parameter.The solution of Eq.( 18) is found by applying the Brent-Dekker hybrid root-finding algorithm.The second step concerns the computation of the ideal Lindhard density response for all Matsubara frequencies from Eqs. (19,20).For state points close to the ground state, Θ 1, l = 500 Matsubara frequencies are required for truncated series convergence.For larger values of the degeneracy parameter, the number of Matsubara frequencies can be safely lowered down to l = 300.The improper integrals are handled with the doubly-adaptive generalpurpose quadrature routine CQUAD of the GSL library; a 0.1 grid resolution and 50 upper cut-off are employed.
The third step concerns the computation of the VS SLFC from Eq.( 21).The STLS SSF is used as the initial value for the SSF.Two initial values are then required for the self-consistency parameter, α num .As a practical rule, the lower the quantum coupling parameter is, the lower the chosen α num initial values should be.The reason behind the necessity for computations for two initial values will become apparent when α th is computed in the fifth step upon the enforcement of the CSR.Moreover, the derivatives with respect to the state point variables (r s , Θ) need to be computed.Therefore, as deduced from the computational stencil that is shown in Fig. 1, the VS scheme needs to be solved simultaneously for the nine state points (r s , Θ), (r s + ∆r s , Θ), (r s − ∆r s , Θ), (r s , Θ + ∆Θ), (r s , Θ − ∆Θ), (r s +∆r s , Θ+∆Θ), (r s +∆r s , Θ−∆Θ), (r s −∆r s , Θ+ ∆Θ), (r s − ∆r s , Θ − ∆Θ).The combination of central, forward and backward difference approximations is such that there is no need to supply boundary conditions outside the stencil in order to close the system of equations.In addition, the x-derivatives are just computed on a grid depending on how fine the normalized wavevector discretization is.To be more specific, we typically have ∆Θ = 0.1, ∆r s = 0.1 and ∆x = 0.1.With this method, the different state points can be easily stored in a twodimensional matrix G[i][j] with the x values stored as the i-components and the (r s , Θ) grid as the j−components.Since all computations are done up to x max = 50, i cor- FIG. 1: 2D grid for the finite difference approximation of the state point derivatives.
responds to 500 and j to 9. For the r s state points, which are (r s , Θ), (r s , Θ + ∆Θ), (r s , Θ − ∆Θ) and with j = {3, 4, 5}, the central second order difference is used: For the r s +∆r s state points, which are (r s +∆r s , Θ), (r s + ∆r s , Θ + ∆Θ), (r s + ∆r s , Θ − ∆Θ) and with j = {0, 1, 2}, the forward second order difference is used: For the r s −∆r s state points, which are (r s −∆r s , Θ), (r s − ∆r s , Θ + ∆Θ), (r s − ∆r s , Θ − ∆Θ) and with j = {6, 7, 8}, the backward second order difference is used: The same is done for the Θ derivatives, but now we have j = {1, 4, 7} for the central difference, j = {0, 3, 6} for the forward difference and j = {2, 5, 8} for the backward difference.Finally, the derivatives with respect to x are computed with the central second order difference: It is emphasized that all finite difference approximations are second order which implies that the truncation errors . The combination of central, backward and forward differences was chosen so that the grid extent is minimized.The fourth step concerns the computation of the VS SSF from Eq.( 22) for the nine relevant UEF state points.Then, steps 3 and 4 are repeated until the convergence criterion i {G for the central point (r s , Θ), where (n), (n− 1) denote the current and previous iteration values of G VS .From now on, this will be referred to as the "inner loop".Typically, the number of iterations required to reach convergence for G VS per numerical self-consistency parameter per state point is of the order of ∼ 100.
The fifth step concerns the convergence of the selfconsistency parameter α num so that the CSR is automatically enforced.The interaction energy u int is computed from the SSF, see Eq. (12).The exchange correlation free energy f xc is computed from the adiabatic connection formula, see Eq.( 13).The respective wavenumber and coupling parameter integrations are again performed with the doubly-adaptive general-purpose quadrature routine CQUAD.Various thermodynamic derivatives of u int , f xc need to be computed.All the involved derivatives are approximated with central differences.For the first-order derivatives (below for u int but also for f xc ), we have For the second-order derivatives, we have These allow the computation of α th from Eq.( 23).After convergence of the initial inner loop, the G

VS
that correspond to the guesses α num yield the values α (k) th (r s , Θ), α (k+1) th (r s , Θ).The root of f (α th , α num ) = α th (r s , Θ)−α num = 0, which yields the α num that satisfies the CSR near exactly, is found with the secant method.After the inner loop has converged, At this point, the fifth step, which constitutes the "outer loop" is concluded, and the inner loop is restarted with initial guesses α num .For the full VS scheme to be solved, the inner and outer loops must be repeated until both convergence criteria are satisfied.On average, the outer loop requires less than five iterations before convergence is reached.Hence, the home-made VS solver is very fast.Typically, the solution of the VS scheme per UEF state point requires 5 sec on a conventional laptop.

IV. RESULTS
In this section, the numerical solutions of our finite temperature VS scheme for thermodynamic quantities and static structural properties (SSF, PCF, SDR, SLFC) will be compared with the predictions of the STLS scheme and the VS-SD scheme.
The VS thermodynamic properties will also be compared with highly accurate warm dense UEF equations of state [31,60,61,100,101].The GDSMFB parametrization of the UEF exchange correlation free energy f xc will be preferred [31,60].The GDSMFB parametrization is based on the Ichimaru-Iyetomi-Tanaka (IIT) functional form for u int [60,88,102].The classical analytical Debye-Hückel limit (Θ → ∞) [90] as well as the Perrot and Dharma-wardana closed-form approximation of the Hartree-Fock limit (r s → 0) [103] are exactly incorporated.The GDSMFB parametrization primarily utilizes finite-size corrected QMC results [54] obtained by various novel PIMC methods within 0.1 ≤ r s ≤ 20 and 0.5 ≤ Θ ≤ 8 which are complemented by ground state QMC results within 0.5 ≤ r s ≤ 20 [104].The GDSMFB parametrization also utilizes synthetic data in the range 0.0625 ≤ Θ ≤ 0.25 (inaccessible to simulations owing to the prevalence of the fermion sign problem [50]) which are constructed by combining the ground state QMC results with a small STLS-based temperature correction.For the paramagnetic case, the quasi-exact QMC results employed for the fit concern 58 warm dense UEF state points and 7 ground state UEF state points [54,60,104].
VS-generated structural properties will also be compared with the highly accurate results of the effective static approximation (ESA).The ESA approach is a recently developed semi-empirical scheme of the dielectric formalism that targets the WDM regime (strict applicability limits within 0.7 ≤ r s ≤ 20 and 0 ≤ Θ ≤ 4) [62,63].More specifically, the ESA constructs an analytical SLFC based on the exact long wavelength limit as determined from the CSR [105], on the exact large wavenumber limit as determined from the cusp relation within the static assumption [106] and on a neural net representation of QMC results at the intermediate wavenumbers [107].The quasi-exact incorporation of the asymptotic limits is facilitated by the availability of the aforementioned GDSMFB f xc parametrization and the availability of an on-top PCF (value at zero distance) parametrization [43,62,104].ESA yields remarkable results for most physical quantities in the WDM regime, not only static (such as the SSF and the electronically screened ion potential) but even dynamic (such as the dynamic structure factor and the stopping power) [63].In addition, the computational cost of the ESA scheme is negligible, since it is comparable to that of the RPA scheme owing to the fact that the SLFC is introduced as function of the wavenumber and not as a functional of the SSF.Nevertheless, given its QMC-based construction, the ESA scheme cannot provide microscopic insights on the interplay between thermal excitations, strong correlations and quantum effects.

A. Self-consistency parameter
The complications in the numerical solution of the finite temperature VS scheme primarily originate from the necessity for a double convergence of the inner and outer loops.The availability of a closed-form expression for the self-consistency parameter α(r s , Θ) is equivalent to cancelling the outer loop and considerably simplifies the VS algorithm.The dependence of the self-consistency parameter α on the UEF thermodynamic variables (r s , Θ) is explored in Fig. 2. The following conclusions are due: For any Θ, α is a monotonically increasing function of the quantum coupling parameter that rises abruptly until it reaches an extended plateau.For any r s , α is a monotonically decreasing function of the degeneracy parameter up to Θ ∼ 4, with the dependence becoming non-monotonic at higher temperatures.In contrast to the ground state limit, α is not bound within 1/2 ≤ α 1 [70] and even obtains negative values at very high densities.In contrast to the ground state limit, α ≃ 2/3 [70] does not constitute an accurate approximation, unless Θ 0.5 and r s 1.
The analytic expression for α(r s , Θ) is obtained by solving the correct finite temperature VS scheme for 2250 UEF state points spanning the entire warm dense matter regime.To be more specific, 150 values of the quantum coupling parameter, r s ∈ [0.1, 15] with a step of 0.1, are combined with 15 values of the degeneracy parameter Θ ∈ [0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8].In order to take advantage of the monotonic r s dependence, the isothermal values of α are first least square fitted to The resulting 15 data points for the four Θ−dependent functions κ α = {γ α , δ α , ζ α , η α } are then fitted to even Padé approximants of the order 6/6, The 28 Padé coefficients are provided in Table I together with the mean absolute relative error of each Padé fit.The Padé approximants are plotted in the inset of Fig. 2.
The overall mean absolute relative error of the total fit is 1.4%.The errors stem mostly from the plateau region, where there are weak α fluctuations that are smoothed out by the fitting function, see Fig. 2.
With the aid of this α(r s , Θ) expression, the VS scheme can now be solved by algorithms that simultaneously solve the STLS scheme for nine UEF state points.Given the ∼ 1% fitting error and the 10 −3 convergence criterion of the outer loop, it can be roughly stated that such STLS algorithms will solve the VS scheme with approximately an order of accuracy less than our pure VS algorithm.

B. Exchange correlation free energy
As per usual, the VS SSF is employed for the determination of the UEF interaction energy u int via the fundamental expression of Eq. (12).Then, the UEF interaction energy is employed for the computation of the UEF exchange correlation free energy f xc via the integral version of the adiabatic connection formula of Eq.( 13).Different thermodynamic properties can then be computed from well known expressions [30].The VS-generated u int and f xc have the standard UEF dependencies [31], as seen in Figs.3,4.In what follows, we exclusively focus on f xc .
The f xc (r s , Θ) parametrization is based on the solution of the finite temperature VS scheme for 2250 UEF states that cover the entire WDM regime.Again, 150 values of the quantum coupling parameter, r s ∈ [0.1, 15] with a step of 0.1, and 15 values of the degeneracy parameter Θ ∈ [0.25, 0.5, 0.75, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8] are probed.Our f xc parametrization is inspired by the IIT parametrization originally proposed for u int [31,60,88].The isothermal values are first least square fitted to where α HF (Θ) corresponds to the exact Hartree-Fock coefficient of the non-interacting UEF exchange free energy f (r s , Θ) = −α HF (Θ)/r s that is rigorously defined by the double integral presentation and has been approximated by the very accurate Perrot and Dharma-wardana parametrization [103] α The resulting 15 data points for the four Θ−dependent functions κ f = {γ f , δ f , ζ f , η f } are then fitted to even Padé approximants of the order 6/6, The 28 Padé coefficients are provided in Table II together with the mean absolute relative error of each Padé fit.The overall mean absolute relative error of the total fit is merely 0.36%.
A direct comparison between the VS-generated f xc and STLS-generated f xc can be performed neither graphically (small differences) nor tabularly (too many data points).An indirect comparison will be facilitated by the very accurate GDSMFB parametrization.The absolute relative deviations of the VS and STLS f xc predictions from the GDSMFB f xc value as a function of r s have been plotted in Fig. 5 for 3 Θ values.Raw VS and STLS data have been preferred over the analytic fits for the purposes of com-   II.parison.Thus, the smoothness of the curves is a direct manifestation of the stability and accuracy of our algorithm.This analysis has been pursued for all 15 probed Θ.In the entire WDM regime, the STLS predictions f xc are more accurate.The reason for the thermodynamic superiority of the STLS has been earlier discussed in the literature and will be elucidated in Sec.IV C. As we shall see, this does not translate to structural superiority.

C. Static structure factor
A comparison between the SSF predictions of the STLS, VS and ESA schemes is featured in Fig. 6. (i) The STLSgenerated SSF is consistently characterized by an overestimation at low wavenumbers followed by an underestimation at intermediate and high wavenumbers, with the crossover point lying within 1.5k F − 2.0k F .On the other hand, the VS-generated SSF is consistently characterized by an underestimation at low and intermediate wavenumbers followed by a slight overestimation at high wavenumbers, with the crossover point lying in the vicinity of 3k F .Since the integrated SSF yields the interaction energy, see Eq.( 12), this observation suggests a favorable cancellation of errors for the STLS scheme and weak cancellation of errors for the VS scheme, which explains the thermodynamic supremacy of the STLS scheme.(ii) At sufficiently high densities for which the SSF maximum is either extremely shallow or even nonexistent, the STLS scheme is consistently more accurate than the VS scheme in nearly the entire wavenumber range.The ESA SSF lies in between the STLS SSF and VS SSF at small wavenumbers, while S ESA > S STLS > S VS is valid at intermediate and large wavenumbers.(iii) At stronger coupling parameters for which the SSF maximum is apparent but remains shallow, the accuracy of the STLS and VS schemes becomes comparable.The ESA-generated SSF is nearly overlapping with the VS SSF at small wavenumbers, in better agreement with the STLS SSF at intermediate wavenumbers and slightly closer to the VS SSF near the large wavenumber limit.It should also be noted that the VS scheme leads to an extremely shallow SSF maximum, while the STLS SSF remains bounded below unity.(iv) Near the boundary between the WDM and the strong coupling regimes where the SSF maximum is pronounced and the secondary extrema begin to emerge, the VS scheme is superior with the exception of a narrow 1.5k F − 2.0k F interval located in the neighborhood of the ESA-STLS crossover.More specifically, the VSgenerated SSF is nearly indistinguishable from the ESAgenerated SSF within the small wavenumber range and it is characterized by a well-developed maximum, albeit weaker and broader than the nearly exact prediction of the ESA scheme.Meanwhile, the STLS-generated SSF remains bounded below unity.At strong coupling, as far as static structural properties are concerned, the improved performance of the VS scheme compared to the STLS scheme is expected from a general equilibrium statistical mechanics perspective [108].Let us first recall an exact relationship between successive order reduced density matrices which reads as [95] ̺ where ̺ s (r s ; R s ) = r s |̺ s |R s is the reduced s-particle density matrix within the coordinate representation and In terms of reduced s-particle densities n s (r s ) = ̺ s (r s ; r s ), this relationship remains intact, i.e., n s (r s ) = 1 N − s d 3 r s+1 n s+1 (r s+1 ) .
Within the thermodynamic limit s ≪ N → ∞ and for homogeneous systems, in terms of reduced s-particle correlation functions g s (r s ) = n s (r s )/n s , this relationship becomes [108,109] g s (r s ) = 1 V d 3 r s+1 g s+1 (r s+1 ) .
Setting s = 2, the above implies a connection between the (isothermal) density derivative of the pair correlation function and the triplet correlation function.Such connections have been long exploited in the theory of classical liquids [110][111][112].Thus, it is concluded that the indirect inclusion of ternary correlations in the VS scheme is responsible for the improved description of equal-time density correlations compared to the STLS scheme.
The dependence of the VS-generated SSFs on the UEF thermodynamic variables (r s , Θ) is explored in Fig. 7.At constant Θ and for increasing r s ; the small wavenumber behavior of the SSF becomes less steep.When Θ = 0.5, a shallow maximum appears for r s 8 which becomes sharper as the quantum coupling parameter further increases.For all probed isothermal state points, the large wavenumber limit of unity is always reached before 4.5k F .Essentially, the Θ−dependency follows the same reasoning as the r s −dependency, provided that an effective coupling parameter Γ eff is introduced which considers the contribution of exchange effects on an equal footing with thermal effects.In view of the classical coupling parameter, one now defines . For constant r s and increasing Θ, the small wavenumber behavior of the SSF becomes steeper.When r s = 10, a shallow maximum appears for Θ 2 that becomes sharper as the degeneracy parameter further decreases.For all probed isochoric state points, the large wavenumber limit of unity is reached before 5k F .

D. Static density response function
The linear SDR function is directly obtained from PIMC results for the imaginary time density-density correlation function via the imaginary-time version of the quantum fluctuation-dissipation theorem [113][114][115].The SDR function is known to be more sensitive to quantum effects compared to the SSF [115].Nevertheless, no surprises are expected in the course of the comparison between the VS scheme and the STLS scheme, given the fact that both schemes treat quantum effects at the level of the RPA.It is pointed out that the large wavenumber limit of the SDR is dictated by the ideal SDR, i.e., [116,117] and that the long wavelength limit of the SDR is dictated by the perfect screening, i.e. [116,117] A comparison of the SDR predictions of the STLS, VS and ESA schemes is featured in Fig. 8.The asymptotic  limits are indicated by the dashed lines.(i) In accordance the SSF observations, in the small wavenumber region, the quasi-exact ESA SDR lies in-between the VSgenerated and STLS-generated SDRs at high densities.When the quantum coupling parameter r s increases, the VS-generated SDR nearly overlaps with the ESA SDR for k 1.5k F , while the STLS-generated SDR exhibits observable deviations from the ESA SDR up to k/k F ∼ 0.7.It is noted that the perfect screening limit is reached for smaller wavenumbers, as the coupling becomes stronger.(ii) Again in accordance with the SSF observations, in the large wavenumber region, the quasi-exact ESA SDRs lie closer to the VS-generated SDRs irrespective of the quantum coupling parameter.Unsurprisingly, as the coupling becomes stronger, the convergence to the ideal SDR is shifted to higher wavenumbers.(iii) Irrespective of the quantum coupling parameter, the STLS scheme yields more accurate predictions for the magnitude of the SDR minimum, while the VS scheme yields more accurate predictions for the position of the SDR minimum.In particular, as the coupling increases; the STLS predictions for the magnitude of the minimum become progressively worse while the offset concerning the position of the minimum remains nearly constant, the VS predictions for the position of the minimum remain very accurate while the discrepancy concerning the magnitude of the minimum remains nearly constant.(iv) Overall, also in line with the SSF observations, the STLS-generated SDRs are more accurate for weak coupling and the VS-generated SDRs are more accurate for moderate coupling.
The dependence of the VS-generated SDRs on the UEF thermodynamic variables (r s , Θ) is explored in Fig. 9.At constant Θ and for increasing r s ; the short wavenumber behavior of the SDR is characterized by a very strong dependence as expected from the normalized unit version lim x→0 χ(x)/(nβ) = −(3π/8λ)(Θ/r s )x 2 of the perfect screening asymptote.On the other hand, all the SDRs exhibit the same large wavenumber behavior given the r s −independence of the ideal SDR with the r s = 15 curve FIG.8: The static density response of the paramagnetic uniform electron fluid in the warm dense matter regime, as predicted by the VS scheme (red solid line), the STLS scheme (black solid line) and the quasi-exact ESA scheme (blue solid line).The ideal non-interacting limit (green dashed line) and the perfect screening limit (orange dashed line) are also indicated.Results for Θ = 1.0 and for rs = 5 (left), rs = 10 (center), rs = 15 (right).attaining its limiting behavior at the highest wavenumber.The position of the omnipresent SDR minimum gets gradually shifted to larger wavenumbers, as the quantum coupling parameter increases.The magnitude of the SDR minimum is monotonically decreasing with r s up to the emergence of the maximum of the SSF.At even higher coupling, the magnitude of the SDR minimum begins to increase with r s [57,118].At constant r s , there is a strong dependence of both asymptotic limits on Θ.It is emphasized that the position of the SDR minimum weakly depends on the degeneracy parameter, whereas the SDR minimum becomes broader as Θ increases.

E. Pair correlation function
The PCF is obtained from the SSF by the inverse Fourier transform of the fundamental microscopic expression S(k) = 1 + nH(k) [97].For isotropic systems such as the UEF, use of spherical coordinates and normalized units leads to It is emphasized that specialized quadrature algorithms are necessary for accurate numerical evaluation of this integral at short distances owing to the rapid oscillatory nature of the sin (xy) factor of the integrand [119][120][121].
It is also pointed out that the truncation of the upper integration limit should be handled with care; the SSF reaches its unity asymptotic limit around y ∼ 5 and thus it is possible that noise of alternating sign is solely integrated at large wavenumbers for which the SSF fluctuates around unity with magnitudes that are lower than the accuracy of the algorithm employed for the solution of the dielectric scheme of interest.
A comparison of the PCF predictions of the STLS, VS and ESA schemes is featured in Fig. 10.It is preferable to distinguish between the near-contact short distance region (which is primarily shaped by the competition between quantum effects that allow the overlap of electrons of opposite spin [122] and strong interaction effects that favor the formation of correlation voids [123,124]) and the intermediate-long distance region that consists of the first coordination shell (which is primarily shaped by the relative strength of the correlations [97,108]) together with the asymptotic decay range (where the PCF reaches its long distance limit of unity [97]).(i) Let us first consider the combined intermediate-long distance region rk F 1. At high densities, the STLS-generated and VS-generated PCFs have comparable accuracy and exhibit small deviations from the quasi-exact ESA PCF that are spread over all distances rk F 4. At stronger coupling, the ESA PCF has a shallow maximum whose position is accurately predicted by the VS scheme, unlike its magnitude which is underestimated by the VS scheme.On the other hand, the STLS scheme leads to an even shallower PCF maximum that is displaced towards larger distances.Thus, STLS-generated PCFs may attain their asymptotic limit at larger distances than VS-and ESAgenerated PCFs.(ii) Let us next proceed with the discussion of the short distance region rk F 1. The STLS and the VS schemes lead to an un-physical negative PCF region near the origin (r = 0).This is a well-known common pathological feature of all non-empirical schemes of the dielectric formalism that originates from the approximate treatment of quantum effects [31,125,126].It is particularly prominent in semi-classical schemes, but it is also present in pure quantum schemes as well as in quantized schemes based on integral equation theories [79,81].On the other hand, mapping approaches that are based integral equation theories for an effective pair interaction potential within a specific classical-quantum state correspondence rule by default lead to a non-negative PCF [125,126].As expected based on their semi-classical classification, both the VS and STLS schemes are characterized by extended negative PCF regions.For most state points, the absolute value of the integrated negative region of the VS-generated PCF is larger than the respective absolute value of the STLS-generated PCF.(iii) Let us also focus on the contact (r = 0).The larger integrated negative PCF region of the VS scheme does not necessarily translate to a more negative on-top PCF, g(r = 0), owing to the existence of a still unphysical negative global minimum of the PCF.The latter is observed for both the VS and the STLS schemes from moderate coupling.More specifically, even though the VS minimum is deeper, the larger VS slope at the contact-side of the minimum yields less negative (or even positive) on-top PCF values compared to the STLS scheme.(iv) For completeness, it is noted that the semi-empirical nature of the ESA scheme allows it to be an exception to the negative PCF pathology of dielectric schemes.It incorporates accurate QMC-based on-top PCF values and, for nearly all state points within its applicability range, yields a PCF that is monotonic prior to the first coordination shell maximum implying that ESA-generated PCFs are always positive.
The dependence of the VS-generated PCFs on the UEF thermodynamic variables (r s , Θ) is explored in Fig. 11.Inside the first coordination shell, the r s −dependence and the Θ−dependence of the PCF can be fully understood on the basis of the earlier SSF discussion at the last paragraph of Sec.IV C. Concerning short distances; the PCF is always positive and monotonic at high electron densities, the onset (but not necessarily the extent) of the negative PCF region is observed to monotonically increase as r s rises and monotonically increase as Θ decreases up to 0.5 Θ 1, the negative PCF minimum is not present for low effective coupling parameters (high Θ and/or low r s ) and the on-top PCF values are positive either at high densities (non-negative PCF) or for large effective coupling parameters (negative PCF region with a minimum).

F. Static local field correction
Provided that the local field correction is static, its large wavenumber limit is algebraically connected to the ontop PCF.After the substitution of G(k, ω) ≡ G(k) in the constitutive relation for the density response function Eq.( 4), the combination with the infinite Matsubara summation of Eq.( 5) and the utilization of the fundamental PCF-SSF connection, an intricate asymptotic analysis yields the general result [68,88] Simultaneously, the cusp condition is exactly valid which states that the derivative of the PCF logarithm at the origin equals the inverse Bohr radius [106,127,128], i.e., ∂g(r) Combining the above, one obtains the exact cusp relation that reads as [68,88] This can be considered as a self-consistency condition [68,88] since it solely originates from the first two exact building blocks of the dielectric formalism, see Eqs. (4,5), and needs to be independently satisfied by the third approximate building block, see the closure functional of Eq.( 6) without the frequency dependence.
It can be easily shown that the STLS scheme satisfies the cusp self-consistency relation.The large wavenumber limit x → ∞ of the bracketed factor of the STLS SLFC integrand, see Eq. (11), is found to be equal to 2, after a Taylor expansion with respect to y/x → 0. This leads to where the state point dependence has been suppressed.The short distance limit x → 0 of the fundamental PCF-SSF connection, see Eq.( 31), since sin (xy) ≃ xy, yields Thus, combining the above, one indeed obtains that It is rather straightforward to prove that the VS scheme cannot satisfy the cusp self-consistency relation.In view of the differential connection between the VS and STLS SLFCs, see Eq. (7) or Eq.( 10), and courtesy of the asymptotic x(∂G STLS (x)/∂x)| x→∞ = 0, we directly obtain [105] lim From general physical arguments that are based on the non-interacting limit and the strongly interacting limit, it can be deduced that the thermodynamic derivatives of the on-top PCF cannot be identically zero.In addition, based on the extended discussion featured at the last paragraph of Sec.IV E, it can be safely concluded that the thermodynamic derivatives of the on-top PCF do not provide negligible contributions to Eq. (36).Therefore, the VS scheme does not even roughly comply with the cusp self-consistency relation.In a sense, enforcement of small wavenumber consistency (satisfaction of isothermal compressibility sum rule) has the undesired side-effect of large wavenumber inconsistency (violation of cusp condition).
A comparison of the SLFC predictions of the STLS, VS and ESA schemes is featured in Fig. 12.It is pointed out that the ESA scheme nearly exactly satisfies the cusp selfconsistency relation and the CSR expression by construction [63].To prevent possible misconceptions, it is also emphasized that the ESA scheme features a frequencyaveraged LFC that does not correspond to the exact static limit of the dynamic LFC as that would emerge from QMC simulations [62,63].(i) Concerning the SLFC long wavelength limit, as expected from its direct connection with the CSR, VS-generated SLFCs are more accurate than STLS-generated SLFCs regardless of the UEF state point.Furthermore, as the quantum coupling parameter increases, the SLFC result exhibits less deviations from the accurate ESA result.This is consistent with the improved agreement of the VS exchange correlation free energy with the GDSMFB parametrization as r s increases.(ii) Concerning the large wavenumber limit of the SLFC, as expected from its connection with the cusp condition, the VS-generated SLFCs are less accurate than the STLS-generated SLFCs regardless of the UEF state point of interest.Regarding the STLS scheme and the ESA scheme, the cusp self-consistency relation also provides an independent check of the numerical accuracy of the respective algorithms; the STLS SLFC large wavenumber limit is indeed larger than unity in accordance with the negative on-top PCF, whereas the ESA SLFC large wavenumber limit is indeed smaller (or nearly equal to) than unity in accordance with the non-negative on-top PCF.(iii) Concerning the intermediate wavenumber region, the deviations from the ESA SLFC are always rather large, but it can be stated that the STLS SLFC is more accurate at high densities and that the VS SLFC becomes more accurate as the coupling increases.The STLS SLFC is characterized by a shallow maximum only near the boundary between the WDM and the strongly coupled regimes, while the VS SLFC only lacks a maximum at high densities.Yet, even at the largest coupling parameters probed, the well-formed VS SLFC maximum could not match the magnitude, position and width of the ESA SLFC maximum.(iv) In contrast to the sharp ESA SLFC transition between the intermediate wavenumber region and the asymptotic value, in case a SLFC maximum is present, the respective VS and STLS SLFC tran- sitions are slow and could even extend over several Fermi wavenumbers.
The dependence of VS-generated SLFCs on the UEF thermodynamic variables (r s , Θ) is explored in Fig. 13.The strong dependence of the large wavenumber limit on (r s , Θ) reflects the strong dependence of the on-top PCF on (r s , Θ), in view of Eq. (36).The aforementioned correlation between the presence of a SLFC maximum and the prolonged asymptotic transition is also apparent.Finally, the weak dependence of the long wavelength limit (r s , Θ), except from high Θ and/or low r s , reflects the relatively weak dependence of the exchange-correlation free energy on (r s , Θ), in view of Eq.( 9).

G. Comparison with the VS-SD scheme
Owing to the erroneous chain rule application in the density derivative, the Sjostrom and Dufty formulation of the finite temperature VS scheme (VS-SD scheme) [71] is considerably simpler than the present corrected formulation of the finite temperature VS scheme.First and foremost, when utilizing second order finite difference schemes to approximate the thermodynamic derivatives, the minimum VS-SD computational stencil features three state points, while the minimum corrected VS computational stencil features nine state points, compare Eqs.(10,16).As a consequence, the inner loop of the VS-SD algorithm needs to be simultaneously run for three state points, while the inner loop of the VS algorithm needs to be simultaneously run for nine state points.Thus, the VS-SD computational cost is at least three times smaller.In addition, the self-consistency parameter equation of the VS-SD scheme involves considerably less thermodynamic derivatives than the self-consistency parameter equation of the VS scheme, compare Eqs.(15,17).This directly suggests less complexity and indirectly implies less computational cost due to the faster convergence of the outer loop.
In this subsection, we compare the predictions of the (corrected) VS scheme with the predictions of the VS-SD scheme.Such a comparison allows us to confirm that the corrected treatment has an meaningful impact on the structural properties of the UEF and is not merely an added complexity.Such a comparison also allows us to identify phase diagram regions wherein the VS scheme is indistinguishable with the VS-SD scheme and thus the simpler version can be employed without introducing errors.After a simple inspection of the respective SLFC closures and respective self-consistency parameter equations, it can be deduced that the two schemes will have maximum differences at high Θ and/or small r s and will nearly overlap at low Θ and/or large r s .On the other hand, the structural properties are not very sensitive to typical values of the quantum coupling parameter that are relevant for WDM applications (r s ∼ 5), see Figs. 7,9.
A comparison of the SSF and PCF predictions of the STLS, VS and VS-SD schemes for Θ = 5 and r s = 5, 10 is featured in Fig. 14.In line with our expectations, the VS-SD scheme leads to unique predictions for both structural properties.In particular, as discerned in the main figures, the VS-generated SSF consistently lies in-between the STLS (upper bound) and VS-SD (lower bound) predictions.Furthermore, as illustrated in the insets, there are substantial PCF differences between the VS and the VS-SD schemes that are mainly confined to short distances.Finally, the inset of the right panel even reveals that the unphysical negative PCF region near the origin can be monotonic for the VS scheme and non-monotonic for the VS-SD scheme.A comparison of the SDR predictions of the STLS, the VS and the VS-SD schemes for Θ = 5 and r s = 5, 10 is featured in Fig. 15.Unsurprisingly, the VS-SD scheme yields distinct predictions for the static limit of the linear density response function.For both coupling parameters, the VS-generated SDR lies in-between the STLS (lower bound) and the VS-SD (upper bound) predictions.In general, this is a pattern that persists within the WDM regime.Therefore, it can be stated that the correct finite temperature formulation of the VS scheme leads to predictions that are closer to the STLS scheme compared to the predictions of the VS-SD scheme.Another comparison of the SDR predictions of the STLS, VS and VS-SD schemes now for Θ = 0.5 and r s = 5, 10 is featured in Fig. 16.The VS and VS-SD results for the SDR are truly indistinguishable, overlapping in the entire wavenumber range.The same applies for the PCF, SSF and SLFC.

A. Summary
In this work, the finite temperature formulation of the Vashishta-Singwi dielectric formalism scheme for the uniform electron fluid was reported correcting for a thermodynamic derivative error that was present in the earlier formulation of Sjostrom and Dufty.This is a truly self-consistent version, where the density expansion parameter α is not locked to an approximate state independent value but is determined by enforcing the compressibility sum rule at each thermodynamic point.An efficient computational scheme was devised for the solution of the VS scheme that was presented in full detail.It features an inner loop that simultaneously solves the scheme for nine state points, which form a minimum extent computational stencil when all finite difference approximations of the thermodynamic derivatives are of the second order.It also features an outer loop that determines the self-consistency parameter α through a non-linear equation of the form α where f 1 contains all possible first-and second-order thermodynamic derivatives of f xc and f 2 contains all possible first-order thermodynamic derivatives of u int , which is solved with the secant method.The VS algorithm is implemented in a home-made hybrid code, where C++ is used for the backend and python is used for the frontend.
The finite temperature VS scheme was numerically solved for 2250 state points of the paramagnetic uniform electron fluid that cover the entire warm dense matter regime.First, this extended parametric study facilitated the extraction of an accurate closed-form expression for the self-consistency parameter α(r s , Θ) whose availability removes the necessity for an outer loop and permits the use of STLS-like algorithms for the solution of the VS scheme.Furthermore, this systematic investigation provided sufficient data for the determination of an accurate closed-form expression for the exchange correlation free energy following the lines of state-of-the-art parametrizations, which adequately represents the thermodynamic predictions of the VS scheme.More important, this allowed a comprehensive thermodynamic and structural comparison with the respective predictions of the ubiquitous STLS scheme as well as the near-exact thermodynamic results of the QMC-based GDSMFB equation of state and the near-exact structural results of the QMCbased effective static approximation.
The primary conclusions drawn from the aforementioned comparison are summarized in the following: (i) In the entire WDM regime, the thermodynamic predictions of the STLS scheme are more accurate than the thermodynamic predictions of the VS scheme benefitting from a very favorable cancellation of errors in the static structure factor integral for the interaction energy.(ii) In the high density range of the WDM regime, r s 5 at the Fermi temperature, the STLS static structure factor is consistently more accurate than the VS static structure factor.The situation reverses at the heart of the WDM regime and towards its boundary with the strongly coupled regime, especially at the long wavelength limit and at the vicinity of the static structure factor maximum.This observation can be explained by the fact that the VS scheme indirectly includes ternary correlations through the density derivative of the pair correlation function.(iii) As far as the static density response is concerned, both long wavelength and large wavenum-ber asymptotic limits are automatically satisfied by both schemes, but the STLS scheme yields more accurate predictions for the magnitude of the omnipresent minimum (which progressively worsen as the coupling increases), while the VS scheme yields more accurate predictions for the position of the omnipresent minimum.(iv) The structural superiority of the VS scheme at the strongly coupled region of the WDM regime and the structural superiority of the STLS scheme at the high density region of the WDM regime are also confirmed from the pair correlation function and the static local field correction.(v) An unphysical negative region in the pair correlation function near contact is present for both schemes.This does not necessarily translate to negative on-top values, because, at strong coupling, the pair correlation function can exhibit a negative minimum.For most state points, the absolute value of the integrated negative region is larger in the VS scheme than in the STLS scheme.This short distance drawback of the VS scheme most probably originates from the violation of the cusp condition, which is exactly satisfied by the STLS scheme.
Finally, our correct formulation of the VS scheme collapses to the earlier Sjostrom and Dufty formulation of the VS scheme at low values of the degeneracy parameter and/or large values of the quantum coupling parameter.This was expected by the nature of the thermodynamic derivative error of the Sjostrom and Dufty analysis and has been confirmed numerically.Thus, there is still some utility in the less computationally costly and easier to numerically implement Sjostrom and Dufty formulation.

B. Future work
One of the main drawbacks of the finite temperature and ground state versions of the VS scheme concerns its semiclassical nature, i.e., the treatment of quantum mechanical effects at the level of the random phase approximation.However, the standard VS closure can be employed to truncate the first member of the quantum BBGKY hierarchy, which would allow the incorporation of quantum effects beyond the Vlasov-Lindhard ideal response substitution that would then be represented by a dynamic LFC (DLFC).We shall refer to this dielectric scheme as the quantum VS (qVS).The standard perturbative analysis of the quantum kinetic equation within the Wigner representation leads to the qVS DLFC that reads as In these expressions; G qSTLS (k, ω) denotes the DLFC of the quantum STLS (qSTLS) scheme, which emerges by truncating the first member of the quantum BBGKY hierarchy with the standard STLS closure [72][73][74], while χ 0 (k, k ′ , ω) denotes the three-argument ideal density response function which collapses to the ideal (Lindhard) density response function when k ′ = k, χ 0 (k, k, ω) = χ 0 (k, ω) [72][73][74].Note that the correspondence between the qVS DLFC, qSTLS DLFC is exactly the same with the correspondence between the VS SLFC, STLS SLFC, compare the above with Eqs.(7,8).It is also important to emphasize that the qVS DLFC can be directly obtained from the VS SLFC after the substitution (k • k ′ )/k 2 → χ 0 (k, k ′ , ω)/χ 0 (k, ω) within the wavenumber integrand; this substitution rule was first observed in Ref. [81] where it was utilized to rapidly quantize semi-classical nonperturbative schemes [77,79,80].To our knowledge, the ground state version of the qVS scheme has only been considered with a locked parameter α = 1/2 in the literature [73,129], which implies that the CSR rule is slightly violated.In the fully self-consistent case, the transition from the present VS scheme to the envisaged qVS scheme will be accompanied by a (manageable) blow-up in the computational cost of the inner loop.In particular, the dynamic nature of the LFC implies that the qVS closure equation depends on the Matsubara frequencies, while the presence of the three-argument ideal density response implies that the qVS closure equation is now described by a triple (instead of a single) integral.
Another main drawback of the finite temperature version of the VS scheme concerns the violation of the cusp self-consistency relation, see Eqs. (34,36), which should be partly responsible for the extended unphysical negative pair correlation function region in the proximity of contact.The importance of simultaneously satisfying the compressibility sum rule and the cusp relation has been manifested by the semi-empirical ESA scheme [62,63].Thus, it is worth trying to enforce both self-consistency relations in a dielectric scheme that does not utilize QMC results.With this objective in mind, we bring forth the possibility for an extended VS (eVS) scheme of semiclassical nature that features two self-consistency parameters.The truncation at the first member of the classical BBGKY hierarchy will be based on the eVS closure f 2 (r, p, r ′ , p ′ , t) = f (r, p, t)f (r ′ , p ′ , t)g(r, r ′ , t) where the non-equilibrium pair correlation function is now given by g(r, r ′ , t) = g eq (|r − r ′ |; n, T ) + α [δn(r, t) + δn(r ′ , t)] ∂g eq (|r − r ′ |; n, T ) ∂n + γ [δT (r, t) + δT (r ′ , t)] ∂g eq (|r − r ′ |; n, T ) ∂T .
The eVS closure is characterized by the addition of temperature perturbations on an equal footing with the density perturbations.This addition will have profound effects in the mathematical treatment of the linearized kinetic equation and thus in the inner loop of the eVS scheme's algorithm.In particular, the linear temperature response function χ T (k, ω) = δT (k, ω)/δU ext (k, ω) would need to be studied in parallel with the linear density response function χ(k, ω) = δn(k, ω)/δU ext (k, ω), which would require that the zeroth and second moments of the linearized kinetic equation are considered simultaneously, forming a set of equations that allows the determination of both χ T (k, ω) and χ(k, ω).An ideal temperature response, i.e. the second moment equivalent of the ideal (Lindhard) density response, would also naturally emerge.The two self-consistency parameters α, γ would be determined in the outer loop of the algorithm by a 2×2 set of equations that is formed when simultaneously imposing the cusp relation and compressibility sum rule.When γ = 0, it is evident that the envisaged eVS scheme collapses to the standard VS scheme.Finally, it is also worth noting that the coupling between density and temperature fluctuations in the eVS scheme would lead to a density response function expression that does not comply with the general result of the polarization potential approach [96], see Eq.( 4).On the other hand, since the zero frequency moment sum rule and the quantum fluctuation-dissipation theorem will remain intact, this coupling between density and temperature fluctuations will have no effect on the infinite Matsubara series expression that connects the static structure factor with the density response function, see Eq.( 5), which constitutes the backbone of all schemes of the dielectric formalism.
p and use of the functional derivative definition of the linear density response function χ(k, ω) = δn k,ω /δU ext k,ω yield the VS density response function

FIG. 5 :
FIG.5:The absolute relative deviations of the VS fxc predictions (blue) and the STLS fxc predictions (orange) from the very accurate fxc value of the GDSMFB parametrization as a function of rs for Θ = 0.5 (top), Θ = 2.0 (centre) and Θ = 4.0 (bottom).In the case of Θ = 0.5, the abrupt monotonicity change for the STLS scheme that occurs at rs ≃ 5.5 corresponds to the switch of the sign of the difference.

FIG. 6 :
FIG.6:The static structure factor of the paramagnetic uniform electron fluid in the warm dense matter regime, as predicted by the VS scheme (red), the STLS scheme (black) and the quasi-exact ESA scheme (blue).Results for Θ = 1.0 and for rs = 5 (left), rs = 10 (center), rs = 15 (right).

FIG. 14 :
FIG.14: Static structure factor (main) and pair correlation function (inset) of the paramagnetic uniform electron fluid in the warm dense matter regime, as predicted by the VS-SD scheme (blue), VS scheme (red) and STLS scheme (black).Results for Θ = 5.0 and rs = 5 (left), rs = 10 (right).

TABLE I :
(24,25)efficients for the self-consistency parameter of the VS scheme, see Eqs.(24,25).The mean absolute relative error of each Padé fit is reported in the last row.

TABLE II :
(26,27)efficients for the exchange correlation free energy of the VS scheme, see Eqs.(26,27).The mean absolute relative error of each Padé fit is reported in the last row.