Computer simulation is one of the main tools in studies of thermodynamic properties of many-particle strongly coupled systems under extreme conditions. Some of the most powerful numerical methods for simulation of quantum systems are Monte Carlo methods, based on path integral (PIMC) formulation of quantum mechanics  . These methods use path integral representation for partition function and thermodynamic values such as average energy, pressure, heat capacity etc. PIMC methods are widely used for studying dense hydrogen plasma  , electron gas in metal    , electron-hole plasma  in semiconductors, superfluidity    and even quark-gluon plasma    .
Unfortunately PIMC methods cannot cope with problem of calculation of average values of arbitrary quantum operators in phase space or momentum distribution function, while this problem may be central in studies of kinetic properties of matter. Quantum effects may strongly disturb the momentum Maxwellian distribution function  and are important in studies of kinetic properties of matter at low temperatures and under extreme conditions, when particles are strongly coupled and perturbative methods cannot be applied. The mentioned before PIMC methods or even recently developed “Configuration Path Integral Monte Carlo” approach    up to now cannot also solve kinetic problem.
The Wigner formulation of quantum mechanics in phase space allows more natural considering not only kinetic but also thermodynamic problems. Methods for phase space treatment of single-particle quantum dynamics in Wigner approach in microcanonical ensemble were proposed in      . Recently a generalization of these one-body Wigner Monte Carlo methods to the many-body distinguishable particles has been done in  . Wigner dynamics of relativistic particle was studied in  . However for theoretical studies of kinetic properties the ab initio many particle approaches in phase space are required.
In this paper we continue developing the ab initio path integral Monte Carlo approach in phase space by using the basic ideas suggested before in  . Here we have generalized the path integral representation for Wigner function to three-dimensional many-particle quantum system with Boltzmann and Fermi statistics. Here we present also alternative single momentum approach for calculations of momentum distributions and quantum observables. This approach allows avoiding harmonic approximation at consideration of degenerate systems of interacting fermions and partially overcoming the well-known sign problem for degenerate Fermi systems. These two methods were tested on some simple models like particles in external potential field and ideal Fermi gas.
2. Wigner Function in Canonical Ensemble
Average value of arbitrary quantum operator can be written as Weyl’s symbol averaged over phase space with Wigner function  :
where the Weyl symbol of operator is :
Weyl symbols for common operators like , , , , , etc. can be easily calculated directly from definition (2). The Wigner function of the many particle system in canonical ensemble is defined as a Fourier transform of the off-diagonal matrix element of density matrix in coordinate representation:
where is partition function of canonical ensemble of N particles, is inverted temperature, V―volume, , , etc. Hamiltonian of the system consists of kinetic and potential parts: , where , ―potential energy of particles, interacting with each other and/or external field. For simplicity we will consider system containing identical particles of mass . Generalization on systems with several kinds of particles is direct. The Wigner function being the analogue of the classical distribution function in the phase space has a wide range of applications in quantum mechanics and statistics.
For spinless bosons or fermions the density matrix in canonical ensemble looks like
where symbol denotes particle permutations, corresponds to Bose statistics, ―to Fermi statistics. Since operators of kinetic and potential energy in Hamiltonian do not commutate, the exact explicit analytical expression for Wigner function does not exist. To overcome this difficulty let us represent Wigner function similarly to path integral representation of the partition function   . So we represent statistical operator as product of large number of operators formally related to high temperature and write down Wigner function in form of multiple integral of the product of the high temperature density matrices:
In the limit the high temperature density matrices can be easily calculated with accuracy up to  :
where we use abbreviation , and is the thermal wavelength at high temperature . Thus expression for Wigner function takes the form:
where is the thermal wavelength at temperature . In continuous limit this expression turns into path integral over trajectories starting at point and ending at , where is dimensionless parameter, which may be interpreted as “imaginary time”  :
Measure of the path integral depends on Fourier variables and positions . To get rid this dependence we change variables from trajectories to closed dimensionless trajectories :
Taking into account new boundary conditions we obtain the desired path integral representation of Wigner function:
where is given by (9), symbol denotes permutation of identical particles and means integration over trajectories. In fact, this expression means that Wigner function is Fourier transform over of symmetrized or antisymmetrized density matrix, which is represented by path integral. Unfortunately this Fourier transform cannot be calculated analytically in general case even for non-interacting bosons or fermions.
3. Harmonic Approximation
Momenta in expression (10) for Wigner function are connected with other variables through -dimensional Fourier transform, which is not integrable analytically or numerically in general case. Exclusions are the linear or harmonic potentials, when power of variable is not more than two. So let us take the approximation for potential arising from its expansion into series. For simple estimations let us consider only term related to identical permutation describing system with Boltzmann statistics. Thus in harmonic approximation we expand potential energy into Taylor series up to second order in :
where we denote summation over repeated indices from 1 to and over from 1 to 3. If we take expansion (11), the expression for Wigner function (10) takes form of generalized Gaussian integral over variable and expression for Wigner function in harmonic approximation can be written in the following form:
Here we introduced scalar functionals and , as well as functionals and of first and second potential derivatives having form of - dimensional vector and matrix correspondingly. They depends on trajectories and positions :
Note that the first term in exponent (12) looks similarly to Maxwell distribution in classical statistics. The major difference is in matrix and cosine which are equal to identity matrix and unit for non-interacting particles. This matrix and cosine provide correlation of particle momenta with each other and with their positions.
The expression for Wigner function (12) is obtained under assumption that potential energy is expandable in Taylor series of second order with a good accuracy. Let us discuss the legality of such assumption in the simplest 1D case. Primarily, note that the exponent (10) contains variable in three positions. The first one is in Fourier term , which makes momenta correlated with other dynamical variables. The second one is in gaussian-like term . The third one is in integral term, where is argument of potential:
. Gaussian term provides rapid decaying when
increases, so main contribution comes from near . Argument of potential function contains multiplied on , which is modulo less than 0.5. Using mean value theorem, we can roughly get symbolic estimation of these integrals in exponent (10) as
where is certain point of trajectory. Numerical value of this integral rapidly decreases: when it equals , and . Thus we expect negligible contribution of high order Taylor terms in potential expansion. Numerical calculations for some potentials, done in this work, confirm this assumption.
For calculation of we are going to use the harmonic approximation of the Wigner function in path integral representation (12). For making use of the Monte Carlo method (MC)  we have to represent path integrals in discrete form of multiple integral like (7). As a result we obtain expression for MC calculations in the following form:
Here brackets denote averaging of any function with positive weight :
4. Single-Momentum Approach
Now we are going to introduce another approach to calculation of momentum distribution functions and quantum averages. As it was discussed above, precise path integral representation (10) is useless for practice due to -dimensional Fourier transform, which cannot be calculated even for a few particles. However we can overcome this difficulty in many important cases.
Let us consider some operator , whose Weyl symbol is a sum of equal single-particle symbols: . Due to the fact that all particles are identical, average values on ensemble of for arbitrary particles and are equal: . Therefore we can use Weyl symbol instead of and integrate Wigner function over other momenta. Thus in single-momentum formalism:
Here we introduce single-momentum Wigner function as integral of over all momenta except , which we will denote by further:
As momenta in definition of Wigner function appear only in Fourier transform , integrals on them gives delta functions for . This exterminates all except , which now is denoted as simple . Thus formula for single-momentum Wigner function is
i.e. it is 3-dimensional Fourier transform of single-momentum density matrix:
Here is the same functional on trajectories as in (13), P and D are potential functional and determinant of permutations:
Such reasoning would not change if we calculate average value of operators of the form or . To summarize, in single-momentum approach we are able to deal with quantum operators of particular form only:
so average value of it is:
with Weyl symbol depending only on first particle’s momentum . In fact many quantum operators being interesting for statistical physics have such single-momentum form: Hamiltonian, kinetic and potential energy, pair correlations and correlations between position and momentum, momentum distribution function etc.
Now let us apply the basic ideas of this approach to the degenerate system of interacting fermions. In degenerate system average distance between fermions is lesser than the thermal wavelength and trajectories in path integrals (see (10), (9)) are strongly entangled. This is the reason that permutations can not strongly affect the potential energy in comparison with the case of identical permutation.
So can be presented as energy related to identical permutation and small difference , which, for example, for 1D Coulomb system is of order
, where is 1D density of particles. Then
for the first terms of the perturbation series on this parameter one can obtain the following simplification: .
Now all permutations in (10) are connected with variables and and can be beared out of the path integral . Moreover, one can gather all permutations into Slater determinant of matrix:
In conclusion of this section we consider practical application of the approach. Making use of the Monte Carlo again we have to represent path integral for single-momentum density matrix (20) in discrete form, i.e. in form of multiple integral over “intermediate quantum coordinates” . Firstly, by standard Monte-Carlo procedure we are able to calculate Fourier preimage of averaged single-momentum operator of form (22):
which depends on momentum by Weyl symbol and only. Here brackets denote averaging of any function with positive weight :
So we obtain in form of distribution function. Then we can do Fourier transform over variable numerically, using Fast Fourier Transform algorithms. After that it is easy to average expression over momenta to obtain desired average value .
To calculate average values of simpler quantum operators like Hamiltonian containing single-momentum operator as separate term we can use easier way. The terms which do not depend on momentum can be calculated by averaging (25) over without Fourier transform. The terms which are functions of can be calculated from single-momentum distribution function:
obtained by Fourier transform of single-momentum density matrix integrated over all positions .
5. Numerical Results
Let us consider some results obtained by both approaches described above. To test the limits of applicability of harmonic approximation we have carried out calculation of thermodynamic values for single quantum particle in some strongly unharmonic potential fields. We have considered and cases and applied single momentum approach to calculations of momentum distributions. Finally, we tested single momentum approach on systems of non-inte- racting strongly degenerated fermions.
Single particle in one dimension As the first example we consider one particle in -power potential of the fourth order, i.e. unharmonic oscillator:
where and are adjusted parameters of potential. When this potential describes harmonic oscillator and formula (12) is exact. If does not equal zero, oscillator is unharmonic and can be used for verification of the harmonic approximation. The Figure 1 shows potentials (29) for and different values of parameter when anharmonicity is significant.
Left plot of Figure 2 shows dependence of full energy on inverted temperature for potential . Results of harmonic approximation are marked by symbols, while results of usual PIMC method are represented by curves. Results of both methods agree well with each other. When temperature is low (from to ) the energy tends to constant value, namely to the ground state energy .
Figure 1. Potential for and different values of unharmonic parameter (1), (2), (3) and (4).
Right plot of Figure 2 shows dependence of squared energy on inverted temperature for potential . From analysis of this figure it follows that harmonic approximation gives correct enough values of even for the ground state. Weyl’s symbol of is polynomial of the 8 order, so the main contribution to comes from “wings” of Wigner function calculated with larger statistical error. This is the reason of some deviation of obtained results especially at low temperature. However ground state is reached earlier and this deviations are not important.
Dependencies of average kinetic and potential energies on inverted temperature are represented on the left and right plots of Figure 3. Note that usual PIMC method has encountered difficulties at calculation of these val-
Figure 2. Left plot: Dependence of average energy on inverted temperature for potential . Parameters for all cases, (1,2), (3,4), (5,6) and (7,8). Symbols refer to harmonic approximation calculations, curves―to usual PIMC method. Right plot: The same for average square of energy .
Figure 3. Left plot: Dependence of average kinetic energy on inverted temperature for potential . Parameter for all cases, (1), (2), (3) and (4). Right plot: The same for average potential energy . Results have obtained by harmonic approximation method and PIMC.
ues as it calculates derivatives of partition function .
As a second example we consider -potential of the fourth order:
where is characteristic energy of potential. Dimensionless parameter specifies contribution of the cubic term. Figure 4 shows potentials (30) for different values of this parameter. With gain of the depth of the potential well rapidly increases and the potential well becomes asymmetric.
Left plot of Figure 5 shows dependence of average energy on inverted temperature. As before symbols relate to harmonic approximation, while curves relate to usual PIMC method. Good agreement between results of these methods is evident. With increasing parameter the average energy becomes lower and approaches ground state energy, which is of order of the potential minimum (see Figure 4).
Right plot of Figure 5 shows dependence of average squared energy on inverted temperature. Dependencies of average kinetic and potential energies on inverted temperature are represented on left and right plots of Figure 6. One can see weak influence of parameter on kinetic energy, while the potential energy rapidly decrease when gains. Thus contribution of potential energy into the full energy is dominant at low temperature. At very low temperatures some points on left plot of Figure 5 are missed for large values of parameter when expression below for the most trajectories becomes negative:
Figure 4. Power potential for different values of parameter (1), (2), curve (3)― (3), curve (4) ― (4).
Figure 5. Left plot: Dependence of average full energy on inverted temperature for different values of parameter (1, 2), (3, 4), (5, 6) and (7, 8). Right plot: The same for average square of energy . Results have been obtained by harmonic approximation method and PIMC.
Figure 6. Left plot: Dependence of average kinetic energy on inverted temperature for different values of parameter (1), (2), (3) and (4). Right plot: The same for potential energy . Results have obtained by harmonic approximation method and PIMC.
and functional can be imaginary. This situation occurs at low temperatures for potentials with negative second derivative like in potential (30). However this happens at very low temperature when the system is already in the ground state.
As the third example we consider one-dimensional soft-core Coulomb potential (SCC):
Here is electrical charge, while parameter is characteristic length of potential. Sometimes the SCC potential is used in numerical calculations instead of true Coulomb potential to avoid Coulomb divergency  . Figure 7 shows SCC potential for different values of . Curves refer to , and , where is Bohr radius. It coincides with Bohr radius in three dimension  . Left plot of Figure 8 shows dependence of average energy on inverted temperature for SCC potential. Results of both methods are in good agreement with each others. The same agreement takes place for average squared energy represented on right plot of Figure 8. At high-temperatures functions and are not affected by parameter , while at low- temperature the strong dependence on is evident.
Finally, dependencies of average kinetic and potential energies on inverted temperature are represented on left and right plots of Figure 9. To the contrary of the potential energy the kinetic energy is not significantly affected by parameter .
It should be noted, that some points are missed when . The reason is the same as for potential . This happens at very low temperature, when system is already in ground state. However when harmonic approximation fails to calculate ground state, as the second derivative of potential takes a large value near origin of coordinates.
Single particle in three dimensions Let us consider one particle in spherically symmetrical fourth order potential, i.e. unharmonic oscillator:
where and are adjusted parameters of potential, is length of radius-vector . When this potential describes symmetrical - harmonic oscillator and harmonic approximation (12) is exact. If does not equal zero, oscillator is unharmonic and can be used for verification of the devel-
Figure 7. Soft-core Coulomb potential for different values of (1), (2), (3). Here is Bohr radius, -Hartry energy.
Figure 8. Left plot: Dependence of average full energy on inverted temperature for potential when: (1, 2), (3, 4), (5, 6). Here is Bohr radius, ―Hartry energy. Symbols refer to harmonic approximation calculations, curves―to usual PIMC method. Right plot: The same for average square of energy .
Figure 9. Left plot: Dependence of average kinetic on inverted temperature for potential when (1), (2), (3). Here is Bohr radius, ―Hartry energy. Results were obtained by harmonic approximation method. Right plot: The same for average potential energy.
oped approach. The Figure 10 shows potentials (33) for and different values of parameter . Solid curve corresponds to , i.e. harmonic oscillator. Another curves correspond to values of equals to and respectively. As it follows from analysis of this figure an harmonicity can be significant.
The Figure 11 shows dependence of full energy on inverted temperature. As before reference results of usual PIMC method are represented by curves, while results of harmonic approximation are marked by turned triangles. Let us stress that results of proposed above single-momentum method are marked by straight triangles. Results of all methods agree well with each other.
Next Figure 12 presents results obtained by the single-momentum method discussed above. Momentum distributions for unharmonic oscillator with (3) calculated by single-momentum method are shown on Figure 12. When , distribution is quite Maxwellian. However when temperature is lower, distribution is much wider than Maxwell one. Results of harmonic approximation for momentum distributions are the same and are not shown here.
Figure 10. Spherically symmetrical potential for and different values of unharmonic parameter: (1), (2), (3).
Figure 11. Dependence of average energy on inverted temperature for potential . Parameter for all cases, (1), (2), (3). Curves refer to usual PIMC method, up-triangles― to harmonic approximation method, down-triangles―to single momentum method.
Figure 12. Momentum distribution for spherically symmetrical potential with on different values of temperature. Dotted yellow line corresponds to Maxwell distribution.
Ideal Fermi gas Finally we consider application of single-momentum approach to system of non-interacting degenerated spinless fermions, i.e. ideal Fermi gas with Fermi-Dirac momentum distribution  :
where is energy, is chemical potential, is statistical weight equal to unity in spinless case. Chemical potential connects implicitly density and temperature by integral equation  , which can be easily numerically solved:
Thermodynamical state of ideal Fermi gas is defined by single dimensionless parameter , where is density and is the thermal wavelength. Degree of degeneracy is reflected by this parameter called as degeneracy parameter. When thermal wavelength is much lesser than average distance between particles, so gas is far from degeneracy. In opposite case degeneracy and exchange effects are significant. We have carried out our calculations for degeneracy parameter from 1 to 7.
Left plot of Figure 13 shows dependencies of spherically symmetric single-momentum density matrix (20) on Fourier variable . When density of Fermi gas is low, density matrix has form of Gaussian exponent, i.e. Fourier preimage of Maxwell distribution. When gas is degenerate, density matrix becomes oscillating function of with weak decaying. Momentum distribution functions shown on the right plot of Figure 13 by dashed lines can be numerically
Figure 13. Left plot: Dependence of single-momentum density matrix (20) on for different values of degeneracy: (1), (2), (3), (4), (5). Right plot: Fermi momentum distributions . solid lines―analytics (34), dashed lines―single-momentum numerical method, dotted line―Maxwell density matrix and distribution.
obtained through sine Fourier transform of single-momentum density matrix. Solid semi-transparent lines correspond to theoretical Fermi distribution (34). When degeneracy parameter is small then distribution is close to Maxwellian. However it tends to the “shelf” form in degenerate Fermi gas.
One can note that in cases agreement of single-momentum Monte- Carlo with explicit theoretical result (34) is very good with error about some percents. However it increase significantly for greater values of , especially at large momenta. The main source of this discrepancy is finite size of basic Monte Carlo cell. In simulations the thermal wavelength has to be much lesser than size of Monte-Carlo cell as in this case the surface effects are negligible . For example, when in our simulations number of particles in Monte-Carlo cell is about 100, for the ration and influence of boundary conditions is small. However for this ratio is more than 0.5, so influence of surface effects becomes significant. One has to increase Monte-Carlo cell and number of particles in it. Left plot of Figure 14 shows influence of boundary conditions on result for Fermi gas at . We change the number of particles in basic Monte-Carlo cell while density was constant. One can see that result for small numbers differs significantly from the others. With increasing the momentum distributions tend asymptotically to unit not depending on more. This behavior is expected as increasing at constant density leads to increasing size of Monte-Carlo cell. So to obtain better results we have to increase number of particles in Monte Carlo cell and, as a consequence, to increase the available computing power or improve algorithm of calculations. This work is now in progress.
However this limitation affects mainly on momentum distribution at high momentum values, but less essential for calculation of integrals such, for example, as average full energy , which is shown on right plot of Figure 14. Results of single momentum method agree well with analytic approaches.
Figure 14. Left plot: Single-momentum density matrix: influence of boundary conditions on result for Fermi gas with under constant density. is number of particle in Monte-Carlo cell. Right plot: Energy of ideal Fermi gas with different values of . Up-triangles corresponds to single-momentum calculation, down-triangles ―to analytical formula 34.
In this paper we have generalized path integral representation of Wigner function for canonical ensemble on many-particle non-relativistic quantum systems. Using this representation we have obtained explicit expression of Wigner function in harmonic approximation resembling the Maxwell distribution on momentum variable but with corrections depending on the second derivatives of potential field. We also propose new single-momentum approach based on Wigner function integrated over momenta of some particles.
We have developed quantum Monte-Carlo method based on harmonic approximation for Wigner function for calculating average values of arbitrary quantum operators and momentum distributions for non-ideal many-particle systems. We have tested this method on some simple systems in potential field. Obtained results have shown very good agreement with results of usual path- integral method Monte-Carlo and proved that harmonic approximation gives practically exact results even for potentials, which have no matter with harmonic potential.
Also we have developed new quantum Monte-Carlo method based on single-momentum approach. It is able to calculate momentum distributions and average values of single-particle operators without approximation. This method is suitable for non-ideal systems of fermions. We have tested it on degenerate ideal Fermi gas. Results are in good agreement with available analytical and numerical data.
We acknowledge stimulating discussions with Profs. M. Bonitz and V.I. Man’ko. This work has been supported by the Russian Science Foundation via grant 14-50-00124.
 Brown, E., Clark, B., DuBois, J. and Ceperley, D. (2013) Path-Integral Monte Carlo Simulation of the Warm Dense Homogeneous Electron Gas. Physical Review Letters, 110, Article ID: 146405.
 Groth, S., Schoof, T., Dornheim, T. and Bonitz, M. (2016) Ab Initio Quantum Monte Carlo Simulations of the Uniform Electron Gas without Fixed Nodes. Physical Review B, 93, Article ID: 085102.
 Dornheim, T., Groth, S., Schoof, T., Hann, C. and Bonitz, M. (2016) Ab Initio Quantum Monte Carlo Simulations of the Uniform Electron Gas without Fixed Nodes II: Unpolarized Case. Physical Review B, 93, Article ID: 205134.
 Filinov, V., Fehske, H., Bonitz, M., Fortov, V. and Levashov, P. (2007) Correlation Effects in Partially Ionized Mass Asymmetric Electron-Hole Plasmas. Physical Review E, 75, Article ID: 036401.
 Dornheim, T., Filinov, A. and Bonitz, M. (2015) Superfluidity of Trapped Quantum Systems in Two and Three Dimensions. Physical Review B, 91, Article ID: 054503.
 Dornheim, T., Groth, D., Filinov, A. and Bonitz, M. (2015) Permutation Blocking Path Integral Monte Carlo: A Highly Efficient Approach to the Simulation of Strongly Degenerate Non-Ideal Fermions. New Journal of Physics, 17, Article ID: 073017.
 Filinov, V., Ivanov, Y.B., Bonitz, M., Fortov, V. and Levashov, P. (2012) Color Path-Integral Monte Carlo Simulations of Quark-Gluon Plasma. Physics Letters A, 376, 1096-1101.
 Filinov, V., Ivanov, Y.B., Bonitz, M., Fortov, V. and Levashov, P. (2013) Color Path-Integral Monte-Carlo Simulations of Quark-Gluon Plasma: Thermodynamic and Transport Properties. Physical Review C, 87, Article ID: 035207.
 Filinov, V., Bonitz, M., Ivanov, Y.B., Ilgenfritz, E.-M. and Fortov, V. (2015) Color Path Integral Equation of State of the Quark-Gluon Plasma at Nonzero Chemical Potential. Plasma Physics and Controlled Fusion, 57, Article ID: 0440041.
 Eletskii, A.V., Starostin, A.N. and Taran, M.D. (2005) Quantum Corrections to the Equilibrium Rate Constants of Inelastic Processes. Physics-Uspekhi, 48, 281-294.
 Schoof, T., Groth, S. and Bonitz, M. (2014) Introduction to Configuration Path Integral Monte Carlo. In: Michael Bonitz, K., Lopez, B.J. and Thomsen, H., Eds., Introduction to Complex Plasmas: Scientific Challenges and Technological Opportunities, Springer Series on Atomic, Optical, and Plasma Physics 82, Springer, Berlin 153-194.
 Shifren, L. and Ferry, D. (2002) A Wigner Function Based Ensemble Monte Carlo Approach for Accurate Incorporation of Quantum Effects in Device Simulation. Journal of Computational Electronics, 1, 55-58.
 Sellier, J.M., Nedjalkov, M., Dimov, I. and Selberherr, S. (2014) A Benchmark Study of the Wigner Monte-Carlo Method. Monte Carlo Methods and Applications, 20, 43-51.
 Sellier, J.M. and Dimov, I. (2015) On the Simulation of Indistinguishable Fermions in the Many-Body Wigner Formalism. Journal of Computational Physics, 280, 287-294.
 Larkin, A., Filinov, V. and Fortov, V. (2016) Path Integral Representation of the Wigner Function in Canonical Ensemble. Contributions to Plasma Physics, 56, 151-360.
 Zamalin, V. and Norman, G. (1973) The Monte-Carlo Method in Feynman’s Formulation of Quantum Statistics. USSR Computational Mathematics and Mathematical Physics, 13, 408-420.