ABSTRACT

Quantum interference and exchange statistical effects can affect the momentum distribution functions making them non-Maxwellian. Such effects may be important in studies of kinetic properties of matter at low temperatures and under extreme conditions. In this work we have generalized the path integral representation for Wigner function to strongly coupled three-dimensional quantum system of particles with Boltzmann and Fermi statistics. In suggested approach the explicit expression for Wigner function was obtained in harmonic approximation and Monte Carlo method allowing numerical calculation of Wigner function, distribution functions and average quantum values has been developed. As alternative more accurate single-momentum approach and related Monte Carlo method have been developed to calculation of the distribution functions of degenerate system of interacting fermions. It allows partially overcoming the well-known sign problem for degenerate Fermi systems.

Keywords:

Wigner Function, Path Integral, Quantum Monte Carlo, Momentum Distribution

1. Introduction

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 [1] . 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 [2] , electron gas in metal [3] [4] [5] , electron-hole plasma [6] in semiconductors, superfluidity [7] [8] [9] and even quark-gluon plasma [10] [11] [12] .

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 [13] 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 [4] [5] [14] 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 [15] [16] [17] [18] [19] . Recently a generalization of these one-body Wigner Monte Carlo methods to the many-body distinguishable particles has been done in [20] . Wigner dynamics of relativistic particle was studied in [21] . 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 [22] . 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 $\hat{A}$ can be written as Weyl’s symbol $A\left(p,q\right)$ averaged over phase space with Wigner function $W\left(p,q;\beta ,V\right)$ [23] :

$\langle \hat{A}\rangle ={\displaystyle \int }\frac{{\text{d}}^{3N}p{\text{d}}^{3N}q}{{\left(2\text{π}\hslash \right)}^{3N}}A\left(p,q\right)W\left(p,q;\beta ,V\right),$ (1)

where the Weyl symbol of operator $\hat{A}$ is :

$A\left(p,q\right)={\displaystyle \int }\frac{{\text{d}}^{3N}s}{{\left(2\text{π}\hslash \right)}^{3N}}{e}^{isq/\hslash }\langle p+s/2\left|\hat{A}\right|p-s/2\rangle .$ (2)

Weyl symbols for common operators like $\hat{p}$ , $\hat{q}$ , ${\hat{p}}^2$ , ${\hat{q}}^2$ , $\hat{H}$ , ${\hat{H}}^2$ 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:

$W\left(p,q;\beta ,V\right)=Z{\left(\beta ,V\right)}^{-1}{\displaystyle \int }{\text{d}}^{3N}\xi e^{ip\xi /\hslash }\rho \left(q-\xi /2,q+\xi /2;\beta \right),$ (3)

where $Z\left(\beta ,V\right)$ is partition function of canonical ensemble of N particles, $\beta =1/kT$ is inverted temperature, V―volume, $p=\left(p_1\mathrm{,}{p}_2\mathrm{,}\cdots \mathrm{,}{p}_N\right)$ , $q=\left(q_1\mathrm{,}{q}_2\mathrm{,}\cdots \mathrm{,}{q}_N\right)$ , $pq={\displaystyle {\sum }_{a=1}^N}{p}_a{q}_a$ etc. Hamiltonian of the system consists of kinetic and potential parts: $\hat{H}=\hat{K}+\hat{U}$ , where $\hat{K}={\displaystyle {\sum }_{a=1}^N}{\hat{p}}^2/2m_a$ , $\hat{U}=U\left({\hat{q}}_1\mathrm{,}{\hat{q}}_2\mathrm{,}\cdots \mathrm{,}{\hat{q}}_N\right)$ ―potential energy of particles, interacting with each other and/or external field. For simplicity we will consider system containing identical particles of mass $m$ . Generalization on systems with several kinds of particles is direct. The Wigner function $W\left(p,q;t\right)$ 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

$\rho \left(q-\xi /2,q+\xi /2;\beta \right)={\displaystyle \underset{P}{\sum }}{\left(\pm 1\right)}^P\langle q-\xi /2\left|e^{-\beta \hat{H}}\right|P\left(q+\xi /2\right)\rangle ,$ (4)

where symbol $P$ denotes particle permutations, $+1$ corresponds to Bose statistics, $-1$ ―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 [1] [24] . So we represent statistical operator $e^{-\beta \hat{H}}$ as product of large number $\left(M\right)$ of operators $e^{-ϵ\hat{H}}$ formally related to high temperature $\left(ϵ=\beta /M\right)$ and write down Wigner function in form of multiple integral of the product of the high temperature density matrices:

$\begin{array}{c}W\left(p,q;\beta ,V\right)=Z{\left(\beta ,V\right)}^{-1}{\displaystyle \int }{\text{d}}^{3N}\xi e^{ip\xi /\hslash }{\displaystyle \int }{\text{d}}^{3N}{q}^1\cdots {\text{d}}^{3N}{q}^{M-1}\\ \times {\left[{\displaystyle \underset{m=0}{\overset{M-1}{\prod }}}\langle q^m\left|e^{-ϵ\hat{H}}\right|q^{m+1}\rangle \right]|}_{q^0=q-\xi /2}^{q^M=P\left(q+\xi /2\right)}.\end{array}$ (5)

In the limit $ϵ\to 0$ the high temperature density matrices can be easily calculated with accuracy up to $O\left({ϵ}^2\right)$ [1] :

$\langle q^m\left|e^{-ϵ\hat{H}}\right|q^{m+1}\rangle \approx {\lambda }_{ϵ}^{-1}\exp \left[-\text{π}{\lambda }_{ϵ}^{-2}{\left(q^{m+1}-q^m\right)}^2-ϵU\left(q^m\right)\right],$ (6)

where we use abbreviation ${\left(q^{m+1}-q^m\right)}^2={\displaystyle {\sum }_{a=1}^N}{\left(q_a^{m+1}-q_a^m\right)}^2$ , and ${\lambda }_{ϵ}=\sqrt{2\text{π}{\hslash }^2ϵ/m}$ is the thermal wavelength at high temperature $M\cdot T$ . Thus expression for Wigner function takes the form:

$\begin{array}{c}W\left(p,q;\beta ,V\right)=Z{\left(\beta ,V\right)}^{-1}{\displaystyle \int }{\text{d}}^{3N}\xi e^{ip\xi /\hslash }\times {\displaystyle \underset{P}{\sum }}{\left(\pm 1\right)}^P{\displaystyle \int }{\text{d}}^{3N}{q}^1\cdots {\text{d}}^{3N}{q}^{M-1}\\ \times \exp {\left\{-\frac{1}{M}{\displaystyle \underset{m=0}{\overset{M-1}{\sum }}}\left[\text{π}\frac{M^2}{{\lambda }^2}{\left(q^{m+1}-q^m\right)}^2+\beta U\left(q^m\right)\right]\right\}|}_{q_0=\left(q-\xi /2\right)}^{q_M=P\left(q+\xi /2\right)},\end{array}$ (7)

where $\lambda =\sqrt{2\text{π}{\hslash }^2\beta /m}$ is the thermal wavelength at temperature $T$ . In continuous limit $M\to \infty$ this expression turns into path integral over trajectories $q\left(\tau \right)$ starting at point $\left(q-\xi /2\right)$ and ending at $P\left(q+\xi /2\right)$ , where $\tau$ is dimensionless parameter, which may be interpreted as “imaginary time” [1] :

$\begin{array}{c}W\left(p,q;\beta ,V\right)=Z{\left(\beta ,V\right)}^{-1}{\displaystyle \int }{\text{d}}^{3N}\xi e^{ip\xi /\hslash }{\displaystyle {\sum }_P{\left(\pm 1\right)}^P}{\displaystyle {\int }_{q\left(0\right)=\left(q-\xi /2\right)}^{q\left(1\right)=P\left(q+\xi /2\right)}}{D}^{3N}q\left(\tau \right)\\ \times \exp \left\{-\text{π}{\displaystyle {\int }_0^1}\text{d}\tau \left[\frac{{\stackrel{\dot{}}{q}}^2\left(\tau \right)}{{\lambda }^2}+\frac{\beta }{\text{π}}U\left(q\left(\tau \right)\right)\right]\right\}.\end{array}$ (8)

Measure of the path integral depends on Fourier variables $\xi$ and positions $q$ . To get rid this dependence we change variables from trajectories $q\left(\tau \right)$ to closed dimensionless trajectories $z\left(\tau \right)$ :

$q\left(\tau \right)=\lambda z\left(\tau \right)+\left(1-\tau \right)\left(q-\xi /2\right)+\tau P\left(q+\xi /2\right).$ (9)

Taking into account new boundary conditions $z\left(0\right)=z\left(1\right)=0$ we obtain the desired path integral representation of Wigner function:

$\begin{array}{c}W\left(p,q;\beta ,V\right)=Z{\left(\beta ,V\right)}^{-1}{\displaystyle \int }{\text{d}}^{3N}\xi e^{ip\xi /\hslash }{\displaystyle {\sum }_P{\left(\pm 1\right)}^P}{e}^{-\frac{\text{π}}{{\lambda }^2}{\left[P\left(q+\xi /2\right)-\left(q-\xi /2\right)\right]}^2}\\ \times {\displaystyle {\int }_{z\left(0\right)=z\left(1\right)=0}}{D}^{3N}z\left(\tau \right)e^{-\text{π}{\displaystyle {\int }_0^1}\text{d}\tau \left[{\stackrel{\dot{}}{z}}^2\left(\tau \right)+\frac{\beta }{\text{π}}U\left(q\left(\tau \right)\right)\right]},\end{array}$ (10)

where $q\left(\tau \right)$ is given by (9), symbol $P$ denotes permutation of identical particles and $Dz\left(\tau \right)$ means integration over trajectories. In fact, this expression means that Wigner function is Fourier transform over $\xi$ 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 $p$ in expression (10) for Wigner function are connected with other variables through $3N$ -dimensional Fourier transform, which is not integrable analytically or numerically in general case. Exclusions are the linear or harmonic potentials, when power of variable $\xi$ is not more than two. So let us take the approximation for potential $U\left(q\right)$ 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 $\xi$ :

$\begin{array}{c}U\left(\lambda z\left(\tau \right)+q+\xi \left(\tau -1/2\right)\right)\approx U\left(\lambda z\left(\tau \right)+q\right)+\left(\tau -1/2\right){\xi }_{a,i}\frac{\partial U\left(\lambda z\left(\tau \right)+q\right)}{\partial q_{a,i}}\\ +\frac{1}{2}{\left(\tau -1/2\right)}^2{\xi }_{a,i}{\xi }_{b,j}\frac{{\partial }^{2}U\left(\lambda z\left(\tau \right)+q\right)}{\partial q_{a,i}\partial q_{b,j}},\end{array}$ (11)

where we denote summation over repeated indices $a\mathrm{,}b$ from 1 to $N$ and $i\mathrm{,}j$ over from 1 to 3. If we take expansion (11), the expression for Wigner function (10) takes form of generalized Gaussian integral over variable $\xi$ and expression for Wigner function in harmonic approximation can be written in the following form:

$\begin{array}{c}W\left(p,q;\beta ,V\right)\\ =Z{\left(\beta ,V\right)}^{-1}{\displaystyle {\int }_{z\left(0\right)=z\left(1\right)=0}}{D}^{3N}z\left(\tau \right)\exp \left\{-\text{π}K\left[z\left(\tau \right)\right]-\text{π}P\left[q,z\left(\tau \right)\right]\right\}\\ \times \exp \left\{-\frac{{\lambda }^2}{4\text{π}{\hslash }^2}{p}_{a,i}{\chi }_{ai,bj}\left[q,z\left(\tau \right)\right]p_{b,j}+\text{π}{J}_{a,i}\left[q,z\left(\tau \right)\right]{\chi }_{ai,bj}\left[q,z\left(\tau \right)\right]J_{b,j}\left[q,z\left(\tau \right)\right]\right\}\\ \times \det {\left|{\chi }_{ai,bj}\left[q,z\left(\tau \right)\right]\right|}^{1/2}\cos \left\{\frac{\lambda }{\hslash }{p}_{a,i}{\chi }_{ai,bj}\left[q,z\left(\tau \right)\right]J_{b,j}\left[q,z\left(\tau \right)\right]\right\}.\end{array}$ (12)

Here we introduced scalar functionals $K$ and $P$ , as well as functionals $J_{a\mathrm{,}i}$ and ${\chi }_{ai\mathrm{,}bj}$ of first and second potential derivatives having form of $3N$ - dimensional vector and matrix correspondingly. They depends on trajectories $z\left(\tau \right)$ and positions $q$ :

$K\left[z\left(\tau \right)\right]={\displaystyle {\int }_0^1}\text{d}\tau {\stackrel{\dot{}}{z}}^2\left(\tau \right),$

$P\left[q,z\left(\tau \right)\right]=\frac{\beta }{\text{π}}{\displaystyle {\int }_0^1}\text{d}\tau U\left(q+\lambda z\left(\tau \right)\right),$

$J_{a,i}\left[q,z\left(\tau \right)\right]=\frac{\beta \lambda }{2\text{π}}{\displaystyle {\int }_0^1}\text{d}\tau \left(\tau -1/2\right)\frac{\partial U}{\partial q_{a,i}}\left(q+\lambda z\left(\tau \right)\right),$

${\chi }_{ai,bj}\left[q,z\left(\tau \right)\right]={\left[{\delta }_{ab}{\delta }_{ij}+\frac{\beta {\lambda }^2}{2\text{π}}{\displaystyle {\int }_0^{\beta }}\text{d}\tau {\left(\tau -1/2\right)}^2\frac{{\partial }^{2}U}{\partial q_{a,i}\partial q_{b,j}}\right]}^{-1}.$

Note that the first term in exponent (12) looks similarly to Maxwell distribution in classical statistics. The major difference is in matrix ${\chi }_{ai,bj}\left[q,z\left(\tau \right)\right]$ 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 $U$ is expandable in Taylor series of second order $\xi$ 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 $\xi$ in three positions. The first one is in Fourier term $\left(ip\xi /\hslash \right)$ , which makes momenta correlated with other dynamical variables. The second one is in gaussian-like term $\left(\text{π}{\xi }^2/{\lambda }^2\right)$ . The third one is in integral term, where $\xi$ is argument of potential:

${\displaystyle \underset{0}{\overset{1}{\int }}}\text{d}\tau \beta U\left(q+\lambda z\left(\tau \right)+\xi \left(\tau -1/2\right)\right)$ . Gaussian term provides rapid decaying when $\xi$

increases, so main contribution comes from $\xi$ near ${\text{π}}^{-1/2}\approx 0.6$ . Argument of potential function contains $\xi$ multiplied on $\tau -1/2$ , which is modulo less than 0.5. Using mean value theorem, we can roughly get symbolic estimation of these integrals in exponent (10) as

$\frac{1}{n!}\frac{{\partial }^{n}U\left(q_0\right)}{\partial q^n}{\displaystyle \underset{0}{\overset{1}{\int }}}\text{d}\tau {\left(\tau -\frac{1}{2}\right)}^n=\frac{1}{n!}\frac{\left(1+{\left(-1\right)}^n\right)}{\left(n+1\right)2^{n+1}}\frac{{\partial }^{n}U\left(q_0\right)}{\partial q^n},$ (13)

where $x_0$ is certain point of trajectory. Numerical value of this integral rapidly decreases: when $n=2,4,6$ it equals $1/24$ , $1/1920$ and $1/322560$ . 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 $\langle \hat{A}\rangle$ 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) [3] 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:

$\langle \hat{A}\rangle =\frac{{\langle A\left(p,q\right)\cdot f\left(p,q,z^1,\cdots ,z^{M-1}\right)\cdot h\left(p,q,z^1,\cdots ,z^{M-1}\right)\rangle }_w}{{\langle f\left(p,q,z^1,\cdots ,z^{M-1}\right)\cdot h\left(p,q,z^1,\cdots ,z^{M-1}\right)\rangle }_w}.$ (14)

Here brackets ${\langle g\left(p,q,z^1,\cdots ,z^{M-1}\right)\rangle }_w$ denote averaging of any function $g\left(p,q,z^1,\cdots ,z^{M-1}\right)$ with positive weight $w\left(p,q,z^1,\cdots ,z^{M-1}\right)$ :

#Math_106# (15)

while

$\begin{array}{l}w\left(p,q,z^1,\cdots ,z^{M-1}\right)=\exp \left\{-\frac{{\lambda }^2}{4\text{π}{\hslash }^2}{p}_{a,i}{\chi }_{ai,bj}{p}_{b,j}+\text{π}{J}_{a,i}{\chi }_{ai,bj}{J}_{b,j}-\text{π}K-\text{π}P\right\},\\ f\left(p,q,z^1,\cdots ,z^{M-1}\right)=\cos \left\{\frac{\lambda }{\hslash }{p}_{a,i}{\chi }_{ai,bj}{J}_{b,j}\right\},\\ h\left(p,q,z^1,\cdots ,z^{M-1}\right)=\det {\left|{\chi }_{ai,bj}\right|}^{1/2}.\end{array}$ (16)

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 $3N$ -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 $A$ , whose Weyl symbol is a sum of equal single-particle symbols: $A\left(p,q\right)={\displaystyle {\sum }_{a=1}^N}a\left(p_a\mathrm{,}{q}_a\right)$ . Due to the fact that all particles are identical, average values on ensemble of $\hat{a}$ for arbitrary particles $a$ and $b$ are equal: $\langle {\hat{a}}_a\rangle =\langle {\hat{a}}_b\rangle$ . Therefore we can use Weyl symbol $N\cdot a\left(\left(p_1\right),\left(q_1\right)\right)$ instead of ${\displaystyle {\sum }_{a=1}^N}a\left(p_a\mathrm{,}{q}_a\right)$ and integrate Wigner function over other momenta. Thus in single-momentum formalism:

$\langle \hat{A}\rangle ={\displaystyle \int }\frac{{\text{d}}^3{p}_1}{{\left(2\text{π}\hslash \right)}^3}{\text{d}}^3{q}_1\cdots {\text{d}}^3{q}_{N}N\cdot A\left(p_1\right)W_{sm}\left(p_1\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)\mathrm{.}$ (17)

Here we introduce single-momentum Wigner function as integral of $W\left(p,q;\beta ,V\right)$ over all momenta except $p_1$ , which we will denote by $p$ further:

$W_{sm}\left(p\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)={\displaystyle \int }{\text{d}}^3{p}_2\cdots {\text{d}}^3{p}_{N}W\left(p_1\mathrm{,}\cdots \mathrm{,}{p}_N\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)\mathrm{.}$ (18)

As momenta in definition of Wigner function appear only in Fourier transform $e^{i{p}_a{\xi }_a/\hslash }$ , integrals on them gives delta functions $\delta \left({\xi }_a\right)$ for $a=2,\cdots ,N$ . This exterminates all ${\xi }_a$ except ${\xi }_1$ , which now is denoted as simple $\xi$ . Thus formula for single-momentum Wigner function is

$W_{sm}\left(p\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)={\displaystyle \int }{\text{d}}^3\xi e^{ip\xi /\hslash }\rho \left(\xi \mathrm{<mo>\; ;\; </mo>}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)\mathrm{,}$ (19)

i.e. it is 3-dimensional Fourier transform of single-momentum density matrix:

$\begin{array}{c}\rho \left(\xi \mathrm{<mo>\; ;\; </mo>}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)=Z{\left(\beta \mathrm{,}V\right)}^{-1}D\left[\xi \mathrm{,}q\right]\\ \times {\displaystyle {\int }_{z\left(0\right)=z\left(1\right)=0}}{D}^{3N}z\left(\tau \right)\exp \left\{-\text{π}K\left[z\left(\tau \right)\right]-\text{π}P\left[\xi \mathrm{,}q\mathrm{,}z\left(\tau \right)\right]\right\}\mathrm{.}\end{array}$ (20)

Here $K$ is the same functional on trajectories as in (13), P and D are potential functional and determinant of permutations:

$\begin{array}{l}P\left[\xi ,q,z\left(\tau \right)\right]=\frac{\beta }{\text{π}}{\displaystyle {\int }_0^{\beta }}\text{d}\tau U\left(q+\lambda z\left(\tau \right)+\xi \left(\tau -1/2\right)\right),\\ D\left[\xi ,q\right]=\det \left|-\frac{\text{π}}{{\lambda }^2}{\left({\left(q+\xi /2\right)}_a-{\left(q-\xi /2\right)}_b\right)}^2\right|.\end{array}$ (21)

Such reasoning would not change if we calculate average value of operators of the form $\hat{C}={\displaystyle {\sum }_{a=1}^N}c\left(p_a\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)$ or $\hat{B}=B\left(q_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)$ . To summarize, in single-momentum approach we are able to deal with quantum operators of particular form only:

$\hat{A}={\displaystyle \underset{a=1}{\overset{N}{\sum }}}c\left({\hat{p}}_a\mathrm{,}{\hat{q}}_1\mathrm{,}\cdots \mathrm{,}{\hat{q}}_N\right)+B\left({\hat{q}}_1\mathrm{,}\cdots \mathrm{,}{\hat{q}}_N\right)\mathrm{,}$ (22)

so average value of it is:

#Math_135# (23)

with Weyl symbol depending only on first particle’s momentum $p:A\left(p\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)=N\cdot c\left(p\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)+B\left(q_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)$ . 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 $\lambda$ and trajectories in path integrals (see (10), (9)) are strongly entangled. This is the reason that permutations can not strongly affect the potential energy $\beta U\left(\lambda z\left(\tau \right)+\left(1-\tau \right)\left(q-\xi /2\right)+\tau P\left(q+\xi /2\right)\right)$ in comparison with the case of identical permutation.

So $\beta U\left(\lambda z\left(\tau \right)+\left(1-\tau \right)\left(q-\xi /2\right)+\tau P\left(q+\xi /2\right)\right)$ can be presented as energy related to identical permutation $\beta U\left(\lambda z\left(\tau \right)+q+\xi \left(\tau -1/2\right)\right)$ and small difference $\beta U\left(\lambda z\left(\tau \right)+\left(1-\tau \right)\left(q-\xi /2\right)+\tau P\left(q+\xi /2\right)\right)-\beta U\left(\lambda z\left(\tau \right)+q+\xi \left(\tau -1/2\right)\right)$ , which, for example, for 1D Coulomb system is of order

$\frac{\beta \partial U\left(\lambda \left(n\lambda \right)\right)}{\partial q}\lambda \left(n\lambda \right)\approx \sqrt{\frac{\beta Ry}{\text{π}}}/\left(n\lambda \right)$ , where $n$ is 1D density of particles. Then

for the first terms of the perturbation series on this parameter one can obtain the following simplification: $U\left(\lambda z\left(\tau \right)+\left(1-\tau \right)\left(q-\xi /2\right)+\tau P\left(q+\xi /2\right)\right)\approx U\left(\lambda z\left(\tau \right)+q+\xi \left(\tau -1/2\right)\right)$ .

Now all permutations in (10) are connected with variables $q$ and $\xi$ and can be beared out of the path integral $D^{3N}z\left(\tau \right)$ . Moreover, one can gather all permutations into Slater determinant of $N\times N$ matrix:

$\begin{array}{l}W\left(p,q;\beta ,V\right)=Z{\left(\beta ,V\right)}^{-1}{\displaystyle \int }{\text{d}}^{3N}\xi e^{ip\xi /\hslash }\det \left|-\frac{\text{π}}{{\lambda }^2}{\left({\left(q+\xi /2\right)}_a-{\left(q-\xi /2\right)}_b\right)}^2\right|\\ \times {\displaystyle {\int }_{z\left(0\right)=z\left(1\right)=0}}{D}^{3N}z\left(\tau \right)e^{-\text{π}{\displaystyle {\int }_0^1}\text{d}\tau \left[{\stackrel{\dot{}}{z}}^2\left(\tau \right)+\frac{\beta }{\text{π}}U\left(\lambda z\left(\tau \right)+q+\xi \left(\tau -1/2\right)\right)\right]}.\end{array}$ (24)

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. $\rho \left(\xi \mathrm{<mo>\; ;\; </mo>}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)$ in form of multiple integral over “intermediate quantum coordinates” $z^1,\cdots ,z^{M-1}$ . Firstly, by standard Monte-Carlo procedure we are able to calculate Fourier preimage of averaged single-momentum operator $\hat{A}$ of form (22):

$\langle \hat{A}\rangle \left(p\mathrm{,}\xi \right)=\frac{{\langle A\left(p\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)\cdot f_{sm}\left(\xi \mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)\rangle }_{w_{sm}}}{{\langle f_{sm}\left(\xi \mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)\rangle }_{w_{sm}}}\mathrm{,}$ (25)

which depends on momentum $p$ by Weyl symbol and $\xi$ only. Here brackets ${\langle g\left(p\mathrm{,}\xi \mathrm{,}q\mathrm{,}{z}^1\mathrm{,}\cdots \mathrm{,}{z}^{M-1}\right)\rangle }_{w_{sm}}$ denote averaging of any function $g$ with positive weight $w_{sm}\left(\xi ,q,z^1,\cdots ,z^{M-1}\right)$ :

#Math_159# (26)

while

$\begin{array}{l}{w}_{sm}\left(\xi ,q,z^1,\cdots ,z^{M-1}\right)=\\ =\left|D\left[\xi ,q\right]\right|{\displaystyle {\int }_{z\left(0\right)=z\left(1\right)=0}}{D}^{3N}z\left(\tau \right)exp\left\{-\text{π}K\left[z\left(\tau \right)\right]-\text{π}P\left[\xi ,q,z\left(\tau \right)\right]\right\}\mathrm{,}\\ f_{sm}\left(\xi \mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\right)=\text{sign}D\left[\xi ,q\right].\end{array}$ (27)

So we obtain $\langle \hat{A}\rangle \left(p\mathrm{,}\xi \right)$ in form of $\xi$ distribution function. Then we can do Fourier transform over variable $\xi$ numerically, using Fast Fourier Transform algorithms. After that it is easy to average expression over momenta to obtain desired average value $\langle \hat{A}\rangle$ .

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 $\xi$ without Fourier transform. The terms which are functions of $p$ can be calculated from single-momentum distribution function:

$W_p\left(p\right)={\displaystyle \int }{\text{d}}^{3N}q{W}_{sm}\left(p\mathrm{,}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)\mathrm{,}$ (28)

obtained by Fourier transform of single-momentum density matrix $\rho \left(\xi \mathrm{<mo>\; ;\; </mo>}{q}_1\mathrm{,}\cdots \mathrm{,}{q}_N\mathrm{<mo>\; ;\; </mo>}\beta \mathrm{,}V\right)$ integrated over all positions $q$ .

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 $1D$ and $3D$ 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 $1D$ -power potential of the fourth order, i.e. unharmonic oscillator:

$V_{2-4}\left(x\right)=a{x}^2+c{x}^4,$ (29)

where $a$ and $c$ are adjusted parameters of potential. When $c=0$ this potential describes harmonic oscillator and formula (12) is exact. If $c$ does not equal zero, oscillator is unharmonic and can be used for verification of the harmonic approximation. The Figure 1 shows potentials (29) for $a=0.5m{\omega }^2$ and different values of parameter $c$ when anharmonicity is significant.

Left plot of Figure 2 shows dependence of full energy on inverted temperature for potential $V_{2-4}\left(x\right)$ . 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 $\beta \hslash \omega =1.5$ to $\beta \hslash \omega =5$ ) the energy tends to constant value, namely to the ground state energy $E_0$ .

Right plot of Figure 2 shows dependence of squared energy on inverted temperature for potential $V_{2-4}\left(x\right)$ . From analysis of this figure it follows that harmonic approximation gives correct enough values of $\langle H^2\rangle$ even for the ground state. Weyl’s symbol of ${\hat{H}}^2$ is polynomial of the 8 order, so the main contribution to $\langle H^2\rangle$ 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-

ues as it calculates derivatives of partition function $Z\left(\beta \right)$ .

As a second example we consider $1D$ -potential of the fourth order:

$V_{3-4}\left(x\right)=E\left(b{\left(x/\sigma \right)}^3+{\left(x/\sigma \right)}^4\right),$ (30)

where $E={\hslash }^2/m{\sigma }^2$ is characteristic energy of potential. Dimensionless parameter $b$ specifies contribution of the cubic term. Figure 4 shows potentials (30) for different values of this parameter. With gain of $b$ 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 $b$ 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 $\langle {\hat{H}}^2\rangle$ 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 $b$ on kinetic energy, while the potential energy rapidly decrease when $b$ 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 $b$ when expression below for the most trajectories $z\left(\tau \right)$ becomes negative:

$1+\frac{\beta {\lambda }^2}{\text{π}}{{\displaystyle \underset{0}{\overset{1}{\int }}}\text{d}\tau {\left(\tau -\frac{1}{2}\right)}^2\frac{{\partial }^{2}V\left(y\right)}{\partial y^2}|}_{y=x+\lambda z\left(\tau \right)}\le 0,$ (31)

and functional $\chi \left[x,z\right]$ 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):

$V_{\text{scc}}\left(x\right)=-\frac{e^2}{\sqrt{x_0^2+x^2}}.$ (32)

Here $e$ is electrical charge, while parameter $x_0$ is characteristic length of potential. Sometimes the SCC potential is used in numerical calculations instead of true Coulomb potential to avoid Coulomb divergency [25] . Figure 7 shows SCC potential for different values of $x_0$ . Curves refer to $x_0=a_0$ , $x_0=2a_0$ and $x_0=3a_0$ , where $a_0=\hslash /m{e}^2$ is Bohr radius. It coincides with Bohr radius in three dimension [25] . 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 $\langle H^2\rangle$ represented on right plot of Figure 8. At high-temperatures functions $\langle \hat{H}\rangle$ and $\langle H^2\rangle$ are not affected by parameter $x_0$ , while at low- temperature the strong dependence on $x_0$ 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 $x_0$ .

It should be noted, that some points are missed when $\beta Ha\ge 7$ . The reason is the same as for potential $V_{3-4}$ . This happens at very low temperature, when system is already in ground state. However when $x_0\le 0.8a_0$ harmonic approximation fails to calculate ground state, as the second derivative of potential $V_{\text{scc}}$ 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:

$V_{2-4}\left(r\right)=a{r}^2+c{r}^4,$ (33)

where $a$ and $c$ are adjusted parameters of potential, $r$ is length of radius-vector $\left(x\mathrm{,}y\mathrm{,}z\right)$ . When $c=0$ this potential describes symmetrical $3D$ - harmonic oscillator and harmonic approximation (12) is exact. If $c$ does not equal zero, oscillator is unharmonic and can be used for verification of the devel-

oped approach. The Figure 10 shows potentials (33) for $a=0.5m{\omega }^2$ and different values of parameter $c$ . Solid curve corresponds to $c=0$ , i.e. harmonic oscillator. Another curves correspond to values of $c$ equals to $m^2{\omega }^3/2\hslash$ and $2m^2{\omega }^3/\hslash$ 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 $c=m^2{\omega }^3/2\hslash$ (3) calculated by single-momentum method are shown on Figure 12. When $\beta Ha\le 1$ , 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.

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 [26] :

$n\left(p\right)=\frac{g}{2{\text{π}}^2{\hslash }^3}\frac{1}{e^{\left(ϵ-\mu \right)/kT}+1},$ (34)

where $ϵ=p^2/2m$ is energy, $\mu$ is chemical potential, $g=2s+1$ is statistical weight equal to unity in spinless case. Chemical potential connects implicitly density and temperature by integral equation [27] , which can be easily numerically solved:

$\frac{N}{V}=\frac{g{\left(mkT\right)}^{3/2}}{\sqrt{2}{\text{π}}^2{\hslash }^3}{\displaystyle {\int }_0^{\infty }}\frac{\sqrt{z}\text{d}z}{e^{z-\mu /kT}+1}.$ (35)

Thermodynamical state of ideal Fermi gas is defined by single dimensionless parameter $\chi =n{\lambda }^3$ , where $n$ is density and $\lambda$ is the thermal wavelength. Degree of degeneracy is reflected by this parameter called as degeneracy parameter. When $\chi \ll 1$ 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 $\chi$ from 1 to 7.

Left plot of Figure 13 shows dependencies of spherically symmetric single-momentum density matrix (20) $\rho$ on Fourier variable $\xi$ . 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 $\xi$ with weak decaying. Momentum distribution functions shown on the right plot of Figure 13 by dashed lines can be numerically

obtained through sine Fourier transform of single-momentum density matrix. Solid semi-transparent lines correspond to theoretical Fermi distribution (34). When degeneracy parameter $\chi$ 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 $\chi \le 5$ 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 $\chi$ , 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 $\left(\lambda \ll L\right)$ . For example, when in our simulations number of particles in Monte-Carlo cell is about 100, for $\chi =1.3$ the ration $\lambda /L\approx 0.25$ and influence of boundary conditions is small. However for $\chi =10$ 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 $\chi =3$ . We change the number of particles $N$ in basic Monte-Carlo cell while density was constant. One can see that result for small numbers $N$ differs significantly from the others. With increasing $N$ the momentum distributions tend asymptotically to unit not depending on $N$ more. This behavior is expected as increasing $N$ 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 $\langle \hat{H}\rangle$ , which is shown on right plot of Figure 14. Results of single momentum method agree well with analytic approaches.

6. Conclusions

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.

Acknowledgements

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.

Cite this paper

Larkin, A.S. and Filinov, V.S. (2017) Phase Space Path Integral Representation for Wigner Function. Journal of Applied Mathematics and Physics, 5, 392-411. https://doi.org/10.4236/jamp.2017.52035

References