Momentum Distribution Functions and Pair Correlation Functions of Unpolarized Uniform Electron Gas in Warm Dense Matter Regime

Abstract

In this paper we continued our research of the uniform electron gas in a warm dense matter regime, focusing on the momentum distribution functions and pair correlation functions. We use the single–momentum path integral Monte Carlo method, based on the Wigner formulation of quantum statistics to calculate both momentum- and coordinate-depending distributions and average values of quantum operators for many-fermion Coulomb systems. We discovered that the single-particle momentum distribution function deviates from the ideal Fermi distribution and forms the so-called “quantum tails” at high momenta, if non-ideality is strong enough in both degenerate and non-degenerate cases. This effect is always followed by the appearance of the short-range order on pair correlation functions and can be explained by the tunneling through the effective potential wells surrounding the electrons. Furthermore, we calculated the average kinetic and potential energies in the wide range of states, expanding our previous results significantly.

1. Introduction

Over the past decade, the interest in warm dense matter (WDM) has been steadily growing. A WDM regime can be characterized by high values of density and temperature, exceeding by several orders of magnitude those in solids, so WDM is usually non-ideal and degenerate. The study of WDM is required in many astrophysical applications, such as planet interiors [1,2,3], brown and white dwarfs [4,5], and compact stars [6,7]. Furthermore, WDM can be realized in experiments on inertial confinement fusion [8,9] and interaction of intense laser beams with dense plasma [10,11].

To understand this exotic state of matter, the study of its structure and thermodynamic properties is required. One of the most important structural quantities is a pair correlation function (PCF), characterizing the average spatial distribution of particles [12]. The analysis of PCFs gives information about the ordering of particles and hence allows to understand the state of the system: gas-like, liquid-like, or crystalline. Besides that, PCFs can be used for the calculation of average potential energy, which is an important quantity in thermodynamics. Another important thermodynamic quantity is a momentum distribution function (MDF). Characterizing the probability density to have certain momenta values for the particles, an MDF is essential in calculations of many other thermodynamic quantities such as internal energy, scattering parameters, the rate constants of chemical and nuclear reactions, and their thresholds [13,14,15,16].

Due to a strong non-ideality and degeneracy, studying WDM is very challenging. The most powerful and developed analytical methods are usually based on the perturbative expansions or other asymptotic approaches. In this way, the perturbative approximations of the MDFs for a Coulomb system leads to a “quantum tail” exceeding the Maxwell distribution in the weakly non-ideal and non-degenerate case [17,18,19]. However, these analytical methods are not applicable if any appropriate small parameter does not exist, and exactly this situation is common in the case of WDM. Therefore, numerical simulation is the most promising instrument in the studies of WDM, such as quantum Monte Carlo methods for relatively simple systems (particles in external field [20,21], hydrogen and electron-hole plasma [22,23,24,25], etc.) and the density functional theory (DFT)-based methods for more complex and realistic materials.

Being the most powerful and universal numerical tools for studying WDM, DFT-based methods, however, require accurate knowledge of the exchange–correlation (xc) functional, characterizing the exchange and correlation interactions between electrons. To construct such an xc-functional, one has to study a many-electron subsystem in thermodynamic equilibrium and calculate its xc-energy in a wide range of WDM conditions [26]. One of the most fruitful and popular models of the electron subsystem is a uniform electron gas (UEG), which has to be considered under WDM conditions for DFT applications [27,28]. Moreover, the UEG is a well-known quantum analog of a one-component plasma [29] and a simple model of alkali metals [30], so the study of UEG allows one to understand general features of the plasmas in the WDM regime. The UEG consists of electrons on a neutralizing rigid background with a uniformly distributed density of the electrical charges. The thermodynamic state of an unpolarized UEG can be described with two dimensionless parameters—the Brueckner parameter $r_s={(4\pi n{a}_0^3/3)}^{-1/3}$ and the reduced temperature $\theta =kT/E_F$, where $E_F$ is the Fermi energy, n is the electron density, $a_0$ is the Bohr radius. Furthermore, the coupling strength $\mathsf{\Gamma }=\left(e^2\right)/\left(r_{s}kT\right)$ and the degeneracy parameter $\chi =n{\lambda }^3$, where $\lambda =\sqrt{2\pi {\hslash }^2/\left(kT\right)}$ is the thermal wavelength, are often used.

PCFs and MDFs in a weakly non-ideal case have been studied with standard perturbative methods [31,32]. In particular, first-order approximations leads to the “quantum tail” exceeding the Maxwell distribution in a weakly non-ideal and non-degenerate case. However, the UEG in the WDM regime cannot be studied with such approaches and requires ab initio numerical simulations. Most of the quantum Monte Carlo methods uses the coordinate representation of the partition function and do not allow to calculate the MDFs. The recent results for the MDFs of the UEG have been obtained by using the configuration path integral Monte Carlo (CPIMC) method [33] in the representation of occupation numbers. However, only the states with weak non-ideality have been considered ($r_s\lesssim 2$, $\theta \gtrsim 0.0625$), because of a poor convergence at $\mathsf{\Gamma }\gtrsim 1$.

To overcome these difficulties, we have developed a single-momentum path integral Monte Carlo (SMPIMC) method, based on the path integral representation of the Wigner formulation of quantum mechanics. Dealing with the phase space allows one to calculate MDFs, PCFs, and other thermodynamic quantities directly from their definitions. The method has already been applied to many quantum systems: a particle in external potential fields [21], dense hydrogen plasma, and electron–hole plasma [22,25]. In particular, we have shown that the MDFs in plasma differs from the Fermi distribution and have the “quantum tail” in the strongly non-ideal and degenerate case. Furthermore, we have calculated average kinetic, potential, and exchange–correlation energies in our recent paper [34].

In this paper, we continue our research of the UEG, started in [34], and study the MDFs and the PCFs in a wide range of parameters $0.5\le \theta \le 4$, $0.2\le r_s\le 36$ covering the thermodynamic states from almost ideal gas to a deep WDM regime. We have analyzed the dependence of the MDFs and the PCFs on the parameters and discovered that in the regime of strong non-ideality, when the short-range order appears and the UEG becomes liquid-like, the MDFs exceed the Fermi distribution at high momenta and have distinct “quantum tails”. Furthermore, we have expanded our table from [34] for the average kinetic and potential energy to the higher values of $r_s$ up to 36.

2. Theoretical Part

2.1. Uniform Electron Gas

Let us consider a UEG with N electrons at temperature T, contained in the cubic cell with volume V. We assume that the system is unpolarized, and do not take into account spin interactions, so the numbers of electrons with two different spin projections are equal and constant: $N_u=N_d=N/2$. Further, we do not taken into account any spin interactions. In order to reduce the finite-size effects, we apply the periodic boundary conditions (PBC). In addition to the Coulomb interactions in the main cell, one has to take into account the interactions with all periodic images using the Ewald summation technique [35]. As a result, the effective Hamiltonian of the UEG is as follows [36]:

$$\begin{array}{c}\hfill \widehat{H}=\sum _{a=1}^N\frac{{\widehat{\mathbf{p}}}_a^2}{2m}+\sum _{a=1}^N\sum _{b=a+1}^N{e}^2\mathsf{\Psi }\left({\mathbf{q}}_a,{\mathbf{q}}_b\right)+\frac{1}{2}{e}^{2}N{\xi }_M,\end{array}$$

(1)

where m and e are the electronic mass and charge, ${\widehat{\mathbf{p}}}_a$ and ${\mathbf{q}}_a$ are the momentum operator and the coordinate of an a-th electron. The periodic Ewald pair potential is:

$$\begin{array}{c}\hfill \mathsf{\Psi }\left(\mathbf{r},\mathbf{s}\right)=\frac{1}{V\pi }\sum _{\mathbf{G}\ne 0}{G}^{-2}{\mathrm{e}}^{-\frac{{\pi }^2{G}^2}{{\kappa }^2}+2\pi \mathrm{i}\mathbf{G}(\mathbf{r}-\mathbf{s})}-\frac{\pi }{{\kappa }^{2}V}+\sum _{\mathbf{R}}\frac{\mathrm{erfc}\left(\kappa |\mathbf{r}-\mathbf{s}+\mathbf{R}|\right)}{|\mathbf{r}-\mathbf{s}+\mathbf{R}|},\end{array}$$

(2)

and the Madelung constant is:

$$\begin{array}{c}\hfill {\xi }_M=\frac{1}{V\pi }\sum _{\mathbf{G}\ne 0}{G}^{-2}{\mathrm{e}}^{-\frac{{\pi }^2{G}^2}{{\kappa }^2}}-\frac{\pi }{{\kappa }^{2}V}+\sum _{\mathbf{R}\ne 0}\frac{\mathrm{erfc}\left(\kappa |\mathbf{r}-\mathbf{s}|\right)}{|\mathbf{r}-\mathbf{s}|}.\end{array}$$

(3)

Here, $\mathbf{R}={\mathbf{n}}_{d}L$ and $\mathbf{G}={\mathbf{n}}_r/L$, with ${\mathbf{n}}_d$ and ${\mathbf{n}}_r$ having integer components. Note that the Ewald parameter $\kappa$ does not affect the values and should be chosen for efficiency reasons.

The MDF of the UEG describes the probability density of the electrons to have certain values of momentum. Due to the kinematic and dynamic identity of the electrons in the state of thermodynamic equilibrium, a single-particle MDF is of the most interest:

$$\begin{array}{c}\hfill F\left(\mathbf{p}\right)=\int \cdots \int \mathrm{d}{\mathbf{p}}_2\dots \mathrm{d}{\mathbf{p}}_{N}F({\mathbf{p}}_1,\dots ,{\mathbf{p}}_N){|}_{\mathbf{p}={\mathbf{p}}_1},\end{array}$$

(4)

where $F({\mathbf{p}}_1,\dots ,{\mathbf{p}}_N)$ is the N-particle MDF, and the integration is taken over all electrons except for the first one (any other electron can also be chosen).

The PCFs $g_{uu}\left(r\right)$, $g_{dd}\left(r\right)$ and $g_{ud}\left(r\right)$ describe the distribution of distances between the electrons with corresponding spin projections:

$$\begin{array}{c}\hfill g_{uu}\left(r\right)\propto \int \cdots \int \mathrm{d}{\mathbf{q}}_{u3}\dots \mathrm{d}{\mathbf{q}}_{uN}\mathrm{d}{\mathbf{q}}_{d1}\dots \mathrm{d}{\mathbf{q}}_{dN}g({\mathbf{q}}_1,\dots ,{\mathbf{q}}_N){|}_{r=|q_{u1}-q_{u2}|},\\ \hfill g_{dd}\left(r\right)\propto \int \cdots \int \mathrm{d}{\mathbf{q}}_{u1}\dots \mathrm{d}{\mathbf{q}}_{uN}\mathrm{d}{\mathbf{q}}_{d3}\dots \mathrm{d}{\mathbf{q}}_{dN}g({\mathbf{q}}_1,\dots ,{\mathbf{q}}_N){|}_{r=|q_{d1}-q_{d2}|},\\ \hfill g_{ud}\left(r\right)\propto \int \cdots \int \mathrm{d}{\mathbf{q}}_{u2}\dots \mathrm{d}{\mathbf{q}}_{uN}\mathrm{d}{\mathbf{q}}_{d2}\dots \mathrm{d}{\mathbf{q}}_{dN}g({\mathbf{q}}_1,\dots ,{\mathbf{q}}_N){|}_{r=|q_{u1}-q_{d1}|},\end{array}$$

(5)

where $g({\mathbf{q}}_1,\dots ,{\mathbf{q}}_N)$ is the N-particle spatial distribution of the electrons, ${\mathbf{q}}_{ua}$ and ${\mathbf{q}}_{da}$ specify the electrons with different spin projections, and the exact values of the normalization factors depending of $N_u$, $N_d$, and V are not required for our purposes.

2.2. Single Momentum Approach

To obtain the expressions for a single-particle MDF and PCFs required for the numerical calculations, we use the “single momentum approach” presented in [34]. This approach is based on the pseudoprobability density $W(p,q)$ in the phase space known as the Wigner function:

$$\begin{array}{c}\hfill W(p,q)=\int {\mathrm{d}}^{3N}\xi {\mathrm{e}}^{\frac{\mathrm{i}}{\hslash }(p,\xi )}\langle q-\xi /2|\widehat{\rho }|q+\xi /2\rangle ,\end{array}$$

(6)

where $\rho =\exp (-\beta H)$, $\beta =1/\left(k_{B}T\right)$, $(p,\xi )={\sum }_{a=1}^N{\mathbf{p}}_a{\mathit{\xi }}_a$. The quantum states $|q\rangle =|{\mathbf{q}}_1,{\mathbf{q}}_2,\dots ,{\mathbf{q}}_N\rangle$ and $|p\rangle =|{\mathbf{p}}_1,{\mathbf{p}}_2,\dots ,{\mathbf{p}}_N\rangle$ are the N-particle states with certain coordinates and momenta, antisymmetrized according to the Fermi–Dirac statistics:

$$\begin{array}{c}\hfill |q\rangle =\frac{1}{N!}\sum _{P_u,P_d}{(-1)}^{P_u+P_d}|{\mathbf{q}}_{p_1}\rangle |{\mathbf{q}}_{p_2}\rangle \dots |{\mathbf{q}}_{p_N}\rangle ,\\ \hfill |p\rangle =\frac{1}{N!}\sum _{P_u,P_d}{(-1)}^{P_u+P_d}|{\mathbf{p}}_{a_1}\rangle |{\mathbf{p}}_{a_2}\rangle \dots |{\mathbf{p}}_{a_N}\rangle .\end{array}$$

(7)

Here, the sums are taken over all permutations $P_u$ and $P_d$ of the electrons with the positive and negative spin projections having the factor $+1$ for each even and $-1$ for each odd permutation.

The Weyl symbol can be associated with a quantum operator $\widehat{A}$:

$$\begin{array}{c}\hfill A(p,q)={(2\pi \hslash )}^{-3N}\int {\mathrm{d}}^{3N}\xi {\mathrm{e}}^{\frac{\mathrm{i}}{\hslash }(p,\xi )}\langle q+\xi /2|\widehat{A}|q-\xi /2\rangle ,\end{array}$$

(8)

so the average value of $\widehat{A}$ over the canonical ensemble can be calculated via the classical-like expression:

$$\begin{array}{c}\hfill \langle \widehat{A}\rangle =\int \int {\mathrm{d}}^{3N}p{\mathrm{d}}^{3N}qA(p,q)W(p,q).\end{array}$$

(9)

The N-particle momentum and coordinate distribution functions can be obtained from the Wigner function via the integration over coordinates and momenta, respectively [37]:

$$\begin{array}{c}\hfill F({\mathbf{p}}_1,\dots ,{\mathbf{p}}_N)=\int {\mathrm{d}}^{3N}qW(p,q),g({\mathbf{q}}_1,\dots ,{\mathbf{q}}_N)=\int {\mathrm{d}}^{3N}pW(p,q).\end{array}$$

(10)

The integration over ${\mathbf{p}}_2,\dots ,{\mathbf{p}}_N$ leads to the product of delta-functions ${\prod }_{a=2}^N{\delta }^{\left(3\right)}\left({\mathit{\xi }}_a\right)$, so the integrals over ${\mathit{\xi }}_a$ ($a=2,\dots ,N$) disappear. The density matrix with off-diagonal elements for the particles with numbers $b\ne 1$ replaced with zeros is known as the “single–momentum density matrix for particle 1” [34]:

$$\begin{array}{c}\hfill {\rho }_{SM}(q,{\mathit{\xi }}_1)=\langle q-\xi /2|\widehat{\rho }{|q+\xi /2\rangle }_{{\mathit{\xi }}_2=\cdots ={\mathit{\xi }}_N=0}.\end{array}$$

(11)

Integrating it over coordinates q, one obtains the “single–particle $\xi$-distribution function”:

$$\begin{array}{c}\hfill f\left(\mathit{\xi }\right)=\int \mathrm{d}q{\rho }_{SM}(q,{\mathit{\xi }}_1){|}_{\mathit{\xi }={\mathit{\xi }}_1}.\end{array}$$

(12)

The Fourier transform of ${\rho }_{SM}$ results in the single-particle MDF:

$$\begin{array}{c}\hfill F\left(\mathbf{p}\right)=\int \mathrm{d}\mathit{\xi }{\mathrm{e}}^{\frac{\mathrm{i}}{\hslash }\mathbf{p}\mathit{\xi }}f\left(\mathit{\xi }\right).\end{array}$$

(13)

As the macroscopic system is isotropic, the single-particle $\xi$-distribution function and MDF depend only on the modules of the vectors $\mathit{\xi }$ and $\mathbf{p}$, and the 3-dimensional Fourier transform can be reduced to the 1-dimensional sine transform:

$$\begin{array}{c}\hfill F\left(p\right)={\int }_0^{\infty }\mathrm{d}\xi 4\pi {\xi }^{2}f\left(\xi \right)\frac{sin\left(p\xi /\hslash \right)}{p\xi /\hslash }.\end{array}$$

(14)

The expressions for PCFs contains multiple integrals over ${\mathbf{p}}_1,\dots ,{\mathbf{p}}_N$ and can be transformed in a similar way. Due to the product of the delta-functions, ${\prod }_{a=1}^N{\delta }^{\left(3\right)}\left({\mathit{\xi }}_a\right)$, only the diagonal elements of the density matrix remain. The density matrix with non-diagonal elements replaced with zeros is known as the “diagonal density matrix” [34]:

$$\begin{array}{c}\hfill {\rho }_{DG}\left(q\right)=\langle q|\widehat{\rho }|q\rangle .\end{array}$$

(15)

Further integration over $N-2$ coordinates gives the corresponding PCF. Due to the identity of the electrons with the same spin projection, one can choose any pair of them with $\sigma =+1/2$ for $g_{uu}$ or $\sigma =-1/2$ for $g_{dd}$ and the one with $\sigma =+1/2$ and the other with $\sigma =-1/2$ for $g_{ud}$. To improve the convergence of the numerical calculations, we consider each pair of electrons and average the PCF over all pairs:

$$\begin{array}{c}\hfill g_{uu}\left(r\right)\propto \frac{2}{N_u(N_u-1)}\sum _{ua=1}^{N_u}\sum _{ub=ua+1}^{N_u}\left[\prod _{uc\ne ua,ub}\mathrm{d}{\mathbf{q}}_{uc}\right]\left[\prod _{dc}\mathrm{d}{\mathbf{q}}_{dc}\right]{\rho }_{DG}\left(q\right){|}_{r=|q_{ua}-q_{ub}|},\\ \hfill g_{dd}\left(r\right)\propto \frac{2}{N_d(N_d-1)}\sum _{da=1}^{N_d}\sum _{db=da+1}^{N_d}\left[\prod _{uc}\mathrm{d}{\mathbf{q}}_{uc}\right]\left[\prod _{dc\ne da,db}\mathrm{d}{\mathbf{q}}_{dc}\right]{\rho }_{DG}\left(q\right){|}_{r=|q_{da}-q_{db}|},\\ \hfill g_{ud}\left(r\right)\propto \frac{1}{N_u{N}_d}\sum _{ua=1}^{N_u}\sum _{db=1}^{N_d}\left[\prod _{uc\ne ua}\mathrm{d}{\mathbf{q}}_{uc}\right]\left[\prod _{dc\ne db}\mathrm{d}{\mathbf{q}}_{dc}\right]{\rho }_{DG}\left(q\right){|}_{r=|q_{ua}-q_{db}|},\end{array}$$

(16)

where we use the same notations as in (5).

2.3. Path Integrals

For the calculation of the diagonal and single-momentum density matrices by Equations (11) and (15), we use the method of path integrals [38]. Let us consider the density matrix $\rho (q^A,q^B)=\langle q^B|{\mathrm{e}}^{-\beta \widehat{H}}|q^A\rangle$. Decomposing the statistical operator into the product of M high-temperature operators and using $M-1$ complete sets of q-states, one can represent the density matrix in the form of a multiple integral:

$$\begin{array}{c}\hfill \rho (q^A,q^B)=\int \cdots \int {\mathrm{d}}^{3N}{q}^1\dots {\mathrm{d}}^{3N}{q}^{M-1}\prod _{k=0}^{M-1}\langle q^{k+1}|{\mathrm{e}}^{-ϵ\widehat{H}}{|q^k\rangle }_{\begin{array}{c}{q}^0=q^A\\ q^M=q^B\end{array}},\end{array}$$

(17)

where $ϵ=\beta /M$. Using the symbol $\left|\right\{{\mathbf{q}}_a\}\rangle$ for the non-antisymmetrized q-state $|{\mathbf{q}}_1\rangle$$|{\mathbf{q}}_2\rangle$…$|{\mathbf{q}}_N\rangle$ and the symbol $Pq$ for the permutation, one can take into account the Fermi statistics and rewrite the expression (17) as follows [39]:

$$\begin{array}{c}\rho (q^A,q^B)=\int \cdots \int {\mathrm{d}}^{3N}{q}^1\dots {\mathrm{d}}^{3N}{q}^{M-1}\sum _{P_u,P_d}{(-1)}^{P_u+P_d}\hfill \\ \hfill \times \prod _{k=0}^{M-1}\langle \{{\mathbf{q}}_a^{k+1}\}|{\mathrm{e}}^{-ϵ\widehat{H}}|\left\{{\mathbf{q}}_a^k\right\}\rangle {|}_{\begin{array}{c}{q}^0=q^A\\ q^M=P{q}^B\end{array}}.\end{array}$$

(18)

The non-antisymmetrized high-temperature matrix elements are well-known [34]:

$$\begin{array}{cc}\hfill \langle \{{\mathbf{q}}_a^{k+1}\}|{\mathrm{e}}^{-ϵ\widehat{H}}|& \left\{{\mathbf{q}}_a^k\right\}\rangle ={\lambda }_{ϵ}^{-3N}\hfill \\ & \times exp\left\{-\frac{\pi }{{\lambda }_{ϵ}^2}\sum _{a=1}^N{({\mathbf{q}}_a^{k+1}-{\mathbf{q}}_a^k)}^2-\frac{ϵ}{2}\left(U\left(q^{k+1}\right)+U\left(q^k\right)\right)\right\}+O\left(M^{-2}\right),\hfill \end{array}$$

(19)

where ${\lambda }_{ϵ}=\sqrt{2\pi {\hslash }^2ϵ/m}$. Thus, the approximation for the density matrix with an accuracy of $O\left(M^{-1}\right)$ is:

$$\begin{array}{c}\rho (q^A,q^B)\approx {\lambda }^{-3NM}\int \cdots \int {\mathrm{d}}^{3N}{q}^1\dots {\mathrm{d}}^{3N}{q}^{M-1}\sum _{P_u,P_d}{(-1)}^{P_u+P_d}\hfill \\ \hfill \times exp\left\{-\sum _{k=0}^{M-1}\sum _{a=1}^N\frac{m}{2}{\left(\frac{{\mathbf{q}}_a^{k+1}-{\mathbf{q}}_a^k}{ϵ\hslash }\right)}^2-\sum _{k=0}^{M-1}\frac{ϵ}{2}\left(U\left(q^{k+1}\right)+U\left(q^k\right)\right)\right\}{|}_{\begin{array}{c}{q}^0=q^A\\ q^M=P{q}^B\end{array}},\end{array}$$

(20)

where $\lambda =\sqrt{2\pi \hslash \beta /m}$ is the thermal wavelength of electrons and $C\left(M\right)$ is a constant depending on the number of high-temperature terms. At $M\to \infty$ the multiple integral turns into a path integral over all $3N$-dimensional trajectories, and the expression becomes exact:

$$\begin{array}{c}\rho (q^A,q^B)={\lambda }^{-3N}\sum _{P_u,P_d}{(-1)}^{P_u+P_d}{\int }_{\begin{array}{c}q\left(0\right)=q^A\\ q(\beta \hslash )=P{q}^B\end{array}}{\mathrm{D}}^{3N}q\left(t\right)\hfill \\ \hfill \times exp\left\{-\frac{1}{\hslash }{\int }_0^{\beta \hslash }\mathrm{d}t\left[\frac{m}{2}\sum _{a=1}^N{\dot{\mathbf{q}}}_a^2\left(t\right)+U\left(q\left(t\right)\right)\right]\right\}.\end{array}$$

(21)

(Note that for some singular attractive potentials the continuous limit must be considered more carefully [40]).

The Formula (21) contains the fermionic sign problem (FSP) due the sign-alternating permutations and, thus, cannot be used in Monte Carlo simulations directly. To avoid this problem we substitute the variables:

$$\begin{array}{c}\hfill {\mathbf{q}}_a\left(t\right)={\mathbf{z}}_a\left(t\right)+\left(1-\frac{t}{\beta \hslash }\right){\mathbf{q}}_a^A+\frac{t}{\beta \hslash }{\mathbf{q}}_{Pa}^B(a=1,\dots ,N),\end{array}$$

(22)

obtaining the path integral with zero boundary conditions:

$$\begin{array}{c}\rho (q^A,q^B)={\lambda }^{-3N}\sum _{P_u,P_d}{(-1)}^{P_u+P_d}{\int }_{\begin{array}{c}z\left(0\right)=0\\ z(\beta \hslash )=0\end{array}}{\mathrm{D}}^{3N}z\left(t\right)\hfill \\ \hfill \times exp\left\{-\frac{1}{\hslash }{\int }_0^{\beta \hslash }\mathrm{d}t\left[\frac{m}{2}\sum _{a=1}^N{\dot{\mathbf{z}}}_a^2\left(t\right)+U\left(q\left(t\right)\right)\right]\right\}exp\left\{-\frac{\pi }{{\lambda }^2}\sum _{a=1}^N\left[{\left({\mathbf{q}}_{Pa}^B\right)}^2-{\left({\mathbf{q}}_a^A\right)}^2\right]\right\}.\end{array}$$

(23)

Now we assume that all permutations in the potential function can be substituted with the identical one in the WDM regime [23,34]. This simplification allows one to move all permutations into the product of the exchange determinants:

$$\begin{array}{c}\hfill D_{\lambda }^{\sigma }(q^A,q^B)=det\left|exp\left\{-\frac{\pi }{{\lambda }^2}{\left({\mathbf{q}}_a^B-{\mathbf{q}}_b^A\right)}^2\right\}\right|,\end{array}$$

(24)

where $\sigma =u,d$ is the index of the spin projection, $a,b$ are the indices of the electrons with a corresponding $\sigma$. Finally, the path integral representation of the density matrix is:

$$\begin{array}{c}\rho (q^A,q^B)\approx {\lambda }^{-3N}{D}_{\lambda }^u(q^A,q^B)D_{\lambda }^d(q^A,q^B)\hfill \\ \hfill \times {\int }_{\begin{array}{c}z\left(0\right)=0\\ z(\beta \hslash )=0\end{array}}{\mathrm{D}}^{3N}z\left(t\right)exp\left\{-\frac{1}{\hslash }{\int }_0^{\beta \hslash }\mathrm{d}t\left[\frac{m}{2}\sum _{a=1}^N{\dot{\mathbf{z}}}_a{\left(t\right)}^2+U\left(q\left(t\right)\right)\right]\right\}.\end{array}$$

(25)

For numerical applications one have to use a discrete approximation of the path integral. Replacing each continuous trajectory $z\left(t\right)$ with the poly-line $\{z^0,z^1,\dots ,z^{M-1}\}$ having the vertices called “beads“, and the path integral with the multiple integral, we obtain the final expression for the density matrix:

$$\begin{array}{c}\rho (q^A,q^B)\approx {\lambda }^{-3N}{D}_{\lambda }^u(q^A,q^B)D_{\lambda }^d(q^A,q^B)\int \cdots \int {\mathrm{d}}^{3N}{z}^1\dots {\mathrm{d}}^{3N}{z}^{M-1}\hfill \\ \hfill \times exp\left\{-ϵ\sum _{k=0}^{M-1}\left[\frac{m}{2}\sum _{a=1}^N{\left(\frac{{\mathbf{z}}_a^{k+1}-{\mathbf{z}}_a^k}{\hslash ϵ}\right)}^2+\frac{U\left(q^{k+1}\right)+U\left(q^k\right)}{2}\right]\right\}{|}_{\begin{array}{c}{z}^0=0\\ z^M=0\end{array}},\end{array}$$

(26)

where ${\mathbf{q}}_a^k={\mathbf{z}}_a^k+\left(1-k/M\right){\mathbf{q}}_a^A+k/M{\mathbf{q}}_{Pa}^B$.

In the case of the diagonal density matrix ${\rho }_{DG}$, one should set ${\mathbf{q}}_a^A={\mathbf{q}}_a^B={\mathbf{q}}_a$ for $a=1,2,\dots ,N$, and in the case of the single-momentum density matrix ${\rho }_{SM}$, one should make the same choice for $a=2,\dots ,N$ and ${\mathbf{q}}_1^{A,B}={\mathbf{q}}_1\pm {\mathit{\xi }}_1/2$ for $a=1$.

For numerical calculation of MDFs and PCFs we use the developed SMPIMC method, based on the discrete path integral representation (26). The basic ideas and detailed algorithm are presented in Appendix A.

3. Simulation Results

We studied the unpolarized UEG in the range of densities and temperatures $0.2\le r_s\le 36$, $0.5\le \theta \le 4$, that corresponds to the coupling parameter $0.03\lesssim \mathsf{\Gamma }\lesssim 20$ and the degeneracy parameter $0.2\lesssim \chi \lesssim 4$. Therefore, a wide range of states has been covered: from almost ideal to strongly non-ideal and, simultaneously, from an almost classical to degenerate system.

We carried out our simulations for $N=66$ electrons in the cubic cell with PBCs. Such number of electrons has been chosen for conformity with the papers [34,36,41] and is quite enough for $\theta \gtrsim 0.5$. Furthermore, we simulated the unpolarized ideal Fermi gas (IFG) under the same conditions, i.e., the UEG with “turned-off” Coulomb interactions at the same $\theta$ and $r_s$. Comparison of these results with the analytical Fermi distributions allows us to take into account the finite-size effects and control the related systematic errors in the MDF normalization factors. The influence of the finite-size effects on the MDFs manifests itself in the deviation of the calculated MDF for the IFG from the analytical Fermi distribution, which appears at some high value of momentum $p_{FS}$, so the interval $p\lesssim p_{FS}$ contains the reliable data. We estimated roughly the value of $p_{FS}$ from the plot of $F\left(p\right)$ in the logarithmic scale. One can find the values of $p_{FS}$ for each plot in

Table 1. Some parameters and data describing the MDFs and PCFs from the Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5. The values are shown with two significant digits.

${\mathit{r}}_{\mathit{s}}$	$\mathit{\theta }$	$\mathsf{\Gamma }$	$\mathit{\chi }$	$\mathit{\lambda },{\mathit{a}}_0$	$\mathit{L},{\mathit{a}}_0$	${\mathit{p}}_{\mathit{F}\mathit{S}},\mathit{\hslash }/{\mathit{a}}_0$	${\mathit{r}}_{\mathit{m}\mathit{a}\mathit{x}},{\mathit{a}}_0$	${\mathit{r}}_{\mathit{m}\mathit{i}\mathit{n}},{\mathit{a}}_0$
$4.0$	$0.5$	$4.3$	$4.3$	10	26	$0.85$	—	—
$4.0$	$1.0$	$2.2$	$1.5$	$7.4$	26	$1.2$	—	—
$4.0$	$4.0$	$0.54$	$0.19$	$3.7$	26	>1.9	—	—
12	$0.5$	13	$4.3$	31	78	$0.28$	20–22	30–32
12	$1.0$	$6.5$	$1.5$	22	78	$0.42$	20–22	—
12	$4.0$	$1.6$	$0.19$	11	78	$0.65$	—	—
16	$0.5$	17	$4.3$	42	100	$0.21$	26–28	40–42
16	$1.0$	$8.7$	$1.5$	30	100	$0.30$	26–28	—
16	$4.0$	$2.2$	$0.19$	15	100	$0.60$	—	—
28	$1.0$	15	$1.5$	52	180	$0.17$	45–50	70
28	$2.0$	$7.6$	$0.53$	37	180	$0.26$	45–50	70–75
28	$4.0$	$3.8$	$0.19$	26	180	$0.38$	50	—
36	$1.0$	20	$1.5$	66	230	$0.14$	60–65	90
36	$2.0$	$9.8$	$0.53$	47	230	$0.21$	60–65	90–95
36	$4.0$	$4.9$	$0.19$	33	230	$0.28$	60–65	—

. This table also contains the approximate positions $r_{max}$ and $r_{min}$ of the first maximum and minimum of the PCFs $g_{uu}$, $g_{dd}$, and $g_{ud}$.

For low values of $r_s$ ($r_s\lesssim 1$), the UEG is characterized by relatively weak non-ideality ($\mathsf{\Gamma }\lesssim 1$ in the considered range of $\theta$), and the results are quite trivial and do not deserve to be shown in figures. The MDFs completely coincide with the related Fermi distribution for $p\lesssim p_{FS}$. The PCFs are also similar to the ideal case: $g_{uu}$ and $g_{dd}$, and are determined mainly by the exchange repulsion and slowly rise from 0 to 1 with growing r, while $g_{ud}$ is determined only by relatively weak Coulomb repulsion and rises much faster.

The results for $r_s=4.0$ are shown in Figure 1. In cases of $\theta =0.5$ and $\theta =1$, the UEG is significantly non-ideal ($\mathsf{\Gamma }\approx 4.3$ and $\mathsf{\Gamma }\approx 2.2$), while at $\theta =4$ it remains weakly coupled ($\mathsf{\Gamma }\approx 0.54$). Although the difference between the MDFs and the related Fermi distribution does not exceed the statistical error in all cases, the PCFs demonstrate completely non-ideal behavior: $g_{ud}$ rises only slightly faster than $g_{uu}$ and $g_{dd}$ (especially at $\theta =4.0$) because of the same order of the Coulomb and exchange interactions.

The results for $r_s=12.0$ are shown in Figure 2. In the case of $\theta =0.5$, the coupling is very strong ($\mathsf{\Gamma }\approx 13$), and the MDF significantly differs from the Fermi distribution demonstrating an extensive “tail” at high momenta, which is not camouflaged by the finite-size effects. All PCFs have maxima exceeding 1 at $r_{max}$, and $g_{ud}$ grows almost with the same rate as $g_{uu}$ and $g_{dd}$, because the Coulomb repulsion plays a major role. In the case of $\theta =1.0$, the coupling is quite strong ($\mathsf{\Gamma }\approx 6.5$). The difference between MDF and the Fermi distribution does not exceed the statistical error. However, there are weak maxima on the PCFs, and $g_{ud}$ almost coincides with $g_{uu}$ and $g_{dd}$. At $\theta =4.0$, the UEG is significantly non-ideal ($\mathsf{\Gamma }\approx 1.6$), the MDF does not differ from the Fermi distribution, and all the PCFs are the same and do not have any maxima.

The results for $r_s=16.0$ are shown in Figure 3. In the case of $\theta =0.5$ ($\mathsf{\Gamma }\approx 17$), the MDF is quite similar to the previously considered one with $r_s=12$, but the difference between it and the Fermi distribution is more significant. The PCFs are also analogous to the case of $r_s=12$, but the maxima are more pronounced and the minima (at $r_{min}$) begin to appear. In the case of $\theta =1.0$, the coupling is also strong ($\mathsf{\Gamma }\approx 8.7$), and the MDF begins to differ from the Fermi distribution. The PCFs have distinct maxima at $r\approx 28a_0$. In the case of $\theta =4.0$, the coupling strength $\mathsf{\Gamma }\approx 2.2$, and the MDF coincides with the Fermi distribution within the statistical error. The PCFs do not have any maxima yet, and there are no differences between $g_{uu}$, $g_{dd}$, and $g_{ud}$, so the exchange interaction is negligible in comparison with the Coulomb one.

The results for $r_s=28.0$ are shown in Figure 4. In the case of $\theta =1$, the UEG is strongly non-ideal ($\mathsf{\Gamma }\approx 15$), and the MDF differs from the Fermi distribution very distinctly. All PCFs are almost the same and have maxima and minima. In the case of $\theta \approx 2$, the coupling strength $\mathsf{\Gamma }\approx 7.6$ and non-ideality is quite strong. The MDF begins to deviate from the Fermi distribution, while the PCFs $g_{uu}$, $g_{dd}$, and $g_{ud}$ are exactly the same and have distinct maxima. In the case of $\theta \approx 4.0$, the coupling strength $\mathsf{\Gamma }\approx 3.8$, and the MDF, slightly differ from the Fermi distribution. The maxima of the PCFs are small but distinguishable.

The results for $r_s=36.0$ are shown in Figure 5. In the case of $\theta =0.5$, the coupling strength is very strong ($\mathsf{\Gamma }\approx 20$). The difference between the MDF and the Fermi distribution becomes more significant. The first maxima and minima of the PCFs become more pronounced, also, the second maxima starts to appear at $r\approx 120a_0$; however, the number of electrons should be increased for better resolution. In the case of $\theta \approx 2$ ($\mathsf{\Gamma }\approx 9.8$), the MDF has quite a distinct “tail”. The PCFs $g_{uu}$, $g_{dd}$, and $g_{ud}$ are completely the same and have maxima and hardly distinguishable minima. In the case of $\theta \approx 4.0$, the coupling strength $\mathsf{\Gamma }\approx 4.9$, the MDF slightly differs from the Fermi distribution and there are maxima on the PCFs.

In addition, we calculated the average kinetic energy $E_{kin}$ and potential energy $E_{pot}$ of the UEG, supplementing our results from [34] with the strongly coupled states of UEG. These results are presented in

Table 2. $E_{IFG0}$—internal (kinetic) energy of the macroscopic IFG obtained from the analytical Fermi distribution; $E_{IFG}$—internal (kinetic) energy of the finite-size IFG obtained from the SMPIMC; $E_{kin}$ and $E_{pot}$—kinetic and potential energy of the UEG. Two significant digits of the statistical error $3\sigma$ are written in brackets. All energies are given per one electron.

${\mathit{r}}_{\mathit{s}}$	$\mathit{\theta }$	${\mathit{E}}_{\mathit{IFG}0}$, Ha	${\mathit{E}}_{\mathit{IFG}}$, Ha	${\mathit{E}}_{\mathit{kin}}$, Ha	${\mathit{E}}_{\mathit{pot}}$, Ha
$4.0$	$0.5$	$0.118$	$0.118\left(02\right)$	$0.121\left(02\right)$	$-0.22386\left(12\right)$
$4.0$	$1.0$	$0.195$	$0.197\left(05\right)$	$0.195\left(05\right)$	$-0.21480\left(10\right)$
$4.0$	$4.0$	$0.702$	$0.702\left(18\right)$	$0.699\left(18\right)$	$-0.17822\left(14\right)$
12	$0.5$	$0.0131$	$0.0131\left(01\right)$	$0.0141\left(03\right)$	$-0.081073\left(40\right)$
12	$1.0$	$0.0217$	$0.0216\left(03\right)$	$0.0224\left(04\right)$	$-0.079157\left(35\right)$
12	$4.0$	$0.0780$	$0.0793\left(18\right)$	$0.0789\left(18\right)$	$-0.070567\left(24\right)$
16	$0.5$	$0.00735$	$0.00735\left(05\right)$	$0.00841\left(12\right)$	$-0.061467\left(25\right)$
16	$1.0$	$0.0122$	$0.0123\left(02\right)$	$0.0126\left(02\right)$	$-0.060535\left(18\right)$
16	$4.0$	$0.0439$	$0.0438\left(10\right)$	$0.0446\left(09\right)$	$-0.054942\left(17\right)$
28	$1.0$	$0.00399$	$0.00402\left(02\right)$	$0.00464\left(03\right)$	$-0.035732\left(12\right)$
28	$2.0$	$0.00738$	$0.00735\left(08\right)$	$0.00755\left(09\right)$	$-0.034815\left(10\right)$
28	$4.0$	$0.0143$	$0.0143\left(02\right)$	$0.0149\left(02\right)$	$-0.033396\left(10\right)$
36	$1.0$	$0.00241$	$0.00240\left(00\right)$	$0.00314\left(00\right)$	$-0.028246\left(08\right)$
36	$2.0$	$0.00446$	$0.00446\left(04\right)$	$0.00480\left(04\right)$	$-0.027628\left(10\right)$
36	$4.0$	$0.00867$	$0.00869\left(11\right)$	$0.00894\left(12\right)$	$-0.026520\left(05\right)$

and Figure 6 and compared with the internal (kinetic) energy of the IFG, calculated with the SMPIMC method ($E_{IFG}$) and obtained from the analytical Fermi distribution directly ($E_{IFG0}$). The kinetic energy rapidly decreases by many orders of magnitude as the parameter $r_s$ increases, similar to the ideal case (for the Boltzmann gas $E_{kin}\propto r_s^{-2}$). $E_{kin}$ starts to exceed $E_{IFG}$ at high values of $r_s$ due to the effect of “quantum tails”; for example, $E_{kin}$ becomes almost $30\%$ higher than $E_{IFG}$ at $r_s=36$, $\theta =0.5$. The potential energy is negative and rapidly falls as $r_s$ grows, while the dependence of $E_{pot}$ on $\theta$ at fixed $r_s$ is more sharp at low values of $r_s$ and becomes weaker at high values (for the Boltzmann gas $E_{pot}\propto -r_s^{-1}$).

4. Discussion

First of all, let us discuss the phenomenon of “quantum tails” of the MDFs and make some suggestions about its physical nature. As it follows from Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5, the single-particle MDF of the UEG starts to differ from the Fermi distribution when the system becomes strongly non-ideal. The value of the coupling parameter $\mathsf{\Gamma }$, at which the “tail” appears, depends on the degeneracy rate: in the case of $\theta =0.5$, (significant degeneracy) it appears at $\mathsf{\Gamma }$ from 8 to 9 ($r_s\approx 8$), but in the case of $\theta =4.0$, the $\mathsf{\Gamma }\approx 5$ is quite enough ($r_s\approx 36$).

At the same time, the deviation of the MDFs from the Fermi distribution is always followed by the appearance of the visible maxima on the PCFs, and with increasing maxima, the deviation also grows. This observation points at a possible physical reason of “quantum tails”: with increasing non-ideality, the UEG becomes more liquid-like and the short-range order appears. This results in the potential wells surrounding the electrons at the average distance $r_{max}$, so the spatial localization of the electrons increases. According to the basic principles of quantum mechanics, the localization of the electrons in the momentum space is weakened and the MDFs become wider. When the coupling increases, the potential wells become deeper and the effect grows.

Furthermore, this hypothesis can explain the fact that the “quantum tail” of the MDF appears at lower values of $\mathsf{\Gamma }$ when the degeneracy is low and at higher values in the opposite case. In fact, at low $\theta$, the thermal wavelength $\lambda$ becomes of the order of the interparticle distance, so the electron can escape from the potential well more easily via the quantum tunneling and becomes less localized. Note that the definition of the coupling strength $\mathsf{\Gamma }$ is based on the classical value of the average kinetic energy $kT$ and gives an overestimated value of the actual non-ideality in the degenerate case, which should be characterized by the quantum average kinetic energy between $kT$ and $E_F$. Thus, at the same values of $\mathsf{\Gamma }$, a non-degenerate system (high $\theta$) is more non-ideal than a degenerate one (low $\theta$), and the ”quantum tail” effect is more significant.

Secondly, the difference between PCFs $g_{uu}$, $g_{dd}$, and $g_{ud}$ becomes lower with increasing $r_s$ even at low values of $\theta$ (degenerate case). As the exchange repulsion is applied only to the electrons with the same spin projections, while the Coulomb interaction is universal, we can conclude that the PCFs of the degenerate UEG are basically conditioned by the last one. This reasoning is confirmed by the fact that the value of $r_{max}$ is almost independent of $\theta$ when $r_s$ is fixed. However, the value of the maxima depends strongly on $\theta$, following the coupling strength $\mathsf{\Gamma }$.

5. Conclusions

In this paper, we have continued our research of the UEG and study the MDFs and the PCFs using the SMPIMC method in a wide range of temperatures and densities: $0.5\le \theta \le 4$, $0.2\le r_s\le 36$. Thus, different thermodynamic states from an almost ideal gas (Boltzmann of Fermi) to a strongly non-ideal system with weak and significant degeneracy have been considered.

Firstly, we have discovered that the single-particle MDFs start to exceed the Fermi distribution at high momentum values, when the UEG is strongly non-ideal. In the degenerate case $\theta =0.5$, this occurs at $r_s\gtrsim 8$ ($\mathsf{\Gamma }\gtrsim 9$), in the semi-degenerate case $\theta =1$—at $r_s\gtrsim 16$ ($\mathsf{\Gamma }\gtrsim 8$), in the almost classical case $\theta =4$—at $r_s=36$ ($\mathsf{\Gamma }\gtrsim 5$). The deviation has the form of the “quantum tail” at high momenta and is reliable at $p\lesssim p_{FS}$, where the value of $p_{FS}$ is conditioned by the finite-size effects.

We explain the phenomena of “quantum tails” as follows. The deviation of the MDFs from the Fermi distribution is always followed by the appearance of the short-range order on the PCFs; so, with increasing coupling strength, the UEG becomes more liquid-like and the potential wells surrounding the electrons appear. As a result, the spatial localization of the electrons increases, so the localization in the momentum space becomes weaker and the MDFs become wider. Besides, at high degeneracy, the electrons are tunneling from the potential wells more easily, so the “quantum tail” appears at higher values of $\mathsf{\Gamma }$ in comparison with the non-degenerate case.

Secondly, we have obtained that the difference between the PCFs of the electrons with the same and opposite spin projections becomes negligible at high values of $\mathsf{\Gamma }$ in the degenerate case. Thus, the spatial structure of a strongly non-ideal UEG is basically conditioned by the Coulomb interaction rather than the exchange repulsion.

Finally, we have calculated the average kinetic and potential energies for $r_s$ up to 36, supplementing our results from [34] significantly.