Bound States of the Exchange—Correlation Excitons in the Uniform Electron Gas by the Monte Carlo Simulations

: The modiﬁed path integral representation of Wigner functions and the new Monte Carlo approach has been suggested to account for the impact of the interparticle interaction on the Pauli exclusion principle of fermions. This approach also allows to calculate the momentum distribution functions and to reduce the “sign problem” that is inaccessible to the standard path integral Monte Carlo methods. The obtained pair electron–electron distribution functions for the “uniform electron gas” demonstrate the short-range quantum ordering of electrons associated with exchange–correlation excitons. The exchange–correlation exciton is caused by the interaction of electrons with positively charged exchange holes and the excluded volume effect. The developed approach allows one to study the density–temperature range of the exciton arising, existence, and decay. Using the potential of the mean force and semiclassical Bohr–Sommerfeld quantization condition, we have demonstrated the existence of bound states disturbing the Maxwellian distribution and estimated their average energy levels. The exchange–correlation excitons have not been observed earlier in the standard path integral Monte Carlo simulations.


Introduction
Computer simulation is one of the main tools in quantum statistics today, as it allows to explain thermodynamic and transport properties even under extreme conditions, starting from fundamental potentials of interparticle interactions (ab initio simulations), when perturbative methods cannot be applied. The most widespread numerical methods in quantum statistics are quantum Monte Carlo methods (QMC), usually based on the representation of the quantum partition function in the form of path integrals in the coordinate space of particles [1,2]. For example, the Path Integral Monte Carlo (PIMC) method is used for consideration of thermodynamical properties of liquid helium [2], dense hydrogen plasma, electron-hole gas, quark-gluon plasma, etc. [2][3][4].
The main disadvantages of the PIMC method for simulations of Fermi systems are the "fermionic sign problem" [1] and restrictions on the calculation of average values of quantum operators in the phase space. An alternative approach to quantum mechanics based on Wigner functions in the phase space [5,6] allows to overcome these disadvantages. Moreover, it allows to calculate quantum momentum distribution functions and transport properties of matter [3].
Generally speaking, the real sign-variable Wigner function of a pure quantum state can be very irregular due to the principal uncertainty, but the averaging over a phase space cell can result in its smooth behavior and even positive definiteness [6]. In this article, we use the Wigner function averaged in the canonical ensemble to reduce irregularities. Another issue is the "sign problem" arising due to the antisymmetrization of the Wigner function. To account for the Fermi statistics and to overcome the "sign problem" we have modified the path integral representation of Wigner functions and the Monte Carlo method to account for the impact of the interparticle interaction on the Pauli exclusion principle. Our main new contribution to solve this problem is the development of a new approximation for the Wigner function (we used a similar approximation for the density matrix in our previous work [7]), in which exchange determinants are calculated with high accuracy by robust standard algorithms of linear algebra.
As an example, we have considered the one-component plasma consisting of a single sort of charged particles, immersed in a rigid neutralizing background. For electrons, the one-component plasma is often called "uniform electron gas" (UEG) or "jellium". The UEG is a model of simple metals [8,9] and is also the source of data for exchange-correlation functionals in the density functional theory (DFT) [10]. Being a strongly coupled quantum system, the UEG can not be studied by analytical methods of the perturbation-based theories. Therefore, ab initio PIMC approach has been used in this article.
We present pair distribution functions, which demonstrate the short-range quantum ordering of electrons associated in the literature [11,12] with exchange-correlation excitons. The physical reasons for this phenomenon is that the exclusion principle pushes electrons with the same spin away from a reference electron, leaving a positive charge. The observed electron ordering is caused by the interaction of electrons with positively charged exchange holes and the excluded volume, resulting in the formation of exchange-correlation excitons. The obtained charge estimation of the exchange hole shows that its average positive charge is equal to the elementary charge, making the exchange-correlation exciton neutral. So, this short-range screening validates a simple Fermi-liquid model of metals, in which nearly-free electrons are assumed.
The presented results allow us to estimate the density-temperature range of arising, existence, and decay of exchange-correlation excitons. The electron momentum distribution functions (the Wigner functions) demonstrate the disturbance of the Maxwellian distribution by narrow sharp peaks due to the arising of the exciton bound states. These physical effects were not observed earlier in the standard PIMC simulations.
A physical analysis of the bound states was done in the framework of the semiclassical approximation, which is very effective for the consideration of many issues of quantum mechanics and mathematical physics. We started from the analysis of the pair electronelectron distribution functions and related potentials of the mean force [13]. Our analysis, based on the semiclassical Bohr-Sommerfeld quantization condition [14,15], allows us to estimate the energy levels averaged over the canonical ensemble. We also demonstrate the stochastic nature of the energy levels disturbing the Maxwell distribution.

Path Integral Representation of the Wigner Function
In this section, we are going to consider a new representation of the Wigner function, which allows for the numerical simulations of strongly coupled quantum systems of particles in the canonical ensemble. As an example, we consider a 3D two-component Coulomb system of particles consisting of N e electrons and N p heavier positive particles (N e = N p ). The Hamiltonian of the systemĤ =K +Û c contains kinetic energyK and Coulomb interaction energyÛ c =Û c pp +Û c ee +Û c ep contributions. Since the operators of kinetic and potential energy do not commutate, the exact explicit analytical expression for the Wigner function does not in general exist, but can be constructed using a path integral approach [1,2,16] based on the operator identity e −βĤ = e − Ĥ · e − Ĥ . . . e − Ĥ , where = β/M. The Wigner function of the multiparticle system in the canonical ensemble is defined as the Fourier transform of the off-diagonal matrix element of the density matrix in the coordinate representation: [1,2,[4][5][6][16][17][18][19]]. An antisymmetrized Wigner function can be written in the form: The partition function Z for a given temperature T and fixed volume V is defined by the expression: where ρ(x, σ; β) denotes the diagonal matrix elements of the density operatorρ = e −βĤ . Here x = {x e , x p } and σ = {σ e , σ p } are the spatial coordinates and spin degrees of freedom of the electrons and positive particles, i.e., x a = {x 1,a . . . x l,a . . . x N a ,a } and σ a = {σ 1,a . . . σ t,a . . . σ N a ,a }, λ a = 2πh 2 β m a is the thermal wavelength with a, b = e, p and l, t = 1, . . . , N a . The integral in Equation (2) can be rewritten as: where we imply that momentum and coordinate are dimensionless variables p l,aλa /h, and x l,a /λ a related to the electron and positive particle wavelength at a temperature ∼ 1/ (λ a = 2πhβ/(m a M)). The spin gives rise to the spin part of the density matrix (S) with exchange effects accounted for by the permutation operatorsP e acting on coordinates of electrons x (M) and spin projections σ . The sum is taken over all permutations with parities κ P e . In Equations (1) and (3), index m = 0, . . . , M − 1 labels the off-diagonal high-temperature density matrices With the error of the order of 1/M 2 each high temperature factor can be presented in the form , arising from neglecting the commutator 2 /2[K, U c ]. In the limit M → ∞ the error of the whole product of high temperature factors is equal to zero (∝ 1/M) and we have an exact path integral representation of the partition function, in which each particle is represented by a trajectory consisting of a set of M coordinates ("beads"):x ≡ {x  N p ,h }. In the thermodynamic limit, the main contribution in the sum over spin variables of electrons comes from the term related to the equal numbers (N e /2) of electrons with the same spin projection [3,4]. For a heavy positive charge, we account for only the identical permutation, as a respective contribution of exchange interaction is negligible, while for electrons the sum over permutations gives the product of determinants: det (N e /2+1) . On the other hand, the high-temperature density matrix ρ (m) = x (m) |e − Ĥ |x (m+1) can be expressed in terms of two-particle density matrices [3] ρ ab (x l,a , x l,a , This results from the factorization of the density matrix into kinetic and potential parts, ρ ab ≈ ρ K 0 ρ U ab . The off-diagonal density matrix element (4) involves an effective pair interaction, which is expressed approximately via its diagonal elements, where e a and e b are the charges of particles, is the reduced mass of the (ab)-pair of particles and the error function is defined by erf( The discussion of other effective potentials can be found elsewhere [22][23][24]. So, for more accurate accounting of quantum effects, the potential energy U in (1) and (3) will be taken as the sum of pair interactions given by Kelbg potentials and ))/2 [3,4].

Modifications of the Path Integral Measure and Exchange Determinant
In general, the complex-valued integral over ξ in the definition of the Wigner function (1) cannot be calculated analytically and is inconvenient for Monte Carlo simulations. The second disadvantage is that Equations (1) and (3) contain the sign-altering determinant det Ψ(x) , which is the reason for the "sign problem", worsening the accuracy of the PIMC simulations.
To overcome these problems, let us replace the variables of integration x (m) by q (m) (q (0) = q (M) = 0) for any given permutation P by the substitution [18,25]: and make use of the Taylor series of the potential U(x) in the variables ξ. An approximation for the Wigner function arising from accounting for the first two terms of the Taylor expansion around ξ = 0 takes the form of a Gaussian integral, and can be calculated analytically [18,25]. Here, for simplicity let us consider an approximation accounting for the first two terms in the Taylor expansion. Inclusion of the quadratic term is obvious. Below, taking into account the interference of the Coulomb and exchange interactions, we are going to consider two approaches. The first one uses the pair pseudopotentials depending on the phase space variables to satisfy the Pauli exclusion principle without any antisymmetrization. In the second approach, the interparticle interaction is included in the exchange determinant.
The main idea of the first approach can be explained using the system of ideal fermions. We are going to consider the sum over permutations in a pair approximation by accounting for only identical and pair permutations as in [26]. For example, for electrons at M = 1 we have: where the exchange pseudopotential is defined asṽ ). This formula shows that the first corrections accounting for the antisymmetrization of the density matrix result in the appearance of the additional pair exchange pseudopotentialṽ a lt of interparticle interaction [26] As a consequence of replacing the variable, the matrix elements of the density matrix can be rewritten in the form of a path integral over "closed" trajectories {q (0) , . . . , q (M) } with q (0) = q (M) = 0 after the integration over ξ [18,25] and some additional transformations (see [7,27] for details). The Wigner function can be written in the form of the Maxwell distribution with a quantum correction (if we take into account the sum over only identical and pair permutations): where To extend the region of applicability of the pair approximation, the parameter α in pseudopotential (9) was used as an adjustable function of the universal degeneracy parameter of ideal fermions nλ 3 and was chosen as α 2 = 0.00505 + 0.056nλ 3 e . As was shown in [18,25] an analogous repulsive Pauli blocking pair pseudopotential in the phase space prevents the occupation of the same phase space cell by two identical fermions. The comparison of a PIMC momentum distribution function with the Fermi distribution for ideal fermions made in [18,25] has demonstrated that the Pauli blocking pseudopotential accounts for the Fermi statistical effects and provides good agreement with the analytical distribution in wide ranges of fermion degeneracy and momenta.
The difference of the Pauli blocking pseudopotential v e lt (9) developed here, from the one from the papers [18,25], is that it is more accurate, accounting for the impact of the pair interparticle interaction on the Fermi statistics effects by the third exponential function in v e lt (9). On the other hand, to avoid restrictions of the pair exchange approximation, the Wigner function can be written in a form that allows for the sum over all permutations. Presented below is the approximation based on the approach developed in [7,27], which is the direct generalization of (7) [26]: (N e /2+1) and Let us stress that the approximate implementation of the Wigner function used in our simulation and specified in Equation (10) accounts for the contribution of all permutations in the form of the determinant det Φ(x) . Moreover, approximation (10) has the correct limits to the cases of weakly and strongly degenerate fermionic systems. In fact, in the classical limit, due to the factor exp[−π|P e x − x| 2 /M], the main contribution comes from the identical permutation and the differences of potential energies in the exponents, which are equal to zero. At the same time, when the thermal wavelength is of the order of the average interparticle distance and the trajectories are highly entangled, the potential energy depends weakly on permutations, which enables us to replace all permutations P by the identical one E and the differences of potential energies in the exponents are equal to zero.
So, the approximations discussed above show that the Pauli exclusion principle can be realized in two ways: by the Pauli blocking pseudopotential in the phase space (8) or by the determinant in the configurational space (10).
The basic idea of the Monte Carlo method is to replace the integration in Equation (11) with the averaging over samples {x 1 ,x 2 , . . . ,x M } of a random vectorx: where the random quantities x i = (p, x, q (1) , . . . , q (M−1) ) i are drawn from any distribution w(x)/Q (Q = Ω w(x)dx). According to the law of large numbers, if random vectors x i are not correlated, the statistical error is proportional to 1/ √ N and can be estimated using the 3σ-rule.
The samples {x 1 ,x 2 , . . . ,x M } of a random vectorx with a probability density w(x) can be obtained using the Metropolis algorithm. The Metropolis algorithm originates a Markov process by constructing the transition probabilities P(x → x ), so that its stationary distribution is equal to w(x). This algorithm consists of sequential steps, each of them divided into two sub-steps: proposal and acceptance. Suppose the system is in a statē x i , i.e., the random vectorx has a valuex i . On the proposal step, a new random vector x i is proposed. On the acceptance step, this new state can be accepted with a probability A(x i →x i ), thenx i+1 =x i , or if rejected, thenx i+1 =x i . The acceptance probability A(x i →x i ) must be set to satisfy the detailed balance equation, and the most common choice is Since the number of particles in Monte Carlo simulations is limited by available computer resources, it is necessary to consider a finite number of electrons N in a finite volume V, for a given density N/V. To reduce the finite-size effects and reproduce the actual macroscopic system as accurately as possible, the main cubic cell with a side length L = V 1/3 and periodic boundary conditions are usually used. If any "bead" q a,k leaves the main cell during the simulation, it must be replaced with the periodic image entering the cell: where [x] means the floor integer value of a real variable x. In calculations of the distances between "beads" ∆q ab,k the nearest images must be taken: The Monte Carlo algorithm for the calculation is as follows.
Select the type of step randomly: δx-step with a probability P x or δq-step with a probability P q = 1 − P x . If the δx-step is chosen, select a random particle a, and assign x a → x a + ∆x where ∆x is a uniformly distributed random vector in the volume ∼ L 3 /N. If the δq-step is chosen, select a random particle a and a random "bead" k and assign q → q + ∆q where ∆q is a uniformly distributed random vector in a volume ∼ λ 3 /(4πM) 3/2 . The resulting state is set as the proposed state x i ; 3.
Accept the proposed statex i with a probability A(x i →x i ) from (16), or reject. If the proposed state is accepted, setx i+1 =x i . If it is rejected, setx i+1 =x i ; 4.
Calculate the average values for the sample of N steps states via the averaging of Weyl's symbols over the sample with w(x) as a weight function: Repeat all the previous steps for l = 1, 2, . . . , N runs , but instead of the initialization x 1 =x 0 use the last state from the previous run:x 1 = x N steps prev ;

7.
As a result, a sample of average values A l , l = 0, 1, . . . , N runs is obtained. Considering the 0-th run as idle in order to "forget" the initial state, calculate the resulting average value and its statistical error as The probabilities of steps P x and P q = 1 − P x are arbitrary parameters, but they should be selected in order to optimize the convergence of the Metropolis algorithm. We use P q = 1/(M) for the same rate of δx-steps and δq-steps for each particle and each "bead". We set the variations ∆x and ∆q in order to the rate of the accepted steps being 30-70% in the whole range of the considered thermodynamic states. Analogous sampling of the momentum is trivial as it does not require periodic boundary condition.
In our simulations we varied both the particle number in the range N e = 30, . . . , 50, the number of high-temperature factors in the range M = 20 . . . 30, the number of steps N steps ∼10 6 and the number of runs N runs ∼10 2 in order to obtain the convergence of calculated results and reliable estimations of statistical errors. For example, small irregular peaks in momentum distribution functions (MDF) characterize the statistical error of WPIMC calculations.

Simulation Results
Here we present the results of the path integral Monte Carlo simulations in the Wigner approach (WPIMC) of the unpolarized UEG. With that in mind let us start from a two-   Let us stress that for a fixed θ with increasing density, the physical temperature increases and the MDF is approaching the Maxwell distribution (MD) (line r s = 0.25, W Mxw ∼ exp[−(pλ/h) 2 /4π]). Moreover, at large momentum in asymptotic regions the WPIMC MDF and Maxwell momentum distributions practically coincide. At smaller momentum and intermediate density the WPIMC MDFs show two types of peaks disturbing the Maxwell momentum distribution function: small irregular peaks characterizing the statistical error of WPIMC calculations and the narrow high individual peaks, whose physical meaning must be analyzed.

Pair Distribution Functions and Potential of the Mean Force
We are going to consider the peaks in the electron pair distribution functions (PDF) associated with the exchange-correlation excitons [27] and the sharp high picks in the momentum distribution functions (MDF) to check the existence of the corresponding electron bound states. To do this, let us start from the analysis of a potential of the mean force (PMF) [13] defined by the electron-electron PDF. The potential of the mean force ought to reproduce the correct radial distribution function in PIMC simulations.
The potential of mean force w (n) [13] of a system with N classical particles is defined as Above, −∇ j w (n) is the averaged force, i.e., the "mean force" acting on particle j. For n = 2 w (2) is related to the radial distribution function of the system as Here, the PDF g(r) for a binary mixture of electrons and positive particles can be written in the form [13,28], where r = R 1 − R 2 , a, and b label the particle species. Figure 2a presents the WPIMC calculations of the electron-electron PDFs, while Figure 2b shows the corresponding PMFs in atomic units. As previously, for convenient comparison the origins of the distributions and PMFs are shifted along both axes. With increasing density for a fixed θ, the characteristic PDF height becomes smaller while the PDF width becomes narrower. Respectively, the depth of the PMF is increasing while the width of the PMF is decreasing. Peaks in PDFs can be associated with the exchangecorrelation excitons [11,12] arising due to the occurrence of the exchange positively charged hole and the excluded volume effect [29]. This effect is a consequence of the interference effects of Coulomb and exchange interactions of electrons.

Bound States and the Energy Levels
Let us analyze the physical consequences of the above-mentioned behavior of the PMF. The semiclassical approximation used here turned out to be very effective for many issues of quantum mechanics and mathematical physics.
The fundamental difference between a one-dimensional potential well from a threedimensional system is that for a one-dimensional potential well there is always at least one energy level with an even wave function. In the case of a spherically symmetric potential well w 2 , this may not be the case. Let us consider a simple criterion for the existence of bound states [15]: where w (2) min is the minimal value of w (2) and both w (2) min and r 0 are taken in atomic units. Here r 0 is the radius of a spherical potential well. Under this condition, not a single level of the discrete spectrum will appear (as the depth of the potential well is too shallow). Values B presented in the caption to Figure 2b demonstrates that at Θ = 4, bound states can exist only at low enough density corresponding to r s ≥ 0.5. At higher density (r s ≤ 0.5), bound states can be destroyed by "pressure ionization" similar to plasma.
For more accurate semiclassical estimations let us use the Bohr-Sommerfeld condition for a particle in a spherically symmetric field w (2) . The Schrödinger equations for the radial part R(r) of the wave function looks like: where E l is the total energy of the particle, l = 0, 1, 2, . . . is the orbital quantum number [15]. In the semiclassical approximation for a particle in the field w (2) , the Bohr-Sommerfeld condition takes the form: where for any energy E ≥ w (2) min , there are only two turning points, r 1 and r 2 , and where n r is the number of zeros of the radial function [15]. Figure 3 presents the average energy levels obtained according to the semiclassical Bohr-Sommerfeld condition (24) with n r = 0 in the range of densities related to r s from 0.87 up to 1.50. The obtained results show that the bound states with negative energy tend to zero at increasing density and then the bound states disappear. To check the disturbance of the momentum distribution functions by the bound states, let us estimate the related values of momenta. To do this, we use the equation for the spherical potential well pa 0 /h = 2(|E|/Ha)/(r 0 /a 0 ). We can see that the sharp peaks of the MDF at pa 0 /h ∼ 0 and pa 0 /h ∼ 1 correspond to the bound states with |E|/Ha ∼ 0, |E|/Ha ∼ 1 − 2 and r 0 /a 0 ∼ 1 (see Figure 2), which are demonstrated in Figures 1 and 4 for Θ = 4, Θ = 2, Θ = 1.   Figure 4c there are two series of high similar peaks presented by lines related to r s = 0.5, r s = 0.75, r s = 0.87 and r s = 1, r s = 1.12, r s = 1.25, respectively. Let us note that lines r s = 1, and r s = 1.25 confirm clearly the existence of the bound states with eigenvalues E of the order of two Hartrees. Here the electron momentum distribution functions are affected by a strong interaction of charges.
The appearance of the bound states is the physical reason of the peaks in the electron MDFs. These results also reveal the density-temperature range of the exchange-correlation exciton arising, existence, and decay.
To analyze the influence of interparticle interaction on the exchange interaction and the properties of the UEG we have compared the obtained results for Equation (10) (marked by A1) to those (marked by A2) for: φ kt = exp{−π|r kt | 2 /M}. In Figure 4, the MDF for case A2 has no sharp peaks and demonstrates only smooth small oscillations related to statistical errors of the WPIMC calculations. The expression A1 in comparison with approximation A2 accounts for more accurately the Coulomb and Fermi repulsion of electrons, which result in the occurrence of the peaks in MDFs and PDFs.

Discussion and Conclusions
The representation of the Wigner function developed in this paper allows us to account for the impact of the interparticle interaction on the Pauli exclusion principle of fermions, as well as to reduce the "sign problem" and to calculate the momentum distribution functions that are inaccessible to the standard PIMC methods.
The obtained pair electron-electron distribution functions for the UEG demonstrate the short-range quantum ordering of electrons associated in the literature with the exchangecorrelation excitons. The presented results allow us to estimate the density-temperature range of arising, existence, and decay of the exchange-correlation excitons. This fine physical effect was not observed earlier in standard path integral Monte Carlo simulations.
The physical analysis of the electron bound states in excitons was done in the framework of the semiclassical Bohr-Sommerfeld quantization condition applied to the potentials of the mean force arising from the corresponding pair electron-electron distribution functions. We demonstrated the existence of the electron bound states disturbing the Maxwellian distribution and estimated the averaged energy levels.
The electron MDF of the mixture of free electrons and exchange-correlation excitons is characterized by narrow sharp peaks in contrast to the monotonic Maxwellian distribution. A qualitative explanation of the MDF behavior averaged over all electrons is as follows: the MDF of the free electrons and the center of mass of the excitons is Maxwellian, with the average momentum of the order of (5 − 10)pa 0 /h at considered temperatures. As it follows from the behavior of PMF in Figure 2b and the Bohr-Sommerfeld condition, the electron momentum of the bound states in the center of mass is of the order of pa 0 /h∼1. The total electron momentum of the exciton is the sum of the center of mass momentum and the relative one to the center of the mass. So, the electron MDF averaged over the excitons with the center of mass momentum pa 0 /h 1 have to coincide with the Maxwellian distribution, as the bound state peaks in the distribution will be smoothed by the Brownian motion. At the same time, for slow electrons and excitons with the center of mass pa 0 /h∼1, the Maxwellian distribution has to be explicitly disturbed by the bound state peaks, as demonstrated in Figure 4, for the MDF averaged over all electrons.
The noticeable repetition in position of the series of similar sharp peaks in the MDF for neighboring densities and temperatures in independent WPIMC calculations confirms the reliability of the obtained results and the average structure of energy levels.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/universe8020079/s1, Table S1: Table of shifts of the MDF origins along axes x and y for Θ = 4 ( Figure 1). Table S2: Table of shifts of the PDF origins along axes x and y for Θ = 4 (Figure 2a). Table S3: Table of shifts of the potential of mean force origins along axes x and y for Θ = 4 (Figure 2b). ; Table S4: Table of shifts of the MDF origins along axes x and y for Θ = 4 (Figure 4a). Table S5: Table of shifts of the MDF origins along axes x and y for Θ = 2 (Figure 4b). Table S6: Table of shifts of the MDF origins along axes x and y for Θ = 1 (Figure 4c).
Author Contributions: All authors equally contributed to the article. All authors have read and agreed to the published version of the manuscript.