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 *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. 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 behavior were 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–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–16}

Unfortunately, QMC simulations of quantum degenerate Fermi systems, such as the UEG,^{17} or ultracold Fermi atoms, such as ^{3}He,^{18} are afflicted with the notorious fermion sign problem,^{19,20} which leads to an exponential increase in the computation time with increasing system size or decreasing temperature. Therefore, QMC simulations of many-electron systems are usually restricted to rather small systems [typically $N\u223cO101\u2212102$] even using state-of-the-art algorithms on modern high-performance computing clusters. At the same time, we note that these simulations are often remarkably accurate with as few as *N* = 14 electrons^{17,21–25} and, in combination with dielectric methods (such as 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 long-range 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 milestone in the description of quantum many-body systems. In general, we distinguish three distinct situations:^{26,28}

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–34}The effective interaction between a test change and the surrounding medium, such as an ion embedded in a uniform electron gas.

^{35,36}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–39}or the effective interaction between two bare nuclei in the interior of dense astrophysical objects.^{40,41}

In this work, we present a detailed investigation of the third situation based on *ab initio* PIMC simulations of the UEG without any restriction 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.^{42} In addition, we present unassailable numerical proof for an *effective electronic attraction* between two electrons for certain parameters, without the mediation by an underlying ionic structure. Moreover, we investigate the impact and physical origin of different nonlinear effects^{43} and demonstrate the hands-on utility of our new results for the interpretation of scattering experiments and the description of the dynamic structure factor.

This 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 (Sec. II A), our implementation of the PIMC method (Sec. II B), the histogram estimation of the effective force and numerical computation of the corresponding effective pair interaction potential (Sec. II C), the quantum effective potential model of Kukkonen and Overhauser (KO) based on linear-response theory^{37} (Sec. II D), and the classical potential of mean force model based on distribution function theories (Sec. II E). Section 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 (Sec. III A), proceeding with the investigation of the dependence of the effective interaction on thermodynamic parameters, such as the density and the temperature (Sec. III B), and culminating in the comparison of our PIMC data to different theoretical models (Sec. III C). Moreover, in Sec. IV, we use our new PIMC results to estimate the effective potential (Sec. IV A) and apply it to the interpretation of the previously reported *roton feature* in the dynamic structure factor of the UEG^{44} within the framework of the recently proposed *electronic pair alignment* model^{45} (Sec. 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,46} 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. 47. We note that we do not apply any additional external potential as it has been done in a number of recent works.^{43,48–50}

From a physical perspective, the UEG is fully characterized by two effective dimensionless parameters:^{51} (a) the density parameter (also known as Wigner–Seitz radius in the literature) $rs=a\u0304/aB$, with $a\u0304$ and *a*_{B} being the mean distance to the nearest neighbor 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\u223c1/rs2$ in the degenerate case, respectively. Hence, Γ ∼ *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, such as the random phase approximation (RPA),^{26} are not strictly applicable. For *r*_{s} → 0, the UEG converges toward an ideal (i.e., noninteracting) Fermi gas. Conversely, for increasing *r*_{s}, it forms an electron liquid^{44,52} and, in the limit *r*_{s} ≫ 1, eventually a Wigner crystal.^{53,54}

In addition, the degeneracy temperature *θ* constitutes a straightforward measure of the importance of quantum degeneracy effects, with *θ* ≫ 1 and *θ* ≪ 1 corresponding to the classical limit^{55} 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 and Coulomb coupling and quantum degeneracy effects and naturally occur in a number of astrophysical applications such as giant planet interiors.^{56,57} 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.^{58}

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^{59} and constitutes an important project for future works.

### B. Path integral Monte Carlo

A 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,60,61} need to be taken into account. In this context, the most promising method is the quantum Monte Carlo (QMC) technique,^{62} 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^{63}). At finite temperature, the most widely used QMC method is the path integral Monte Carlo (PIMC) approach,^{4–6} which is based on the exact isomorphism between the quantum system and a classical system of interacting ring-polymers.^{64} 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 follows:

Here, $R=(r1,\u2026,rN)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 $SNi$ (*i* ∈ {*↑*, *↓*}), where $\pi \u0302\sigma i$ is the corresponding permutation operator. This ensures the correct fermionic antisymmetry under the exchange of particle coordinates. The basic idea behind the PIMC method is to perform a Trotter decomposition^{65} of the *a priori* unknown matrix elements of the density operator $\rho \u0302=e\u2212\beta H\u0302$. 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 as follows:

with *ϵ* = *β*/*P*, which becomes exact in the large-*P* limit as $OP\u22122$.^{66} 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 $\rho \u0302$ have been studied in the literature^{67,68} and become important for PIMC simulations in the low-temperature regime.

The PIMC representation of the partition function can, then, be expressed in the following compact form:

where the new meta-variable $X=(R0,\u2026,RP\u22121)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 the imaginary time *τ* ∈ [0, *β*]; this is the origin of the ring-polymers from the celebrated classical isomorphism by Chandler and Wolynes.^{64} The basic idea of the PIMC method is to use a variation of the Metropolis algorithm^{70} 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, such as ^{4}He^{4}) or hypothetical distinguishable 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 a 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 and also ultracold atoms such as ^{3}He^{18}). Here, the configuration weight *W*(**X**) is *negative* for an odd number of pair exchanges in a given PIMC configuration, cf. 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*,^{69} 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 follows:

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 referred to as the *average sign*. The sign changes in both the numerator 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,71} In practice, the latter manifests as an exponential increase in the required computation 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 co-workers^{72–75} 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.*^{76} 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 surprisingly large bias has subsequently been substantiated in independent studies.^{17,77–79}

This unsatisfactory situation has sparked a recent surge of activity regarding the development of fermionic QMC methods at finite temperature;^{77,78,80–88} see Ref. 79 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 $\u223c105$CPUh for a single run, distributed among $\u223c104$ 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^{89} 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 the 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})^{47,84} 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 downward. 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 toward 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. 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 effect 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, for Boltzmann particles, and in the classical limit. Recalling the evaluation of a fermionic expectation value given in Eq. (6) above, Eq. (10) becomes

and 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 $S\u0302$ 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 follows:

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.

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 expression for 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 and 39. Their formulation requires the knowledge of accurate dynamic or static local field corrections (SLFCs)^{90} and provides a convenient reference point for comparison of our PIMC results.

Aside from the uniform electron liquid, applications also concern the theory of superconductivity.^{91–93}

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 nonexistent 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), respectively, $\sigma \u03021$ and $\sigma \u03022$ 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^{94}

where *χ*_{0}(*q*) is the density response of a noninteracting 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}, *θ*)^{95} within the effective static approximation^{96} 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.^{97–100} The effective potential was formally introduced in a seminal article by Kirkwood,^{101} although its physical meaning had also been discussed in earlier works of Onsager, Einstein, and Smoluchowski.^{102} In the canonical ensemble, it is elementary to show that *F*_{cl}(*r*) = ∇_{r} ln[*g*(*r*)], which leads to^{99–101}

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., $\beta w(s)(rs)=lngs(rs)$, with *g*^{(s)}(*r*^{s}) the *s*-correlation function and *r*^{s} = (*r*_{1}, *r*_{2}, …, *r*_{s}) the reduced configuration vector.^{99}

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,^{97,99} the physical interpretation of the Born-Green-Yvon (BGY) hierarchy [the thermodynamic equilibrium version of the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy]—which can be viewed as a force balance equation^{100} as well as the physical interpretation of the Kirkwood charging equation, which can also be viewed as an expression of force balance.^{99} 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;^{103–105} 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\u2212\beta \varphi (r)$, which is exact in the dilute classical limit *n* → 0, to any order in density.^{106,107}

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,^{108} 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*)],^{109} the Jancovici theoretical result for the second-order term,^{110} and the consideration of the interaction-site molecule as an extreme case of binary ionic mixture for the zero-order term.^{111} 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)],^{40,41} when the semiclassical Wentzel–Kramers–Brillouin (WKB) approximation is applicable for the solution of the relevant two-particle Schrödinger equation.^{112,113}

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 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^{114} or the Kirkwood line,^{115} 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 toward 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

### 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 kept constant. This discrepancy leads to the presence of so-called *finite-size errors*, which must be taken into account.^{21,24,25,47,84,117–119}

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.^{120,121} 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,122} 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, such as 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, nonpolarizable positive background. Therefore, the electrons are effectively *near-sighted* and, in practice, do not experience correlations to electrons beyond the remarkably short range of *d* ∼ 3*r*_{s}. We emphasize, though, that the above observation does not imply that the Coulomb interaction can be truncated at *d* ∼ 3*r*_{s} by introducing a sharp or shifted cutoff. 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 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 behavior of the system. The combination of computational methods with theoretical models, such as dielectric theories,^{123–127} then, allows for a highly accurate description of the system of interest everywhere.^{21,23–25,118,119} Likewise, the fast decay of the effective force can be viewed as a theoretical justification for fragmentation methods^{128} that partition the electron density into smaller subsystems.^{129}

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 datasets 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, $rmaxN$, 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.5*r*_{s}. On this scale, we do observe small finite-size effects for both *N* = 4 and *N* = 6, with *F*_{eff}(*r*) decaying very fast, especially for the former case.

In Fig. 3(b), we show analogous results for a lower density, *r*_{s} = 4, which is close to sodium.^{130} Overall, the behavior 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 dataset 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} i.e., as the density increases (i.e., decreasing the density parameter *r*_{s}), larger numbers of QMC simulated electrons are necessary to achieve convergence of wavenumber-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 [Fig. 3(c)], which corresponds to a moderately coupled system at the margin of the electron liquid regime.^{52} 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**, *ω*);^{44,45,131} 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* ≳ 2*r*_{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 interparticle 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 [Fig. 3(d)]. In this case, the attractive minimum in *F*_{eff}(*r*) becomes even more pronounced and is followed by a local maximum around *r* = 3*r*_{s}. In fact, a very shallow yet significant second minimum can even be discerned around *r* = 4*r*_{s}. These features are a direct consequence of the incipient mesoscopic order in the system upon entering the electron liquid regime,^{26,52} 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 rescaled the forces by a factor of $rs2$ to make them comparable between different densities. Evidently, all depicted datasets exhibit a qualitatively similar behavior. At large distances, they quickly decay toward zero and do not exhibit the long-range Coulomb tail, as discussed at length in Sec. III A. 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.^{132} Yet, in that case, the effective attraction is mediated by phonons, i.e., by the underlying ionic structure. In contrast, the rigid, uniform, and nonpolarizable 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* toward 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 toward 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 similar but 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. First, it constitutes a direct measure for the degree of UEG collective behavior, which gives rise to various interesting phenomena such as the spectrum of density fluctuations.^{45,133} 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.^{93} Similar possibilities have been suggested for semiconductors^{134} and metallic hydrogen^{91} 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. 5(a) 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* ∼ 2*r*_{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 with the presence of a peak in *g*(*r*), which we show in Fig. 5(b) 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 *quasi-exact* 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)–(17) by setting *G*(*q*) ≡ 0. Evidently, the RPA results converge toward Coulomb interaction at small distances. This can be seen particularly well in Fig. 6(b), where the ratio of the respective datasets 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 toward 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.^{135,136} The RPA substantially overestimates the drop in *F*_{eff}(*r*) due to the medium for *r* ∼ *r*_{s}, and it exhibits Friedel oscillations^{137} for *r* ≳ 2*r*_{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^{95} within the effective static approximation^{96} 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 toward zero for *r* ≳ *r*_{s}. At the same time, qualitative differences arise, which deserve further exploration. In particular, the KO expression overestimates the magnitude of the effective electronic attraction at *r* ∼ 2*r*_{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 the three-body correlations are connected to the quadratic density response function, and so on.^{138} 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;^{43,48,139} see also the Appendix for the analysis of the classical limit. In addition, the KO force monotonically converges toward 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 antisymmetry of the true fermionic density matrix under the exchange of particle coordinates is completely omitted. The results are shown as a 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 counterintuitive finding, which calls for an explanation. In particular, it is well understood that fermionic exchange effects should substantially influence the behavior of two electrons at sufficiently small distances *r* at the present conditions. Indeed, if that was not 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 behavior 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 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 histogram 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,71} 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 counteracting *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.^{138} 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, gray, 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 datasets due to the counteracting 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 computation time with the system size *N*. The short-range 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* ≲ 3*r*_{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.^{79,86} Our present findings clearly indicate that the impact of fermionic exchange effects, 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*^{42} that has been introduced in the context of density functional theory.^{140} 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; this 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 linear-response 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. 8(b) 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 behavior of the system. Remarkably, though, *ϕ*_{cl}(*r*) nicely captures the attractive minimum around *r* = 2*r*_{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*).^{45} 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}

We note that *l* and *m* denote the eigenstates of the *N*-body Hamiltonian given by 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\u0302(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**, *τ*).^{44,131,141,142} 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* ∼ 2*q*_{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 dashed-dotted blue curve shows the exact PIMC results from Ref. 44 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* ∼ 2*q*_{F}, i.e., when the wavelength *λ* = 2*π*/*q* becomes comparable to the average interparticle distance.

Based on the dramatic short-range nature of *ϕ*_{eff}(*r*), Dornheim *et al.*^{45} 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^{45}

Let us briefly reiterate 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 given by

where *W*_{RPA}(*q*) is obtained by inserting *g*_{RPA}(*r*) and *ϕ*_{RPA}(*r*) into Eq. (20). Since the impact of exchange–correlation effects onto the kinetic energy is negligible at these conditions, we immediately get the ansatz

Here, the screening function^{94}^{,}*α*(*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.^{143}

*In lieu* of the true PIMC results for *ϕ*_{eff}(*r*), Dornheim *et al.*^{45} 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*^{44} (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. 45 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 toward the exact (dashed-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*) toward 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,117,119} 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 a numerical proof for an effective attraction between two electrons around *r* ∼ 2*r*_{s} for *r*_{s} ≳ 5 and *θ* ≲ 2. We reiterate our earlier point that this does *not* arise due to the mediation through phonons or other ionic effects and that the rigid, nonpolarizable 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.

A 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.^{138} 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 further improvement of dielectric theories^{124,127,144,145} and guide the development of explicitly nonlinear potentials;^{146,147} (ii) the possible role of the attractive minimum in *F*_{eff}(*r*) as a *pairing mechanism* in the UEG^{38} and its connection to the potential emergence of superconductivity in the electron liquid regime^{93} 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* ≲ 3*r*_{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^{86} 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 such as superfluidity^{4,14,148} and supersolid behavior.^{149–152} In addition, we suggest to use the *pair alignment* model^{45} to explain the spectrum of density fluctuations in other systems, with the *original roton feature* in ultracold ^{4}He^{14,148,153} and ^{3}He^{18,154–156} being a promising candidate; (v) finally, we will extend the present study beyond the UEG model, starting with hydrogen in the WDM regime.^{157,158} This will allow us to investigate the effective interaction between all individual components and provide us with unprecedented insights into the behavior of WDM at the nanoscale.

## ACKNOWLEDGMENTS

This work was partially supported by the Center for Advanced Systems Understanding (CASUS), which is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon state government out of the state budget approved by the Saxon State Parliament. The PIMC calculations were carried out at Norddeutscher Verbund für Hoch -und Höchstleistungsrechnen (HLRN) under Grant No. shp00026 and on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### 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 and 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

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β*Φ(

**)[1 −**

*q**G*(

**)] = [1/**

*q**S*(

**)] − 1. Thus,**

*q*In order to further simplify the above, the Fourier space connection between the static structure factor and total correlation function *S*(** q**) = 1 +

*nH*(

**) is employed and the Fourier transformed exact Ornstein–Zernike equation**

*q**H*(

**) =**

*q**C*(

**) +**

*q**nH*(

**)**

*q**C*(

**) is utilized, where**

*q**C*(

**) is the Fourier transformed direct correlation function**

*q**c*(

**). The resulting expression can be trivially inverted back to real space and reads as**

*r*Recall that the exact nonlinear closure equation of the integral equation theory of classical liquids is given by *g*(** r**) = exp[−

*βϕ*(

**) −**

*r**c*(

**) +**

*r**h*(

**) +**

*r**b*(

**)],**

*r*^{106,107}with

*b*(

**) the bridge function. Application of the natural logarithm, disposal of the common**

*r**βϕ*(

**) +**

*r**c*(

**) −**

*r**h*(

**) expression, and substitution for the potential of mean force ultimately yield**

*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*(

**) = 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:**

*r*^{26,134}the substitution of the internal electron–hole propagator by noninteracting 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.^{124,144} Given the availability of a closed-form expression for the bridge function of the classical one-component plasma at arbitrary coupling,^{108,144} $\beta \varphi eff(r)=\beta \varphi effKO(r)\u2212b(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.

## REFERENCES

^{4}gas

^{4}

^{3}He without fixed nodes

^{4}He across the normal-superfluid transition

^{4}He. Path integrals and maximum entropy

^{3}He in two dimensions