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
can be written as Weyl’s symbol
averaged over phase space with Wigner function
[23] :
(1)
where the Weyl symbol of operator
is :
(2)
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:
(3)
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
(4)
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 [1] [24] . 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:
(5)
In the limit
the high temperature density matrices can be easily calculated with accuracy up to
[1] :
(6)
where we use abbreviation
, and
is the thermal wavelength at high temperature
. Thus expression for Wigner function takes the form:
(7)
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” [1] :
(8)
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
:
(9)
Taking into account new boundary conditions
we obtain the desired path integral representation of Wigner function:
(10)
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
:
(11)
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:
(12)
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
(13)
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) [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:
(14)
Here brackets
denote averaging of any function
with positive weight
:
#Math_106# (15)
while
(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
-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:
(17)
Here we introduce single-momentum Wigner function as integral of
over all momenta except
, which we will denote by
further:
(18)
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
(19)
i.e. it is 3-dimensional Fourier transform of single-momentum density matrix:
(20)
Here
is the same functional on trajectories as in (13), P and D are potential functional and determinant of permutations:
(21)
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:
(22)
so average value of it is:
#Math_135# (23)
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:
(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.
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):
(25)
which depends on momentum
by Weyl symbol and
only. Here brackets
denote averaging of any function
with positive weight
:
#Math_159# (26)
while
(27)
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:
(28)
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:
(29)
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:
(30)
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:
(31)
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):
(32)
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 [25] . 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 [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
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:
(33)
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 [26] :
(34)
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 [27] , which can be easily numerically solved:
(35)
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.
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.