Solution of the ’sign problem’ in pair approximation

Solution of the 'sign problem' in pair approximation

Abstract

The main difficulty for path integral Monte Carlo studies of Fermi systems results from the requirement of antisymmetrization of the density matrix and is known in literature as the ’sign problem’. To overcome this issue the new numerical version of the Wigner approach to quantum mechanics for treatment thermodynamic properties of degenerate systems of fermions has been developed. The new path integral representation of quantum Wigner function in the phase space has been obtained for canonical ensemble. Explicit analytical expression of the Wigner function accounting for Fermi statistical effects by effective pair pseudopotential has been proposed. Derived pseudopotential depends on coordinates, momenta and degeneracy parameter of fermions and takes into account Pauli blocking of fermions in phase space. The new quantum Monte-Carlo method for calculations of average values of arbitrary quantum operators has been proposed. To test the developed approach calculations of the momentum distribution functions and pair correlation functions of the degenerate ideal system of Fermi particles has been carried out in a good agreement with analytical distributions. Generalization of this approach for studies influence of interparticle interaction on momentum distribution functions of strongly coupled Fermi system is in progress.

1 Introduction

Over the last decades significant progress has been observed in theoretical studies of thermodynamic properties of strongly correlated fermions at non-zero temperatures, which is mainly conditioned by the application of numerical simulations (see review [1]). The reason for this success is the possibility of an explicit representation of the density matrix in the form of the Wiener path integrals [2] and making use of the Monte Carlo method for further calculations. The main difficulty for path integral Monte Carlo (PIMC) studies of Fermi systems results from the requirement of antisymmetrization of the density matrix [2]. As result all thermodynamic quantities are presented as the sum of alternating sign terms related to even and odd permutations and are equal to the small difference of two large numbers, which are the sums of positive and negative terms. The numerical calculation in this case is severely hampered. This difficulty is known in the literature as the ’sign problem’.

To overcome this issue some approaches have been developed [1, 3, 4, 5, 6, 7, 8, 9]. For example approaches [3, 5, 6] get ideal Fermi gas with very good accuracy. For interacting fermions the first results [5] have been obtained for the high densities for 33 particles and were now extrapolated to the thermodynamic limit very accurately in [6].

The ’fixed-node method’ [1, 7, 8, 9] is very widely known. However the main result of work [10] is that the ’fixed–node method’ can not reproduce even the well known ideal fermion density matrix and should be considered as an uncontrolled empirical approach for treatment thermodynamics of fermions. The analogous contradictions have been analytically obtained many years ago in [11] from virial decomposition of the many fermion ’fixed-node’ density matrix’.

In this work to treat the ’sign problem’ the new numerical version of the Wigner approach to quantum mechanics allowing studies of thermodynamic properties of the degenerate systems of fermions has been developed. The new path integral representation of quantum Wigner function in the phase space has been obtained for canonical ensemble. Explicit analytical expression of the Wigner function accounting for Fermi statistical effects by effective pair pseudopotential has been proposed. Derived pseudopotential depends on coordinates, momenta and degeneracy parameter of fermions and takes into account Pauli blocking of fermions in phase space. The new quantum Monte-Carlo method for calculations of average values of arbitrary quantum operators has been proposed. To test the developed approach calculations of the momentum distribution functions of the ideal system of Fermi particles has been carried out. Calculated by Monte Carlo method the momentum distributions and pair correlation functions for degenerate ideal fermions are in a good agreement with analytical distributions. Generalization of this approach for studies influence of interparticle interaction on momentum distribution functions of strongly coupled Fermi system is in progress. First results show appearance of long quantum ’tails’ in the Fermi distribution functions.

2 Wigner function for canonical ensemble

Average value of arbitrary quantum operator can be written as its Weyl’s symbol , averaged over phase space with the Wigner function [12, 13]:

 ⟨^A⟩=∫dp\rmdx2πℏA(p,x)W(p,x;β), (1)

where the Weyl’s symbol of operator is:

 A(p,x)=∫dξ2πℏe−i⟨ξ|p⟩/ℏ⟨x−ξ/2|^A|x+ξ/2⟩. (2)

Weyl’s symbols for usual operators like , , , , , etc. can be easily calculated directly from definition (2). The Wigner function of many particle system in canonical ensemble is defined as a Fourier transform of the off–diagonal matrix element of the density matrix operator in coordinate representation:

 W(p,x;β)=Z(β)−1∫dξei⟨p|ξ⟩/ℏ⟨x−ξ/2|e−β^H|x+ξ/2⟩. (3)

Here is the density matrix operator of a quantum system of particles with the Hamiltonian equal to the sum of kinetic and potential energy operators, while , is partition function.

There are well known difficulties in derivation of exact explicit analytical expression for Wigner function as operators of kinetic and potential energy in Hamiltonian do not commutate. To overcome this obstacle let us represent Wigner function in the form of path integral like in the well known case of the partition function. As example let us consider equilibrium a 3D two-component mass asymmetric electron–hole mixture consisting of electrons and holes () [14]. Here is defined as:

 Z(Ne,Nh,V;β)=1Ne!Nh!λ3Neeλ3Nhh∑σ∫Vdxρ(x,σ;β), (4)

where denotes the diagonal matrix elements of the density operator . In equation (4), and are the spatial coordinates and spin degrees of freedom of the electrons and holes, i.e. and , is the thermal wave length, with .

Of course, the density matrix elements of interacting quantum systems is not known (particularly for low temperatures and high densities), but it can be constructed using a path integral approach [2, 16] based on the operator identity , where , which allows us to rewrite the integral in equation (4) as

 ∑σ∫\rmdx(0)ρ(x(0),σ;β)=∫\rmdx(0)…\rmdx(m)…dx(M−1)ρ(1)⋅ρ(2)…ρ(M−1)× ∑σ∑Pe∑Ph(±1)κPe+κPhS(σ,^Pe^Phσ′)^Pe^Phρ(M)∣∣x(M)=x(0),σ′=σ. (5)

In equation (5) the index labels the off–diagonal high-temperature density matrices .

With the error of order arising from neglecting commutator in each high temperature factor can be presented in the form . In the limit the error of the whole product of high temperature factors is equal to zero and this approach gives exact path integral representation of the partition function.

The spin gives rise to the spin part of the density matrix () with exchange effects accounted for by the permutation operators and acting on the electron and hole coordinates and spin projections . The sum is over all permutations with parity and . So thermodynamic values are equal to very small difference between large (of order ) positive and negative contributions giving by the even and odd permutations. The problem of accurate calculation of this difference is the well known sign problem for degenerate Fermi systems. The aim of this work is to develop simple and accurate approach for calculation this difference.

3 Exchange effects in pair approximation

To explain the basic ideas of our approach is enough to consider the simple system of ideal fermions (electrons and holes), so further . The hamiltonian of the system () contains kinetic energy of electrons and holes respectively. Due to the commutativity of these operators the path integral representation of density matrix (5) is exact at any finite number . For our purpose it is enough to consider the sum over permutations in pair approximation at (see [15]):

 ∑σ∑Pe∑Ph(±1)κPe+κPhS(σ,^Pe^Phσ′a)^Pe^Phρ∣∣x1=x(0),σ′=σ =∑σe∑Pe(±1)κPeS(σe,^Peσ′e)ρe∣∣x1e=x(0)e,σ′e=σe ×∑σh∑Ph(±1)κPhS(σh,^Phσ′h)ρh∣∣x1h=x(0)h,σ′h=σh =∑σe{1−∑l

where

 fa;lt=exp(−π|x(0)l,a−x(0)t,a|2λ2a) ~valt=−kTln(1−fa;ltfa;tl)=−kTln(1−δσl,aσt,aexp(−2π|x(0)l,a−x(0)t,a|2λ2a)) (7)

is exchange potential [15]. This formula shows that the first corrections accounting for the antisymmetrization of the density matrix result in the endowing particles by the pair exchange potential . Below to take into account exchange effects in Wigner functions we are going to use analogous pair potential depending on the phase space variables.

4 Path integral representation of Wigner function

Antisymmetrized Wigner function can be written in the form:

 W(p,x;β)=1Z(β)Ne!Nh!λ3Neeλ3Nhh∑σ∑Pe∑Ph(±1)κPe+κPhS(σ,^Pe^Phσ′)∣∣σ′=σ ×∫dξei⟨ξ|p⟩⟨x−ξ/2|M−1∏m=0e−ϵ^Km|^Pe^Ph(x+ξ/2)⟩ (8)

Now replacing intermediate variables of integration in (8) (see (5)) for any permutation :

 x(m)e=(Pexe−xe)mM+xe+q(m)e−(M−m)ξe2M+mPeξe2M x(m)h=(Phxh−xh)mM+xh+q(m)h−(M−m)ξh2M+mPhξh2M (9)

we obtain

 W(p,x;β)=C(M)Z(β)Ne!Nh!λ3Neeλ3Nhh∑σ∑Pe∑Ph(±1)κPe+κPhS(σ,^Pe^Phσ′)∣∣σ′=σ ∫dξ∫dq(1)…dq(M−1)exp{−π⟨ξ|PePh+E|ξ⟩2M+i⟨ξ|p⟩−π|PePhx−x|2M −M−1∑m=0[π|q(m)−q(m+1)|2+ϵU((PePhx−x)mM+x+q(m)−(M−m)ξ2M+mPePhξ2M)]}

where angle brackets in mean the scalar product of vectors and , is unit matrix, while presenting permutation matrix is equal to unit matrix with appropriately transposed columns. Here and further we imply that momentum and coordinate are dimensionless variables like and , where . Here constant as will be shown further is canceled in calculations of average values of operators. As a result, we have a new representation of Wigner function for canonical ensemble in the finite difference form of path integral. Let us note that integration here relates to the integration over the Wiener measure of all closed trajectories , which start and end at . In fact, a particle is presented by the closed trajectory with characteristic size of order in coordinate space. This is manifestation of the uncertainty principle.

Then the Wigner function can be written in the following form:

 W(p,x;β)=C(M)Z(β)Ne!Nh!λ3Neeλ3Nhh∑σ∑Pe∑Ph(±1)κPe+κPhS(σ,^Pe^Phσ′)∣∣σ′=σ ×exp{−π|PePhx−x|2M}∫dq(1)…dq(M−1)exp{−M−1∑m=0[π|q(m)−q(m+1)|2]} ×∫dξexp{−π⟨ξ|PePh+E|ξ⟩2M+i⟨ξ∣∣∣p⟩}. (11)

In this paper we are going to allow for the exchange effects in the pair approximation by effective pseudopotental like have been discussed above (see (6)). So in this approximation Wigner function can be written as:

 W(p,x;β)≈C(M)Z(β)Ne!Nh!λ3Neeλ3Nhh∫dq(1)…dq(M−1) ×exp{−M−1∑m=0[π|qm−qm+1|2]}exp{−M4π|p|2} ×∑σe{1−∑l

The main idea of deriving expression (4) can be explained on example of two electrons in 1D space. For two electrons the sum over permutations consist of two terms related to identical permutation (matrix is equal to unit matrix ) and non identical permutation (matrix is equal to matrix with transposed columns). To do integration in (4) over let us analyze eigenvalues of matrix . For idetical permutation the eigenvalues are equal to each other and are equal to two, while the eigenvalues of matrix related to non identical permutation are equal to zero and two. Integration over for identical permutation is trivial, while for non identical permutation matrix before integration have to be presented in the form , where is diagonal matrix with zero and two as the diagonal elements. Here matrix and inverse matrix are given by the formulas:

 O=∣∣∣11−11∣∣∣.
 O−1=12∣∣∣1−111∣∣∣.

Replacing variables by relation for each pair we can obtain expression (4).

To obtain the final expression we have to approximate delta-function in (4) by the standard Gaussian exponent with small parameter :

 W(p,x;β)≈C(M)Z(β)Ne!Nh!λ3Neeλ3Nhh∫dq(1)…dq(M−1)exp{−M−1∑m=0[π|q(m)−q(m+1)|2]} ×exp{−M4π|p|2}∑σexp(−β∑l

where

 valt≈−kTln{1−δσl,aσt,aexp(−2π|xl,a−xt,a|2M)√M2παexp(−πM|pl,a−pt,a|2(2πα)2)}

Note that the expression (4) contains explicitly term related the classical Maxwell distribution. The others terms account for the influence of exchange interaction on the momentum distribution function. In the limit of small the rescaling by factor regularizes integration over momenta in (1) and allows to use simplified version of effective pair pseudopotential ( is included in small ):

 valt≈−kTln{1−δσl,aσt,aexp(−2π|xl,a−xt,a|2λ2a)exp(−|(pl,a−pt,a)λa|2(2πℏ)2α2)}

Momenta and coordinates are written here in natural units ().

5 Average values of quantum operators

For calculation of average values of quantum operators we are going to use the Monte Carlo method (MC) [17, 18]. To do this we have to use expression (4) presenting the discrete form of path integrals. As a result we obtain final expressions for MC calculations in the following form:

 ⟨^A⟩=⟨A(p,x)⟩w⟨1⟩w. (14)

Here brackets denote averaging of any function with positive weight :

 ⟨g(p,x,q1,…,qM−1)⟩w =∫dp\rmdx∫dq1…dqM−1g(p,x,q1,…,qM−1)w(p,x,q1,…,qM−1). (15)

while

 w(p,x,q1,…,qM−1)=exp{−M−1∑m=0[π|q(m)−q(m+1)|2]}exp{−M4π∣∣∣p∣∣∣2)}, ×∑σexp(−β∑l

Note that denominator in (14) is equal to nominator with , so in (4) is canceled.

Calculations of the average values of quantum operators depending only on coordinates of particles is more convenient and reasonable to carry out in configurational space by standard path integral Monte Karlo method (PIMC). Within considered above approach it can be done if we change the following functions:

 ~w(x,q(1),…,q(M−1))=exp{−M−1∑m=0[π|q(m)−q(m+1)|2+ϵU(x+q(m))]} ×∑σexp(−β∑l

where is defined by equation (7),

 ⟨g(x,q(1),…,q(M−1))⟩~w=∫dx∫dq(1)…dq(M−1)g(x,q(1),…,q(M−1))~w(x,q(1),…,q(M−1)) (18)

and

 ⟨~^A⟩=⟨~A(x)⋅~h(x,q(1),…,q(M−1))⟩~w⟨~h(x,q(1),…,q(M−1))⟩~w. (19)

6 Results of numerical calculations

We define momentum distribution functions and pair correlation functions for holes () and electrons () by the following expressions:

 wa(|p|)=⟨δ(|p1,a|−|p|)⋅h(p,x,q(1),…,q(M−1))⟩w⟨h(p,x,q(1),…,q(M−1))⟩w gab(r)=⟨δ(|x1,a−x1,b|−r)⋅~h(x,q(1),…,q(M−1))⟩~w⟨~h(x,q(1),…,q(M−1))⟩~w (20)

where is delta function, and are types of the particles.

To test the developed approach we have carried out calculations of the path integral representation of Wigner function in the form (4). To extent the region of applicability of pair approximation we have used the small parameter as adjustable function of the universal degeneracy parameter of ideal fermions , namely . Calculations have been done for two hundred particles each presented by twenty beads. Results have been obtained by averaging-out over one million particle configurations. To simplify calculations we fix the number of electrons and holes with the same spin projection equal to and respectively.

For ideal electron–hole plasma Figure 1 shows the momentum distributions and pair correlation functions scaled by ratio of the Plank constant to the electron thermal wavelength () and Bohr radius () respectively. In left column of Figure 1 results of Monte Carlo calculations for electrons and holes are presented by lines 1 and 3, while lines 2 and 4 shows ideal Fermi distributions. Presented distribution functions are normalized to one. Let us note that holes is these calculations are two times heavier than electrons, so the related parameters of degeneracy is times smaller. As it follows from the analysis of Figure 1 agreement of PIMC calculations and analytical Fermi distribution are good enough up to parameter of degeneracy equal to .

It necessary to stress that one of the reason of increasing discrepancy at large degeneracy of femions is limitation on available computing power allowing calculations with several hundred particles in Monte Carlo cell. When parameter of degeneracy is approaching the thermal electron wave length is of order Monte Carlo cell size and influence finite number of particles and periodic boundary conditions becomes significant as was tested by our calculations.

Right column of Figure 1 presents results of Monte Carlo calculations of pair correlation functions . Influence of Fermi repulsion at distance less than thermal wave length is evident enough. At the same time the electron—hole pair correlation functions are identically equal to one as exchange interaction between particle of different type is missing.

Presented results have been obtained in pair exchange approximation described by introduced above the effective pair pseudopotentials. Figure 2 presents contour plots of exchange pair pseudopotentials for parameter of degeneracy equal to . Momenta and coordinates axises are scaled by the electron thermal wave length with Plank constant and factor ten for momentum. As before holes are two times heavier than electrons. As it follows from analysis of Figure 1 the Pauli blocking of fermions in phase space accounting for by these exchange pseudopotentials provides agreement of PIMC calculations and analytical Fermi distribution in wide ranges of fermion degeneracy and fermion momenta, where decay of the distribution functions is about of five orders of magnitude.

7 Conclusion

The new path integral representation of the quantum Wigner function in the phase space has been developed for canonical ensemble. Explicit analytical expression of the Wigner function accounting for Fermi statistical effects by effective pair pseudopotential has been obtained. Derived pseudopotential depends on coordinates, momenta and degeneracy parameter of fermions. The new quantum Monte-Carlo method for calculations of average values of arbitrary quantum operators has been proposed. To test the developed approach calculations of the momentum distribution function and pair correlation functions of the ideal system of Fermi particles has been carried out. Calculated by Monte Carlo method the momentum distributions for degenerate ideal fermions are in a good agreement with analytical Fermi distribution in a wide range of momentum and degeneracy parameter. Generalization of this approach for studies influence of interparticle interaction on momentum distribution functions of strongly coupled Fermi system is in progress. First results show appearance of long quantum ’tails’ in the Fermi distribution functions.

\ack

We acknowledge stimulating discussions with Profs. A.N. Starostin, Yu.V. Petrushevich, E.E. Son, I.L. Iosilevskii, M. Bonitz and V.I. Man’ko.

References

References

1. McMahon J M, Morales M A, Pierleoni C and Ceperley D 2012 Rev. Mod. Phys. 84 1607
2. Feynman R P and Hibbs A R 1965 Quantum Mechanics and Path Integrals (New York: McGraw-Hill)
3. Schoof T, Bonitz M, Filinov A, Hochstuhl D, and Dufty J W, 2011 Contrib. Plasma Phys. 51 687
4. Dornheim T, Groth S, Filinov A and Bonitz M 2015 New Journal of Physics 17 073017
5. Schoof T, Groth S, Vorberger J, and Bonitz M 2015 Phys. Rev. Lett. 115 130402
6. Dornheim T, Groth S, Sjostrom T, Malone F D, Foulkes W M C and Bonitz M 2016 Phys. Rev. Lett. 117 156403
7. Ceperley D 1991 J.Stat. Phys. 63 1237
8. Ceperley D 1992 Phys. Rev. Let. 69 331
9. Militzer B and Pollock R 2000 Phys. Rev. E 61, 3470
10. Filinov V S 2014 High Temperature, 52, 615
11. Filinov V 2001 J. Phys. A: Math. Gen. 34, 1665
12. Wigner E P 1932 Phys. Rev. 40 749
13. Tatarskii V 1983 Sov. Phys. Uspekhi 26 311
14. Filinov V S et al. 2007 Physical review E 75 036401
15. Huang K 1966 Statistical mechanics (Moscow, Mir)
16. Zamalin V M, Norman G E and Filinov V S, 1977 The Monte – Carlo Method in Statistical Thermodynamics, Nauka, Moscow.
17. Metropolis N et al. 1953 Journal of Chemical Physics, 21 1087
18. Hasting W K 1970 Biometrika 57 97
19. Bosse J., Pathak K. N., Singh G. S. 2011 Phys. Rev. E 84 042101
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters