Effective electronic forces and potentials from ab initio path integral Monte Carlo simulations

The rigorous description of correlated quantum many-body systems constitutes one of the most challenging tasks in contemporary physics and related disciplines. In this context, a particularly useful tool is the concept of effective pair potentials that take into account the effects of the complex many-body medium consistently. In this work, we present extensive, highly accurate \emph{ab initio} path integral Monte Carlo (PIMC) results for the effective interaction and the effective force between two electrons in the presence of the uniform electron gas (UEG). This gives us a direct insight into finite-size effects, thereby opening up the possibility for novel domain decompositions and methodological advances. In addition, we present unassailable numerical proof for an effective attraction between two electrons under moderate coupling conditions, without the mediation of an underlying ionic structure. Finally, we compare our exact PIMC results to effective potentials from linear-response theory, and we demonstrate their usefulness for the description of the dynamic structure factor. All PIMC results are made freely available online and can be used as a thorough benchmark for new developments and approximations.


I. INTRODUCTION
The accurate description of interacting quantum many-body systems [1] constitutes one of the most active frontiers in a number of disciplines, such as theoretical physics, quantum chemistry, and material science. The first groundbreaking insights into their collective behaviour have been obtained on the basis of well established theoretical approximations, such as the model of the dilute Bose gas [2] and the weakly coupled uniform electron gas (UEG) that can be described within the random phase approximation (RPA) [3]. Naturally, the situation becomes incomparably more complicated and interesting in the regime of moderate to strong coupling strength, where the aforementioned approximations break down. These theoretical challenges have facilitated the development of sophisticated numerical methods that are capable of capturing the rich interplay of quantum effects such as diffraction and quantum statistics with correlation effects as well as thermal excitations. A case in point is the ab initio path integral Monte Carlo (PIMC) method [4][5][6], which is a finite-temperature implementation of the quantum Monte Carlo (QMC) paradigm [7]. Specifically, the PIMC approach has been of pivotal importance for our understanding of Bose liquids such as ultracold 4 He [4,8,9]. In fact, state-of-the-art implementations [10,11] allow for the quasi-exact simulation of up to N ∼ 10 3 − 10 4 particles, thereby giving us microscopic insights into important physical effects such as superfluidity [4,12,13] and Bose-Einstein-condensation [14][15][16].
Unfortunately, QMC simulations of quantum degenerate Fermi systems, such as the UEG [17] or ultracold Fermi atoms like 3 He [18] are afflicted with the notorious * t.dornheim@hzdr.de fermion sign problem [19,20], which leads to an exponential increase in compute time with increasing the system size or decreasing the temperature. Therefore, QMC simulations of many-electron systems are usually restricted to rather small systems (typically N ∼ O 10 1 − 10 2 ) even using state-of-the-art algorithms on modern highperformance computing clusters. At the same time, we note that these simulations are often remarkably accurate with as few as N = 14 electrons [17,[21][22][23][24][25] and, in combination with dielectric methods (like the RPA) that become exact in the long-range limit, provide an adequate description over all length scales. This is particularly remarkable in the case of electrons, as the longrange Coulomb interaction is known to cause difficulties and divergencies in a number of approaches [26,27].
In this context, a particularly important concept is the notion of effective pair interactions and forces, which constitute a true mile stone in the description of quantum many-body systems. In general, we distinguish three distinct situations [26,28]: 1. The effective interaction between two test charges in a medium. The most prominent example is given by effective ion-ion interactions that are screened due to the presence of a medium consisting of free, nonideal electrons [29][30][31][32][33][34].
2. The effective interaction between a test change and the surrounding medium, such as an ion embedded in a uniform electron gas [35,36].
3. The effective interaction between two particles that are part of the surrounding medium. For example, the effective force F eff (r) between two electrons in the UEG [37,38], or the effective interaction between two bare nuclei in the interior of dense astrophysical objects [39,40].
In the present work, we present a detailed investigation of the third point based on ab initio PIMC simulations of the UEG without any restrictions on the nodal structure of the thermal density matrix, and without any limitations based on linear-response theory. This gives us a direct insight into the remarkable absence of finite-size errors in QMC calculations of wavenumber resolved properties such as the static structure factor S(q), and opens up new possibilities for future decompositions related to Kohn's celebrated principle of nearsightedness [41]. In addition, we present unassailable numerical proof for an effective electronic attraction between two electrons at certain parameters, without the mediation by an underlying ionic structure. Moreover, we investigate the impact and physical origin of different nonlinear effects [42], and demonstrate the hands-on utility of our new results for the interpretation of scattering experiments, and the description of the dynamic structure factor.
The paper is organized as follows: In Sec. II, we provide the required theoretical background, including a brief introduction to the UEG model and its dimensionless characteristic parameters (II A), our implementation of the PIMC method (II B), the histogram estimation of the effective force and numerical computation of the corresponding effective pair interaction potential (II C), the quantum effective potential model of Kukkonen and Overhauser based on linear-response theory [37] (II D) and the classical potential of mean force model based on distribution function theories (II E). Sec. III is devoted to the in-depth discussion of our new PIMC simulation results, starting with a hands-on discussion of possible finite-size effects (III A), proceeding with the investigation of the dependence of the effective interaction on thermodynamic parameters such as the density and the temperature (III B), and culminating in the comparison of our PIMC data to different theoretical models (III C). Moreover, in Sec. IV, we use our new PIMC results to estimate the effective potential (IV A), and apply it to the interpretation of the previously reported roton feature in the dynamic structure factor of the UEG [43] within the framework of the recently proposed electronic pair alignment model [44] (IV B). We conclude our work with a brief summary of our main findings and a list of possible future investigations in Sec. V. Finally, an exact effective electronic potential connection between quantum linear response theory and classical distribution function theory is established in the appendix.

II. THEORY
A. Uniform electron gas and effective system parameters Throughout this work, we consider an unpolarized electron gas (i.e., N ↑ = N ↓ = N/2, with N being the total number of electrons) in a homogeneous rigid neutralizing background [17,26,45]. The Hamiltonian can then be expressed in the general form (we use Hartree atomic units throughout the main text) where φ E (r a , r b ) is the usual Ewald pair interaction with r a denoting the electron positions as discussed e.g. in Ref. [46]. We note that we do not apply any additional external potential as it has been done in a number of recent works [42,[47][48][49]. From a physical perspective, the UEG is fully characterized by two effective dimensionless parameters [50]: a) the density parameter (also known as Wigner-Seitz radius in the literature) r s = a/a B , with a and a B being the mean distance to the nearest neighbour and the first Bohr radius, respectively, and b) the degeneracy temperature θ = k B T /E F , with E F denoting the Fermi energy [26], T the temperature, and k B the Boltzmann constant. It is common practice to introduce a coupling parameter Γ = V / K . The interaction energy and kinetic energy scale as V ∼ 1/r s and K ∼ 1/r 2 s in the degenerate case, respectively. Hence, it holds Γ ∼ r s , which means that the density parameter plays the role of a quantum coupling parameter in the UEG. The range of metallic densities r s = 1, . . . , 5 can, thus, be characterized as moderately coupled such that weak-coupling expansions like the random phase approximation (RPA) [26] are not strictly applicable. For r s → 0, the UEG converges towards an ideal (i.e., noninteracting) Fermi gas. Conversely, for increasing r s it forms an electron liquid [43,51] and, in the limit r s 1, eventually a Wigner crystal [52,53] . In addition, the degeneracy temperature θ constitutes a straightforward measure for the importance of quantum degeneracy effects, with θ 1 and θ 1 corresponding to the classical limit [54] and fully degenerate case, respectively. In the present work, we are primarily interested in the intermediate regime of r s ∼ θ ∼ 1. These so-called warm dense matter (WDM) conditions are characterized by the complicated interplay of thermal, Coulomb coupling, and quantum degeneracy effects and naturally occur in a number of astrophysical applications such as giant planet interiors [55,56]. Moreover, these extreme states can be realized in experiments at large-scale research facilities using different techniques, see e.g. the topical review by Falk [57].
For completeness, we note that a third relevant parameter is the spin-polarization ξ = (N ↑ −N ↓ )/N , where ξ = 0 and ξ = 1 correspond to the unpolarized (paramagnetic) and fully polarized (ferromagnetic) limits. In the present work, we restrict ourselves to the case of ξ = 0; a detailed investigation of the dependence of the effective force and interaction on the spin-polarization is particularly relevant in the presence of an external magnetic field [58] and constitutes an important project for future works.

B. Path integral Monte Carlo
The rigorous theoretical description of the UEG in the WDM regime, the primary focus of the present work, constitutes a formidable challenge, as a number of complex effects [17,59,60] need to be taken into account. In this context, the most promising method is the quantum Monte Carlo (QMC) technique [61], as it can, in principle, give an exact solution to the full quantum many-body problem of interest without any empirical input (such as the a priori unknown exchange-correlation functional of density functional theory [62]). At finite temperature, the most widely used QMC method is the path integral Monte Carlo (PIMC) approach [4][5][6], which is based on the exact isomorphism between the quantum system and a classical system of interacting ringpolymers [63]. A detailed introduction to PIMC has already been presented elsewhere [4,10]. We, therefore, restrict the presentation here to a concise discussion of the underlying idea.
Let us consider the partition function of N unpolarized electrons in the canonical ensemble (i.e., the inverse temperature β = 1/k B T , number density n = N/Ω and volume Ω are fixed), which can be expressed in coordinate space as Here, R = (r 1 , . . . , r N ) T contains the coordinates of both spin-up and spin-down electrons. The sums run over all possible permutations σ i from the respective permutation group S N i (i ∈ {↑, ↓}), whereπ σ i is the corresponding permutation operator. This ensures the correct fermionic anti-symmetry under the exchange of particle coordinates. The basic idea behind the PIMC method is to perform a Trotter decomposition [64] of the a priori unknown matrix elements of the density operator ρ = e −βĤ . This leads to a product of P density matrices, which have to be evaluated at P -times the original temperature. If this parameter is chosen sufficiently large, one can evaluate the density matrix within a suitable high-temperature approximation. In the present work, we employ the simple primitive factorization with = β/P , which becomes exact in the large-P limit as O P −2 [65]. The convergence with P has been carefully checked (we find P = 200 sufficient at the present conditions), and our results are, therefore, quasi-exact. For completeness, we note that higher-order factorizations ofρ have been studied in the literature [66,67] and become important for PIMC simulations in the lowtemperature regime. The PIMC representation of the partition function can then be expressed in the compact form where the new meta-variable X = (R 0 , . . . , R P −1 ) T includes P sets of coordinates of all N particles, and the integration includes the sum over all permutations. A graphical illustration of Eq. (4) is given in Fig. 1, where we show a configuration of N = 3 particles. Evidently, each particle is now represented by an entire path along Shown is a configuration of N = 3 particles in the τ -x-plane. The horizontal dashed lines correspond to the P time slices with the imaginary-time step = β/P . Each particle is represented by a closed path along the τ -direction. The two particles on the right are involved in a pair-exchange cycle [68], and the corresponding weight is thus negative, W (X) < 0.
the imaginary time τ ∈ [0, β]; this is the origin of the ring-polymers from the celebrated classical isomorphism by Chandler and Wolynes [63]. The basic idea of the PIMC method is to use a variation of the Metropolis algorithm [69] to generate a Markov chain of random configurations {X} that are distributed according to the probability P (X) = W (X)/Z, with W (X) being the corresponding configuration weight. For bosons (such as ultracold atoms like 4 He [4]) or hypothetical distinguish-able quantum particles (often denoted as boltzmannons), this procedure is straightforward. Indeed, using modern Monte Carlo sampling techniques [10,11] quasi-exact results for up to N ∼ 10 4 particles for such system are obtained.
In stark contrast, the situations is incomparably more complicated in the case of fermions (such as the electrons in which we are interested in the present work, but also ultracold atoms like 3 He [18]). Here, the configuration weight W (X) is negative for an odd number of pair exchanges in a given PIMC configuration, c.f. the sign function sgn(σ ↑ , σ ↓ ) in Eq. (2). A corresponding example is shown in Fig. 1, where the two particles on the right are involved in a so-called permutation cycle [68], i.e., a trajectory with more than a single particle in it.
From a practical perspective, the fermionic antisymmetry under the exchange of particle coordinates directly implies that P (X) cannot be interpreted as a probability in the Monte Carlo procedure. To overcome this obstacle, we define a modified configuration space as and randomly generate the paths according to P (X). In particular, Eq. (5) means we carry out a bosonic PIMC simulation; the exact fermionic expectation value of interest is subsequently extracted by evaluating the ratio Here S(X) = W (X)/|W (X)| is the sign of a particular configuration, and the denominator in Eq. (6) is commonly being referred to as the average sign. The sign changes in both the enumerator and the denominator in Eq. (6) lead to a cancellation of positive and negative terms, which is the root cause of the notorious fermion sign problem [19,20,70]. In practice, the latter manifests as an exponential increase in the required compute time with important system parameters such as the system size N or the inverse temperature β [20]. Due to its fundamental nature, a number of strategies have been proposed to tackle the sign problem. Ceperley and others [71][72][73][74] have introduced the fixed-node approximation to the PIMC method, sometimes referred to as restricted PIMC (RPIMC) in the literature. On the one hand, RPIMC formally avoids the exponential bottleneck and, therefore, allows for simulations at parameters that are unfeasible otherwise. Yet, this advantage comes at the cost of an uncontrolled, systematic approximation, effectively breaking the quasi-exact nature of PIMC. Indeed, Schoof et al. [75] have been able to assess the accuracy of RPIMC at WDM conditions based on exact configuration PIMC (CPIMC) simulations. They have found systematic deviations in the exchange-correlation energy between the two methods in excess of 10%. This The black arrows indicate the bare Coulomb forces, and the blue arrow shows the total pair force on the blue reference particle. The green bead depicts a second particle at a distance r, and the green arrow the corresponding effective force F eff (r), i.e., the projection of the total force along the direction of the difference vector r.
surprisingly large bias has subsequently been substantiated in independent studies [17,[76][77][78]. This unsatisfactory situation has sparked a recent surge of activity regarding the development of fermionic QMC methods at finite temperature [76,77,[79][80][81][82][83][84][85][86][87]; see Ref. [78] for a topical discussion of some of these developments. In the present work, we use the standard PIMC method that is heavily afflicted by the sign problem in the WDM regime [20]. As a consequence, our simulations are computationally expensive (we use up to ∼ 10 5 CPUh for a single run, distributed among ∼ 10 4 cores on a modern high-performance computing cluster), but exact within the stated Monte Carlo error bars. All PIMC data of our study are freely available online [88] and can serve as a rigorous benchmark for the development of new methods and approximations.

C. Effective electronic force in an electronic medium
The central goal of the present work is to accurately compute the effective force F eff (r) between two electrons immersed in a UEG, i.e., the average force between two particles at a certain distance r. By definition, this quantity must not include the direct exchange and correlation effects between the two electrons themselves [26]. This is schematically illustrated in Fig. 2, where we show a configuration of N = 6 particles in coordinate space.
The blue bead represents an (arbitrarily chosen) reference particle at a position r a . The black arrows correspond to the individual pair forces on that particle computed from the gradient of the Ewald pair potential φ E (r a , r l ) [46,83] mentioned above. The total force on the reference particle is depicted as the blue arrow in Fig. 2, and is given by To compute the sought-after effective force between a second particle (green bead in Fig. 2) at the position r b , with r ≡ r a − r b and r ≡ |r|, we project the total force F a onto the direction of r, see the green arrow. The relatively short distance between the green and blue particles leads to a comparably strong bare Coulomb force, i.e., the black arrow pointing downwards. The presence of the surrounding electronic medium effectively reduces this force. In other words, the two red particles at the bottom of the figure push the blue reference particle towards the green particle and, thus, reduce the repulsion between them. For bosons, and distinguishable Boltzmann particles, the estimation of this effective force in our PIMC simulations is then straightforward. On each imaginary-time slice, one has to i) compute the full force on each particle F l and ii) project that force along the direction of all N − 1 particle pairs. Note that step ii) is repeated for all N particles on each slice. In practice, this results in a histogram estimator, where one can measure F eff (r) only for certain distances r = r l,k in each configuration.
For fermions, on the other hand, the situation is somewhat more subtle. Specifically, the simple histogram estimator does not take into account fermionic exchange effects, but only implicitly corrects for correlations in the bosonic reference system, cf. Sec. II B above. A more rigorous expression for the effective force is given by where the second term corresponds to the actual measurement of the force for certain configurations, and the first term to the corresponding normalization. It is easy to see that this definition exactly fulfills the notion that F eff (r) must exclude any correlation effects between the two reference particles. We note that Eq. (10) is the most general form and holds for every type of quantum statistics; the histogram estimator discussed above is recovered for bosons, Boltzmann particles, and in the classical limit. Re-calling the evaluation of a fermionic expectation value given in Eq. (6) above, Eq. (10) becomes we have already made use of the fact that the expectation values of the average sign cancel from both terms. Evidently, both the normalization and the actual force measurement in Eq. (11) sensitively depend on fermionic exchange effects via the sign operatorŜ that has to be evaluated in each bosonic reference configuration. Specifically, this consistently fermionic definition of the normalization distinguishes Eq. (10) from the simple histogram, which can be expressed as This has a profound impact on the force between two electrons of the same spin-orientation for r r s , as it is shown in Sec. III C below. Finally, we note that results for F eff (r) can be used to compute a corresponding interaction potential by numerically solving the one-dimensional integral D. Linear response theory and the Kukkonen-Overhauser potential Within linear response theory, the effective potential between any two electrons in the UEG was first derived by Kukkonen and Overhauser (KO) [37], see also the more recent Refs. [26,38]. Their formulation requires the knowledge of accurate dynamic or static local field corrections [89] and provides a convenient reference point for comparison of our PIMC results.
It should be first noted that the general KO treatment includes a contribution that describes the effects of lattice screening on the effective electronic potential; a term that is important for real metals (see the phonon mediated effective electronic attraction that drives a system to a superconducting BCS state), but is non-existent for the UEG [26,37]. For the unpolarized case of interest, the spin-resolved Fourier transformed KO effective potential is given by a two-electron operator in spin space that reads as [26] where χ nn (q, ω) is the density-density response function, χ ss (q, ω) is the longitudinal spin-spin response function, G + (q, ω), G − (q, ω) are the spin symmetric and spin antisymmetric dynamic local field corrections (DLFC),σ 1 , σ 2 are the relevant spin-operators and Φ(q) = 4π/q 2 is the Fourier transformed bare Coulomb interaction.
In the unpolarized case, the latter contribution cancels out when interest lies in the effective potential between any two electrons. Restricting oneself to the static limit, the KO effective potential ultimately reads as [26] with χ nn (q, 0) = χ(q) the static density-density response that describes the linear density response of the full paramagnetic UEG to an external perturbation, i.e., of both spin-up and spin-down electrons, and where G + (q, 0) = G(q) is the static local field correction (SLFC). The static density-density response is conveniently expressed as [93] where χ 0 (q) is the density response of a non-interacting Fermi gas at the same conditions that can be readily computed [26]. Therefore, the effective KO potential is given exclusively as a functional of the SLFC G(q), which contains the full wavenumber-resolved information about static exchange-correlation effects; setting G(q) ≡ 0 thus corresponds to the RPA. In the present work, we use the highly accurate neural-net representation of G(q; r s , θ) by Dornheim et al. [22] that is based on extensive PIMC simulation data. For completeness, we note that the recent analytical representation of G(q; r s , θ) [94] within the effective static approximation [95] would be equally suitable for this purpose. Finally, it is evident that, within linear response theory, the effective force between two electrons in the UEG is given by which can be evaluated both using the SLFC and within the RPA.
E. Classical theory of liquids and the potential of mean force The effective potential between two particles is an important quantity in the classical theory of liquids, where it has been coined as the potential of mean force [96][97][98][99]. The effective potential was formally introduced in a seminal article by Kirkwood [100], although its physical meaning had also been discussed in earlier works of Onsager, Einstein and Smoluchowski [101]. In the canonical ensemble, it is elementary to show that F cl (r) = ∇ r ln [g(r)] which leads to [98][99][100] w (2) where βw (2) (r) is the interaction potential between two particles held apart at a fixed distance r with the contributions of the remaining N − 2 particles ensemble averaged over all canonical configurations. Note that the formal equation and its physical interpretation can be generalized to arbitrary correlation order, i.e., βw (s) (r s ) = ln [g s (r s )], with g (s) (r s ) the s−correlation function and r s = (r 1 , r 2 , ..., r s ) the reduced configuration vector [98]. The potential of mean force βw (2) (r) is invoked in the non-probabilistic interpretation of the Kirkwood superposition approximation which amounts to the assumption that the potential of mean force for a triplet of particles βw (3) (r) is pair-wise additive [96,98], the physical interpretation of the BGY hierarchy (the thermodynamic equilibrium version of the BBGKY hierarchy) which can be viewed as a force balance equation [99] as well as the physical interpretation of the Kirkwood charging equation which can also be viewed as an expression of force balance [98]. Furthermore, the short-distance determination of the potential of mean force from computer simulations is the most difficult step in the indirect extraction of bridge functions [102][103][104]; the key diagrammatic object of the integral equation theory of classical liquids. It is worth noting that Eq. (18) generalizes the Boltzmann relation g(r) = exp [−βφ(r)], that is exact in the dilute classical limit n → 0, to any order in density [105,106].
Moreover, the potential of mean force leads to the definition of a screening potential βH(r) = βφ(r)−βw (2) (r), which is a measure of the residual beyond-pair potential interaction that remains finite even at contact r = 0. The screening potential is a very important quantity for the classical one-component plasma [107], for which its short range behavior can be well-approximated with the aid of the Widom even-power expansion for the cavity distribution function y(r) = g(r) exp [βφ(r)] [108], the Jancovici theoretical result for the second-order term [109] and the consideration of the interaction-site molecule as an extreme case of binary ionic mixture for the zeroorder term [110]. In particular, the leading order term in the enhancement of nuclear reaction rates due to effective interactions in the interior of dense astrophysical objects can be approximated by exp [βH(0)] [39,40], when the semi-classical Wentzel-Kramers-Brillouin (WKB) approximation is applicable for the solution of the relevant two-particle Schrödinger equation [111,112].
Finally, the physical trends of the potential of mean force can be easily deduced from Eq. (18) and the known functional behavior of the pair correlation function. At weak coupling, g(r) exhibits a monotonic behavior that is characterized by a pure exponential asymptotic decay to unity. Therefore, the potential of mean force is monotonic and positive everywhere. At strong coupling, g(r) has a well-defined extended correlation void and then exhibits a non-monotonic profile that is characterized by an exponentially damped oscillatory behavior around unity. As a consequence, the potential of mean force is highly non-monotonic and features alternating regions of effective attraction and repulsion that grow weaker with the distance. Identical to g(r), the transition from monotonic to non-monotonic behavior formally defines a line in the phase diagram known as the Fisher-Widom [113] or the Kirkwood line [114], depending on the mathematical origin of the transition.
As deduced from Eq. (18), the maximum of the first coordination shell corresponds exactly to the global maximum of the attractive effective force. Any particle essentially requires another particle to be quasi-localized therein for energy minimization purposes and the background ensures that another particle is attracted towards that region. This is followed by a local maximum of the repulsive effective force that corresponds exactly to the minimum of the first coordination shell. Essentially, the quasi-localization of a particle at the coordination shell maximum is primary responsible for this repulsive force that prevents particles from occupying the nearby coordination shell minimum. In this manner, the unfolding of short range order generates alternating effective attractions and repulsions.
In the present work, we use highly accurate ab initio PIMC results for the UEG to assess the validity of these classical trends in a strongly coupled quantum system.

III. RESULTS
All PIMC simulation results shown in this work have been obtained on the basis of the extended ensemble scheme that has been introduced in Ref. [115]. A repository of our simulation data is freely available online [88].

A. Finite-size effects
In the preceding sections, we have asserted that our PIMC simulations give us the exact solution to the full quantum many-body problem of interest. While factually true, this statement neglects an important point: PIMC simulations are carried out for a finite system size N . In contrast, practical applications require information about the thermodynamic limit, i.e., the limit of N, Ω → ∞ where the number density n = N/Ω is being kept constant. This discrepancy leads to the presence of so-called finite-size errors, which must be taken into account [21,24,25,46,83,[116][117][118].
In Fig. 3, we thus carefully analyze the dependence of our PIMC simulation data on the number of electrons N for four different values of the density parameter r s at the electronic Fermi temperature, θ = 1. Let us begin our investigation at r s = 2, shown in Fig. 3 a). This is a metallic density that can be probed experimentally for example in aluminum [119,120]. As a reference, we also include the bare Coulomb force F C (r) = 1/r 2 as the dashed blue curves in all panels.
First and foremost, we notice that the effective force F eff (r) decays remarkably fast with r and does not exhibit the long Coulombic tail. This is a rather profound observation, which can hardly be overstated in its importance. In fact, it is well-known that the long-range nature of the Coulomb interaction can easily lead to a number of divergent terms in different theoretical approaches that require special care to be removed or avoided [26,27,121]. Correspondingly, one might have assumed that ab initio simulation methods would require one to perform simulations of mesoscopic system sizes such as N ∼ 10 6 to overcome these difficulties. In practice, this would have precluded the application of state-of-the-art methods like QMC or density functional theory in many cases. Fortunately, this long-range beast is effectively tamed by the screening in charged many-body systems. It is important to note that this does not only apply to effective ion-ion interactions [26,31,32] that are, on average, screened by the electronic medium, but does also occur in a purely electronic medium with a fixed, non-polarizable positive background. Therefore, the electrons are effectively nearsighted and, in practice, do not experience correlations to electrons beyond the remarkably short range of d ∼ 3r s . We emphasize, though, that the above observation does not imply that the Coulomb interaction can be truncated at d ∼ 3r s by introducing a sharp or shifted cut-off. In other words, the long range nature of the Coulomb pair interactions should be fully treated (see the Ewald summation for periodic systems) for the short range nature of the effective electronic interactions to properly emerge.
This empirical finding constitutes the basis for a number of applications in statistical physics, quantum chemistry, and related disciplines. In particular, one can use QMC methods (or less accurate, but computationally cheaper alternatives such as density functional theory) to accurately estimate exchange-correlation properties, that are of short-range nature, but cannot be estimated from any analytical theory. This information concerning the competition between Coulomb correlations, quantum degeneracy, and diffraction effects can then be combined with other methods to model the large scale behaviour of the system. The combination of computational methods with theoretical models like dielectric theories [122][123][124][125][126] then allows for a highly accurate description of the system of interest everywhere [21, 23-25, 117, 118]. Likewise the fast decay of the effective force can be viewed as a theoretical justification for fragmentation methods [127] that partition the electron density into smaller subsystems [128].
Let us now return to the task at hand, which concerns the assessment of finite-size errors in our PIMC data for the effective force F eff (r). Evidently, any N -dependence in the different data sets is so remarkably small, that it can hardly be discerned with the naked eye. In addition, the vertical arrows at the x-axis show the maximum possible particle separation r N max for each N . Yet, F eff (r) has nearly decayed to zero in all depicted cases, as seen particularly well in the inset that focuses on a magnified segment located around r s ≤ r ≤ 2.5r s . On this scale, we do observe small finite-size effects for both N = 4 and N = 6, with F eff (r) decaying too fast especially for the former case.
In panel b), we show analogous results for a lower density, r s = 4, which is close to sodium [129]. Overall, the behaviour of the PIMC data is very similar compared to the previous case of r s = 2. The main difference lies in the fact that finite-size effects for N = 4 and N = 6 are even smaller, with the N = 4 data set exhibiting the correct asymptotic decay. This is a general trend, which is consistent to previous investigations of finite-size errors in the warm dense UEG [17,24,25], namely: as the density increases (i.e., decreasing the density parameter r s ), larger numbers of QMC simulated electrons are necessary to achieve convergence of wave-number resolved quantities such as the static structure factor S(q). In the present work, we observe that this is a direct consequence of the fact that many-body correlations as they are embodied in F eff (r) have a characteristic scale that exceeds the length of the simulation cell L = Ω 1/3 in these cases.
We proceed with r s = 10 [panel c)], which corresponds to a moderately coupled system at the margin of the electron liquid regime [51]. Such conditions are very interesting for many reasons. Quantum liquids allow one to gain new insights into the interplay of quantum effects with correlations and incipient localization effects. In addition, the UEG is known to exhibit a wealth of intriguing phenomena in this regime, such as the roton feature in the dynamic structure factor S(q, ω) [43,44,130]; see Sec. IV A for a detailed discussion of this point. First and foremost, we note that finite-size errors are almost nonexistent over the entire distance range. In addition, we find an attractive minimum in the effective force around r 2r s . While being comparably small in magnitude at these conditions, this effect is significant and clearly not an artifact due to the finite number of electrons in the simulation cell; indeed, only N = 4 cannot correctly resolve this feature as its location exceeds the maximum possible inter-particle separation in this case. The physical origin of the attraction and its potential consequences for the physics of the UEG are discussed in Sec. III B.
Finally, we also show results for even stronger coupling, r s = 20 [panel d)]. In this case, the attractive minimum in F eff (r) becomes even more pronounced, and is followed by a local maximum around r = 3r s . In fact, a very shallow yet significant second minimum can even be discerned around r = 4r s . These features are a direct consequence of the incipient mesoscopic order in the system upon entering the electron liquid regime [26,51], see also the discussion in Sec. II E.

B. Dependence on density and temperature
In what follows, we shall more closely examine the main physical trends in F eff (r). To this end, we show PIMC simulation results for different values of r s at the electronic Fermi temperature θ = 1 in Fig. 4 a). We note that we have re-scaled the forces by a factor of r 2 s to make them comparable between different densities. Evidently, all depicted data sets exhibit a qualitatively very similar behaviour. At large distances, they quickly decay towards zero and do not exhibit the long-range Coulomb tail, as discussed at length in the previous section. Moreover, the effective force is dominated by Coulomb repulsion in the limit of small distances on the depicted scale. The main feature in Fig. 4 is the density-dependence of the effective attraction between the electrons, which we find for r s 5 in our PIMC results. This is a signature of short range order, as discussed in Sec. II E. Naturally, electronic attraction suggests a kind of pairing mechanism such as the formation of Cooper pairs in the BCS theory for superconductivity [131]. Yet, in that case, the effective attraction is mediated by phonons, i.e., by the underlying ionic structure. In contrast, the rigid, uniform, and non-polarizable background of the UEG, by definition, cannot be the origin of this effect. Instead, it is purely due to the electronic medium, and can be intuitively understood by considering the sketch in Fig. 2 shown above. Let us again examine the blue and green beads; clearly, the former is effectively pushed towards the latter by the two red beads at the bottom. While in the depicted example the Coulomb repulsion between the green and blue beads still exceeds this push, this need not be the case. In particular, the UEG undergoes an incipient localization when the density is decreased. As a direct consequence of the role of r s as the quantum coupling parameter, the electrons get 1) more clearly separated and 2) more ordered. In other words, they are, on average, pushed by the medium towards a spatial structure that minimizes their free energy landscape.
To further illustrate this effect, we show PIMC results for the pair correlation function g(r) in Fig. 4 b). Evidently, g(r) exhibits a (small) peak exactly for those r s -values where we find the attractive feature in F eff (r). For completeness, we note that the peak location is sim-ilar, but is not identical to the location of the latter, in contrast to the classical limit, where the F eff (r) and g(r) peak locations coincide, see Eq. (18) and Fig. 8.
Let us briefly touch upon the possible physical consequences of the effective electronic attraction in the UEG. Firstly, it constitutes a direct measure for the degree of UEG collective behaviour, which gives rise to various interesting phenomena such as the spectrum of density fluctuations [44,132]. In addition, it has been speculated that the attraction might constitute a possible pairing mechanism that would eventually lead to superconductivity in the UEG at r s 40 and low temperatures [92]. Similar possibilities have been suggested for semiconductors [133] and metallic hydrogen [90] in the literature. A detailed study of these conditions is beyond the scope of the present work and deserves dedicated exploration in the future.
A further interesting variable is the dependence of the attraction on the reduced temperature θ, which we investigate in Fig. 5a) for r s = 10. Remarkably, the overall dependence of the effective electronic force in the UEG on the temperature is quite small. The most pronounced differences that we find concern the attractive feature around r ∼ 2r s . While the two curves for θ = 0.75 and θ = 1 are in very close agreement, the minimum becomes very shallow for θ = 2 and completely vanishes for θ = 4. This is again directly correlated to the presence of a peak in g(r), which we show in Fig. 5b) at the same conditions.

C. Comparison to theoretical models
A further interesting topic of investigation is the capability of theoretical models to capture the microscopic effective interaction between two electrons in the UEG. This is explored in Fig. 6. Specifically, the left column shows F eff (r) from different theories, and the top row corresponds to r s = 10 and θ = 1. As usual, the dashed blue line shows the bare Coulomb force that has been included as a reference, and the black symbols depict our quasiexact PIMC data [cf. Eq. (11)] with their corresponding error bars. In addition, the dotted green line shows the linear-response prediction within the RPA, which has been obtained from Eqs. (15)(16)(17) by setting G(q) ≡ 0. Evidently, the RPA converges towards Coulomb at small distances. This can be seen particularly well in Fig. 6  b), where the ratio of the respective data sets to the bare Coulomb force has been plotted. Interestingly, we observe that the exact PIMC results start to systematically deviate from the Coulomb force in the limit of very small distances. For completeness, it is worth pointing out that the increasing error bars in the PIMC data towards contact are a direct consequence of the histogram estimation of F eff (r), since pair-particle encounters become exceedingly infrequent within our simulations as the distance decreases. In principle, this bottleneck can be overcome in future works via the implementation of umbrella sampling techniques that were originally introduced in classical Monte Carlo methods to improve the sampling of ultra-rare configurations [134,135]. The RPA substantially overestimates the drop in F eff (r) due to the medium for r ∼ r s , and it exhibits Friedel oscillations [136] for r 2r s that slowly decay with increasing r. We note that these oscillations are not a real feature of the UEG at the present conditions, but are a direct consequence of the insufficient treatment of electronic exchange-correlation effects within the RPA.
In contrast, the solid red curve has been obtained by evaluating the KO potential [26,37] with input from the neural-net representation of the static LFC G(q; r s , θ) from Ref. [22], which is based on accurate PIMC data. We note that the recently introduced analytical representation of the SLFC [94] within the effective static approximation [95] would basically lead to indistinguishable results. In practice, we find that the red curve constitutes a substantial improvement over the RPA, and exhibits the correct qualitative decay towards zero for r r s . At the same time, qualitative differences arise which deserve a further exploration. In particular, the KO expression overestimates the magnitude of the effective electronic attraction at r ∼ 2r s , which can be roughly understood as follows. The minimum in F eff (r) is a direct consequence of many-body effects, and the increasing impact of Coulomb correlations. Specifically, two particles are effectively pushed together by the surrounding medium. Yet, it is well known that the usual linear response does only explicitly include two-body correlations (although all orders are implicitly present in the SLFC), since three-body correlations are connected to the quadratic density response function, and so on [137]. Therefore, we conclude that the linear-response description within the framework of KO cannot fully capture the effective attraction, which would require the incorporation of nonlinear effects [42,47,138], see also the appendix for the analysis of the classical limit. In addition, the KO force monotonically converges towards the bare Coulomb force in the limit of r → 0 [see also Fig. 6 b)], and does not exhibit the relative drop of the exact PIMC results.
To acquire additional insights into the physical mechanism behind the effective force, we carried out PIMC simulations using Boltzmann statistics. More specifically, these calculations include the full information about quantum diffraction and related effects, but the particles are assumed to be distinguishable; the anti-symmetry of the true fermionic density matrix under the exchange of particle coordinates is completely omitted. The results are shown as the yellow curve in Fig. 6, and are in excellent agreement with the true PIMC data over the entire r-range. We stress that this is a highly counter-intuitive finding, which calls for an explanation. In particular, it is well understood that fermionic exchange effects should substantially influence the behaviour of two electrons at sufficiently small distances r at the present conditions. Indeed, if that was not be the case, there would be no need for fermionic PIMC simulations, and the notorious fermion sign problem would have been effectively circumnavigated.
To understand this remarkable behaviour of F eff (r), we need to recall its definition, and the corresponding estimation within the PIMC method discussed in Sec. II C above. Specifically, Eqs. (10) and (11) imply that the effective force does not include exchange and correlation effects between the two electrons themselves; this has been ensured by dividing by the expectation value of the corresponding normalization. In other words, our estimator corrects for fermionic exchange effects, which leads to the observed equality between fermions and Boltzmann particles. Aiming to further illustrate this point, we also implemented the naive histogram estimator discussed in the beginning of Sec. II C. In this way, we only correct for correlation effects in the fermionic simulation, but do not take into account the impact of the fermionic antisymmetry under the exchange of particle coordinates.
The results are included as the light green diamonds in Fig. 6. Evidently, they are in excellent agreement to the exact PIMC and Boltzmann PIMC results for r r s , but strongly deviate for smaller r. In other words, fermionic exchange effects between two electrons only have an impact at very small distances, but do not significantly influence other quantum effects such as diffraction. In practice, the Pauli blocking prevents two identical fermions from occupying the same position, and can be interpreted as an additional effective force that separates two particles (i.e., the "degeneracy pressure"). Yet, this effective repulsion does not directly show up in the histrogram estimator for F eff (r). Indeed, we use the PIMC method to compute the expectation value of Eq. (8), where we only sum over the Ewald forces. In contrast, the Pauli repulsion manifests indirectly due to the cancellation of positive and negative terms during the PIMC simulation [20,70]. In practice, this leads to the following: holding together two identical fermions to short distances r < r s , which is the only case when we sample F eff (r < r s ) in our PIMC simulations, requires the medium to counteract the effective Pauli repulsion. Yet, the latter does not directly manifest in F eff (r); the counter-acting push bringing the electrons together, on the other hand, leads to the observed reduction in the light green curve. On the other hand, the consistent definition given by Eq. (11) completely removes this effect. Let us next discuss the bottom row of Fig. 6, where we show the same information for the metallic density of r s = 4. In this case, RPA is considerably more accurate compared to r s = 10 as electronic exchange-correlation effects are less important. In addition, we find smaller differences between the full PIMC data (black curve) and the KO results over the entire r-range. This, too, follows from the reduced importance of correlations, as nonlinear effects are directly related to higher-order correlation functions known from many-body theory [137]. Finally, we note that the discrepancy between the exact PIMC results and the histogram estimator for small distances is actually increased compared to r s = 10, as fermionic exchange effects are more pronounced at the higher density [17].
We shall continue this investigation by exploring the spin-resolved effective force shown in Fig. 7. As usual, the black, yellow, and dashed blue curves show the full fermionic PIMC results [Eq. (11)], the Boltzmann PIMC results, and the bare Coulomb force, respectively. In addition, the light green, grey and purple symbols show the effective force computed using the naive histogram estimator between any two electrons, between two electrons of the same spin-orientation (e.g. up-up) and between two electrons of opposite spin (up-down). Evidently, the latter closely follows the exact PIMC and Boltzmann PIMC results, as it is not directly influenced by exchange effects. The up-up force, on the other hand, is strongly reduced compared to the other data sets due to the counter-acting of the effective Pauli repulsion discussed earlier. For completeness, we note that no significant spin-dependence can be resolved in our spin-resolved implementation of the true fermionic expectation value Eq. (11), which is not pictured in Fig. 7.
From a practical perspective, the observed trends open up intriguing new possibilities for the description of Fermi systems, in general, and many-electron systems, in particular. As we have mentioned above, the main bottleneck of QMC simulations is the exponential increase in the compute time with the system-size N . The shortrange nature of the effective force between two electrons leads to a great reduction of this problem, as many-body correlation effects are basically limited to the range of r 3r s in the regime of moderate coupling, r s ∼ 1 − 10. While this makes QMC studies feasible in many cases, a considerable part of the interesting density-temperature range remains out of reach [78,85]. Our present findings clearly indicate, that the impact of fermionic exchangeeffects, which is the most challenging part for all QMC methods, is limited to even shorter length scales, and becomes negligible for r r s . In a certain sense, this finding constitutes an empirical confirmation of Kohn's principle of nearsightedness [41] that has been introduced in the context of density functional theory [139]. In particular, it opens up the enticing possibility for a further decomposition: specifically, we propose to use exact fermionic QMC simulations of very small systems to accurately capture the full interplay of Coulomb correlations with fermionic exchange effects for r r s ; the still formidable effects of quantum diffraction are fully captured by the Boltzmann PIMC simulations, which do not suffer from the exponential bottleneck. Indeed, simulations with N ∼ 10 2 − 10 3 distinguishable electrons are feasible; these information can finally be combined with dielectric theories such as the RPA to get an accurate description of the system over all length scales. The practical implementation of this idea will be explored in detail in a future publication.

IV. APPLICATIONS
A. Effective pair potential Let us next utilize our new PIMC results for the effective force F eff (r) to obtain the effective interaction potential between two electrons in the UEG medium. To this end, we numerically solve the simple one-dimensional integral Eq. (13) and compare the results to the RPA potential and the KO potential with the PIMC-based SLFC in Fig. 8. The left panel corresponds to r s = 10 and phenomenologically resembles the respective effective forces shown in Fig. 6 a). Interestingly, the Friedel oscillations in the RPA are less pronounced, whereas the overall systematic errors of RPA and the accurate linearresponse prediction of the KO potential compared to the true PIMC results are similar to our results for F eff (r). Similarly, the KO potential is quite accurate over the entire r-range, whereas the RPA exhibits more substantial deviations; the same trends can be found in Fig. 8b) for the metallic density of r s = 4. Finally, the solid yellow curves show the classical potential of mean force, which we estimate by inserting our PIMC results for g(r) into Eq. (18). Evidently, the classical expression breaks down for small distances r, where quantum effects such as diffraction and exchange predominantly shape the physical behaviour of the system. Remarkably, though, φ cl (r) nicely captures the attractive minimum around r = 2r s for r s = 10 and is substantially more accurate than the KO result.

B. Spectrum of density fluctuations
Very recently, it has been shown that microscopic information about the effective electronic interaction gives straightforward access to the spectrum of density fluctuations ω(q) [44]. The latter is typically estimated in terms of the maximum of the dynamic structure factor S(q, ω), which can be expressed in the exact spectral representation as [26] S(q, ω) = m,l P m n ml (q) 2 δ(ω − ω lm ) .
We note that l and m denote the eigenstates of the Nbody Hamiltonian Eq. (1), ω lm = (E l − E m )/ is their energy difference, and n ml is the corresponding transition element induced by the singe-particle density operator n(q). In other words, Eq. (19) constitutes the sum over all possible transitions between the (time-independent) eigenstates of the system, with P m being the probability of the initial state. Naturally, a direct evaluation of Eq. (19) constitutes a most formidable task and would, in principle, require the complete diagonalization of the many-body Hamiltonian. This is both unfeasible and not necessary; highly accurate results for S(q, ω) of the warm dense UEG have become available based on the analytic continuation of exact PIMC data for the imaginary-time density-density correlation function F (q, τ ) [43,130,140,141]. These results are based on the fully frequency-dependent DLFC G(q, ω), which allows one to resolve the interesting, nontrivial structure of S(q, ω) at intermediate wavenumbers, q ∼ 2q F , for moderate to strong coupling. For r s 10, this leads to the emergence of a roton feature in ω(q), that is shown in Fig. 9 a). Specifically, the dash-dotted blue curve shows the exact PIMC results from Ref. [43] and the dotted green curve the estimation of ω(q) within the RPA. Clearly, the latter does not capture the pronounced down-shift in the dispersion curve around q ∼ 2q F , i.e., when the wave length λ = 2π/q becomes comparable to the average interparticle distance. Based on the dramatic short-range nature of φ eff (r), Dornheim et al. [44] have proposed the electronic pair alignment model. Within this framework, the roton feature is caused by the excitation of electron pairs, which, in the presence of the electronic medium, interact via Eq. (13). In particular, the probability to find two particles at a distance r prior to the excitation, P (r) = ng(r), plays the role of P m in Eq. (19). The difference between the initial and final state is assumed to be restricted to this particle pair, which, after the excitation, are separated by a distance λ along the direction of q. In the regime of the roton feature, such pair excitations will effectively separate the particles, thereby reducing their interaction energy by the amount [44] Let us briefly re-iterate the physical meaning of Eq. (20). The probability of finding a particle pair of distance r in the unperturbed system is P (r) = ng(r). Due to the excitation, the particles will be separated by the wavelength λ after the excitation along the direction of the perturbation; the other coordinates remaining unchanged. The corresponding final separation is denoted as r q . Clearly, this change in particle distance leads to a change in the effective potential energy ∆φ eff (r, r q ) = φ eff (r) − φ eff (r q ).
The average interaction energy change ∆W (q) of a pair excitation of wavenumber q is thus the average over all such transitions.
To estimate the red-shift in ω(q) compared to RPA, we have to estimate the corresponding exchange-correlation contribution where W RPA (q) is obtained by inserting g RPA (r) and φ RPA (r) into Eq. (20). Since the impact of exchangecorrelation effects onto the kinetic energy is negligible at these conditions, we immediately get the ansatz Here the screening function [93] α(q) = χ(q)/χ 0 (q) takes into account the fact that such pair excitations will become impossible for λ r s due to the perfect screening in the UEG [142].
In lieu of the true PIMC results for φ eff (r), Dornheim et al. [44] have used the linear-response expression due to Kukkonen and Overhauser, Eq. (15), and the results are shown as the solid red lines for r s = 10 [ Fig. 9 a)] and r s = 4 [ Fig. 9 b)]. It is evident that these curves capture the main part of the exchange-correlation induced down-shift compared to RPA in both cases. Interestingly, they closely follow the static approximation [43] (dashed black), which has been obtained by setting the DLFC to its exact static limit, i.e., G(q, ω) ≡ G(q, 0). This is expected, as both the static approximation and Eq. (20) constitute a measure for an average energy shift; see the original Ref. [44] for a more detailed discussion of this point. Within the context of the present work, the most important result is illustrated as the purple curves in Fig. 9, which has been obtained by evaluating Eq. (20) using as input our new PIMC results for the true effective potential, Eq. (13). Overall, the results are in very good agreement to the red curves, which further substantiates the validity of the pair alignment mechanism. In addition, we note that the true φ eff (r) leads to a small improvement towards the exact (dash-dotted blue) PIMC curves for both densities; for r s = 10, we even get a shallow roton feature in the purple curve at the correct position that is not present in the results that have been obtained with the linear-response KO potential.

V. SUMMARY AND DISCUSSION
In this work, we have presented extensive new ab initio PIMC results for the effective force F eff (r) and the effective pair potential φ eff (r) between two electrons in the presence of the UEG. First and foremost, we have consistently observed a rapid decay of F eff (r) towards large r, which is in stark contrast to the long-range tail of the bare Coulomb interaction. This has profound implications for computational quantum many-body theory, and explains the observed absence or the small magnitude of finite-size effects in wavenumber resolved properties, such as S(q), that was observed in earlier studies of the UEG [17,21,22,24,25,116,118]. Indeed, finite-size effects are negligible in both F eff (r) and φ eff (r) for as few as N = 14 electrons at WDM conditions. From a physical perspective, our study has provided the numerical proof for an effective attraction between two electrons around r ∼ 2r s for r s 5 and θ 2. We re-iterate our earlier point that this does not arise due to the mediation through phonons or other ionic effects, and that the rigid, non polarizeable background in the UEG model cannot be the origin of this effect. Instead, it is a direct consequence of the incipient localization of the electrons upon entering the electron liquid regime.
The comparison of our new PIMC simulation data to analytical theories has revealed that the widely used random phase approximation does not provide an adequate description of the microscopic interaction between a pair of UEG electrons at metallic densities. In stark contrast, the linear-response expression by Kukkonen and Overhauser [37] allows us to consistently incorporate electronic exchange-correlation effects based on previous PIMC data for the static local field correction. This leads to a dramatic improvement over the RPA, and accurately captures the main trends of the exact PIMC results. At the same time, the KO results overestimate the depth of the attractive minimum in both F eff (r) and φ eff (r). This is a direct consequence of the correlational origin of the effective attraction, which, therefore, requires a more advanced description that incorporates nonlinear terms [137]. In addition, we have found substantial deviations between the true fermionic PIMC estimator for the effective force and the naive histogram estimator for r < r s , where the former approaches the bare Coulomb repulsion, whereas the latter is comparably reduced in magnitude. These deviations have been explained by the fermionic degeneracy pressure, which is consistently removed from the true force, but biases the simple histogram implementation.
Finally, we have used our new PIMC results for the effective electronic pair interaction φ eff (r) to estimate the spectrum of density fluctuation ω(q) within the recently introduced electronic pair alignment model. Our results further validate previous data that were based on the linear-response KO potential, and constitute a small improvement over the latter. We, thus, safely conclude that the studied effective interaction gives one direct microscopic insight into the origin of the roton feature in the warm dense UEG.
We are convinced that our investigation will open up many avenues for important and impactful future research, and briefly sketch a number of promising possibilities: i) our PIMC results give one direct insight into the microscopic impact of nonlinear effects onto the effective interaction between a pair of electrons. This can be used as the basis for the further improvement of dielectric theories [123,126,143,144] and guide the development of explicitly nonlinear potentials [145,146]; ii) the possible role of the attractive minimum in F eff (r) as a pairing mechanism in the UEG, and its connection to the potential emergence of superconductivity in the electron liquid regime [92] deserve closer exploration; iii) our analysis has revealed the remarkably fast decay of fermionic exchange effects for r ∼ r s , which opens up the intriguing possibility for a further decomposition: use of fermionic PIMC for very small N to study short range exchange-correlation effects at r r s , use of Boltzmann PIMC to study diffraction at r s r 3r s , and use of dielectric theories for r r s to get a highly accurate description across the full domain of length scales. This might substantially increase the scope of feasible density-temperature combinations [85], and give us for the first time access to the regime of metallic densities at θ ∼ 0.1, . . . , 0.5; iv) it is straightforward to estimate both F eff (r) and φ eff (r) for a number of other systems that can be simulated with the PIMC approach. This can give new insights into important phenomena like superfluidity [4,14,147] and supersolid behaviour [148][149][150][151]. In addition, we suggest to use the pair alignment model [44] to explain the spectrum of density fluctuations in other systems, with the original roton feature in ultracold 4 He [14,147,152] and 3 He [18,[153][154][155] being a promising candidate; v) finally, we will extend the present study beyond the UEG model, starting with hydrogen in the WDM regime [156,157]. This will allow us to investigate the effective interaction between all individual components, and provide us with unprecedented insights into the behaviour of WDM at the nanoscale.

APPENDIX: CLASSICAL RELATION BETWEEN THE KUKKONEN-OVERHAUSER POTENTIAL AND THE POTENTIAL OF MEAN FORCE
In what follows, we shall derive an exact classical relation between the Kukkonen-Overhauser potential and the potential of mean force. This classical connection is not restricted to Coulomb interactions but it is valid for arbitrary pair interaction potentials (even anisotropic). It can be exploited in order to systematically improve the effective interaction predictions of linear response theory without significant effort.
The general static form of the effective KO potential in the unpolarized case remains intact in the classical limit. Irrespective of the pair interaction potential, it reads as βΦ KO eff (q) = βΦ(q) + 1 n {nβΦ(q)[1 − G(q)]} 2 χ(q) nβ .
The classical version of the fluctuation-dissipation theorem, the zeroth frequency moment rule and the Kramers-Kronig relation lead to the well-known expression S(q) = −χ(q)/(nβ). Together with the static limit of the ideal Vlasov density response χ 0 (q) = −nβ and static limit of the general relation for the density response, Eq.16, the above expression yields nβΦ(q)[1 − G(q)] = [1/S(q)] − 1. Thus, In order to further simplify the above, the Fourier space connection between the static structure factor and total correlation function S(q) = 1 + nH(q) is employed and the Fourier transformed exact Ornstein-Zernike equation H(q) = C(q) + nH(q)C(q) is utilized, where C(q) is the Fourier transformed direct correlation function c(r). The resulting expression can be trivially inverted back to real space and reads as βφ KO eff (r) = βφ(r) + c(r) − h(r) .
Recall that the exact non-linear closure equation of the integral equation theory of classical liquids is given by g(r) = exp [−βφ(r) − c(r) + h(r) + b(r)] [105,106], with b(r) the bridge function. Application of the natural logarithm, disposal of the common βφ(r) + c(r) − h(r) expression and substitution for the potential of mean force, ultimately yield βφ KO eff (r) = βφ cl (r) + b(r) .
Note that, thus far, only exact classical expressions have been employed, which implies that βφ cl (r) is the exact effective interaction. The above expression demonstrates that the effective KO potential is only exact within the hypernetted chain approximation (HNC), which assumes that the bridge function is zero, b(r) = 0. Hence, in the classical limit, it can be deduced that the use of linear response theory is exactly equivalent to the approximation that the contribution of highly connected irreducible diagrams (which constitute the bridge function) to the pair correlation function is negligible. Treatment of the general quantum case with diagrammatic perturbation theory has revealed that the KO derivation implicitly contains two additional approximations [26,133]: the substitution of the internal electron-hole propagator by non-interacting electron-hole propagators and the replacement of the irreducible electron-hole interaction by some kind of an average. In the dilute strongly coupled uniform electron fluid regime, r s 20, Eq. (A2) can be employed to improve the effective potential predictions of the KO expression. At lower densities, fermionic exchange effects are far less pronounced and the same should apply for the consequences of the implicit approximations of the KO derivation. Simultaneously, irreducible diagram effects become significant but they can be treated classically, as deduced from the success of the integral equation theory based dielectric scheme at strong coupling [123,143]. Given the availability of a closed-form expression for the bridge function of the classical one component plasma at arbitrary coupling [107,143], βφ eff (r) = βφ KO eff (r) − b(r) should constitute an improvement over the original KO expression. This approximation will be explored in a future work that will present PIMC results beyond the WDM regime.