Quantum Statistical Mechanics Results for Argon, Neon, and Helium Using Classical Monte Carlo

TL;DR

This paper computes quantum corrections to classical pressure for argon, neon, and helium using Monte Carlo simulations, revealing significant quantum effects especially in helium at liquid-like densities.

Contribution

It introduces a method to calculate quantum corrections up to fourth order in Planck's constant for Lennard-Jones gases using classical Monte Carlo simulations.

Findings

Quantum corrections are about 1% for argon, 1 for neon, and 1000 for helium.

First order wave function symmetrization effects are negligible.

Quantum effects are more pronounced in lighter gases like helium.

Abstract

Quantum corrections to the classical pressure are obtained for Lennard-Jones models of argon, neon, and helium using classical Metropolis algorithm computer simulations. The corrections for non-commutativity are obtained to fourth order in Planck's constant. Compared to the classical virial pressure on all isotherms at liquid-like densities, the quantum correction is found to be ${\cal O}(10^{-2})$ for argon, ${\cal O}(10^{0})$ for neon, and ${\cal O}(10^{3})$ for helium. The first order correction due to wave function symmetrization is also obtained, but this is relatively negligible.

Tables

Table 1. Table 1: Lennard-Jones Parameters. Sciver12

	$\sigma$	$r_{\min }$	$ϵ/k_B$	$m$
	(nm)	(nm)	(K)	(amu)
He	0.2556	0.2869	10.22	4.002602
Ne	0.2789	0.3131	35.7	20.1797
Ar	0.3418	0.3837	124.0	39.948

Equations

$${\varphi }_{\mathbf{n}}^{\pm }=\frac{1}{\sqrt{N!{\chi }_{\mathbf{n}}^{\pm }}}\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p{\varphi }_{\hat{\mathrm{P}}\mathbf{n}}.$$

$${\chi }_{\mathbf{n}}^{\pm }=\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p⟨{\varphi }_{\mathbf{n}}|{\varphi }_{\hat{\mathrm{P}}\mathbf{n}}⟩.$$

$$\sum _{\mathbf{n}}{}^{\prime }{f}_{\mathbf{n}}=\frac{1}{N!}\sum _{\mathbf{n}}{\chi }_{\mathbf{n}}^{\pm }{f}_{\mathbf{n}}.$$

$${\chi }_{11}^{\pm }$$

$$=$$

$${\chi }_{12}^{\pm }$$

$$\sum _{\mathbf{n}}{}^{\prime }{f}_{\mathbf{n}}$$

$$=$$

$${\Xi }^{\pm }$$

$${\zeta }_{\mathbf{p}}(\mathbf{r})$$

$${\zeta }_{\mathbf{q}}(\mathbf{r})$$

$$⟨{\zeta }_{\mathbf{p}}|{\zeta }_{\mathbf{q}}⟩$$

$$⟨{\zeta }_{\mathbf{p}}|{\zeta }_{\mathbf{q}}⟩⟨{\zeta }_{\mathbf{q}}|{\zeta }_{\mathbf{p}}⟩\equiv \frac{(8\pi {\xi }^2{)}^{3N/2}}{V^N}{e}^{-2{\xi }^2{p}^2/{\hslash }^2},$$

$$e^{-\beta \hat{\mathcal{H}}}{\zeta }_{\mathbf{p}}(\mathbf{r})$$

$$\widehat{\tilde{\mathcal{H}}}$$

$$⟨{\zeta }_{\mathbf{q}}|e^{-\beta \hat{\mathcal{H}}}|{\zeta }_{\mathbf{p}}⟩$$

$$\sum _{\mathbf{q},\mathbf{p}}\Rightarrow \frac{1}{({\Delta }_p{\Delta }_q{)}^{3N}}\int \mathrm{d}\mathbf{\Gamma }.$$

$${\Xi }^{\pm }$$

$${\eta }^{\pm }(\mathbf{\Gamma })$$

$${\Xi }_1=\sum _{N=0}^{\infty }\frac{z^N}{N!h^{3N}}\int \mathrm{d}\mathbf{\Gamma }{e}^{-\beta \mathcal{H}(\mathbf{\Gamma })}{e}^{w(\mathbf{\Gamma })}.$$

$$\frac{{\Xi }^{\pm }}{{\Xi }_1}$$

$${\Xi }_{1,0}\equiv {\Xi }_1(w=0)=\sum _{N=0}^{\infty }\frac{z^N}{N!h^{3N}}\int \mathrm{d}\mathbf{\Gamma }{e}^{-\beta \mathcal{H}(\mathbf{\Gamma })}.$$

$$e^{-\beta [{\Omega }_1-{\Omega }_{1,0}]}=\frac{{\Xi }_1}{{\Xi }_{1,0}}={⟨e^w⟩}_{1,0},$$

$$\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p\hat{\mathrm{P}}$$

$${\eta }^{(1)}(\mathbf{\Gamma })\equiv \frac{⟨{\zeta }_{\mathbf{p}}|{\zeta }_{\mathbf{q}}⟩}{⟨{\zeta }_{\mathbf{p}}|{\zeta }_{\mathbf{q}}⟩}=1.$$

$${\eta }_{jk}^{\pm (2)}(\mathbf{\Gamma })$$

$${\eta }_{jkl}^{\pm (3)}(\mathbf{\Gamma })$$

$${\eta }^{\pm }(\mathbf{\Gamma })$$

$$\frac{{\Xi }^{\pm }}{{\Xi }_1}$$

Full text

I Introduction

The Lennard-Jones potential is arguably the simplest inter-particle potential that contains the two features of realistic atoms and molecules that are necessary for the properties of condensed matter. These are a short-range repulsion that reflects the particle size, and a long-range attraction that is responsible for the cohesion of the particles and for the existence of the liquid phase, as well as for the nature of the solid phase. That the attraction decays with the sixth power of separation is due to the correlated electronic fluctuations that induce temporary dipoles according to their polarizability, and again this is a feature of all real atoms and molecules.

Although extensive tabulations of Lennard-Jones parameters for various molecules exist,CRC the potential is most suited for spherical atoms and molecules since it is a function of the separation alone. The noble gas atoms in particular have been fitted to low-density scattering data, amongst other things, to find the appropriate Lennard-Jones parameters.Sciver12

Because of its simplicity and realism, the Lennard-Jones potential has a long history in statistical mechanics. In particular, since the early days of computer simulation it has been used to predict the physical properties of condensed matter systems, to establish the viability of various algorithms, to test the performance of different approximations and expansions, and to obtain benchmark results to which more approximate approaches could be compared. Barker71 ; Barker73 ; Rowley75 ; Barker76 ; Singer84 ; Allen87

One of the early concerns of these computer simulations was quantum effects, and how these modify the classical equation of state of the noble elements. Barker71 ; Barker73 ; Rowley75 ; Barker76 ; Singer84 ; Allen87 The primary quantum correction that was explored was the one given by Wigner and Kirkwood, which arises from the non-commutativity of the position and momentum operators, and which is second order in Planck’s constant.Wigner32 ; Kirkwood33 As is often the case, the first term in an expansion is very much simpler to deal with than the subsequent terms. Although Kirkwood gave a recursion formula for the higher order terms, Kirkwood33 the main focus of this earlier work was on the primary quantum correction.

The present paper also explores the quantum corrections to the classical equation of state using classical Monte Carlo simulations of Lennard-Jones fluids that model the noble elements. One advance on the earlier work is to include the second non-zero quantum correction (ie. the terms that are fourth order in Planck’s constant). Besides the numerical results, the analytic expansion by which this is accomplished is of some interest in its own right. Kirkwood’s expansionKirkwood33 is such that the higher order terms beyond the second order one are non-extensive, which poses difficulties since the grand potential, which the expansion is meant to yield, must be extensive. This means that Kirkwood’s expansion cannot be used beyond the second order without some sort of re-summation. Three related but different expansions to fourth order in Planck’s constant are given in the present paper. Unlike Kirkwood’s expansion, the present expansions are based on a quantity that is strictly extensive.

A second contribution of the present paper is to include fully the effects of wave function symmetrization.Attard16 KirkwoodKirkwood33 gave the primary correction in this series, but his correction differs from the first term used in the present computations by a factor of two due to double counting.Attard16b Although the effects of wave function symmetrization are small for the noble elements in the regime explored here, it is nevertheless of interest to establish that fact. It is also of interest to see how the effects of symmetrization and non-commutativity are combined.

II Formalism

II.1 Partition Function

In statistical mechanics, the partition function is the total number of states of the total system, and its logarithm is in essence the total entropy. Some care needs to be taken with the counting of states, so that forbidden states are not included, and so that the same state is not counted more than once.

The symmetrization of the wave function plays a rôle in this, as states that differ only by the permutation of the particle labels are not distinct. Further, states with repeated label values may be allowed for bosons but forbidden for fermions. The way to properly formulate the partition function to account for these effects is to define a characteristic function for the microstates.

II.1.1 Characteristic Function

A basis wave function formed from the orthonormal set of unsymmetrized wave functions $\{\phi_{\bf n}\}$ has symmetrized form

[TABLE]

Here $\hat{\mathrm{P}}$ is the permutation operator, $p$ is its parity, the upper sign is for bosons, and the lower sign is for fermions. The characteristic function of the state, $\chi^{\pm}_{\bf n}$, is inversely proportional to the number of non-zero distinct permutations of the wave function. Specifically, normalization, $\langle\phi_{\bf n}^{\pm}|\phi_{\bf n}^{\pm}\rangle=1$, gives

[TABLE]

With the characteristic function, the sum over distinct, allowed states of some symmetric function $f_{\hat{\mathrm{P}}{\bf n}}=f_{\bf n}$ can be written as a sum over all states,

[TABLE]

II.1.2 Example: Two Particles in Two States

As an example, consider a two particle system, each of which can exist in one of two one-particle states. One has

[TABLE]

Since two fermions cannot be in the same state, the 11 and 22 states are forbidden. Also

[TABLE]

The state 12 is the same as the state 21.

Hence the sum over distinct, allowed states of some symmetric function $f_{\hat{\mathrm{P}}{\bf n}}=f_{\bf n}$ can be written as the sum over all states by using the characteristic function,

[TABLE]

For fermions, the 11 and 22 states are excluded. For bosons the 11 and 22 states are each counted with weight 1. The $12=21$ state is counted once for both bosons and fermions.

II.1.3 Partition Function

Let $\phi_{\bf n}$ be an entropy (energy) eigenfunction, $\hat{\cal H}|\phi_{\bf n}\rangle={\cal H}_{\bf n}|\phi_{\bf n}\rangle$. Here $\hat{\cal H}$ is the Hamiltonian or energy operator, and the $\{\phi_{\bf n}\}$ form a complete orthonormal unsymmetrized set, and ${\bf n}$ labels entropy (energy) microstates.

As just mentioned, the partition function is the total number of allowed, distinct states of the total system, which is the total number of reservoir-weighted allowed, distinct states of the sub-system. The emphasis is on the words allowed and distinct, because it would be wrong to count forbidden states, or to count the same state more than once. Hence the grand partition function isAttard16 ; Attard16b

[TABLE]

Here $N$ is the number of particles, the fugacity is $z\equiv e^{\beta\mu}$, where $\mu$ is the chemical potential, and $\beta=1/k_{\mathrm{B}}T$ is sometimes called the inverse temperature, with $k_{\mathrm{B}}$ being Boltzmann’s constant and $T$ the temperature.

The standard expression for the partition function is as the trace of the Maxwell-Boltzmann operator or density matrix. Neumann27 ; Messiah61 ; Merzbacher70 It can be derived from the collapse of the sub-system wave function into entropy microstates due to entanglement with the reservoir.QSM1 ; QSM Unlike the standard expression, the present expression explicitly assigns the correct weight to forbidden states and to duplicate states.Attard16 The present expression is dependent upon using entropy microstates and their characteristic function $\chi^{\pm}(\phi_{\bf n})$. It does not appear possible to write the partition function as the trace of an operator without altering the definition of trace.

Consider another orthonormal, complete, unsymmetrized but otherwise at this stage arbitrary basis $\{\zeta_{\bf p}\}$. Since this is complete, $\sum_{\bf p}|\zeta_{\bf p}\rangle\,\langle\zeta_{\bf p}|=\hat{\mathrm{I}}$, the identity operator in this form can be inserted to yield

[TABLE]

The second equality follows from the replacements ${\bf n}\Rightarrow\hat{\mathrm{P}}{\bf n}$ and ${\bf p}\Rightarrow\hat{\mathrm{P}}{\bf p}$, and the fact that $\langle\zeta_{\hat{\mathrm{P}}{\bf p}}|\phi_{\hat{\mathrm{P}}{\bf n}}\rangle=\langle\zeta_{\bf p}|\phi_{\bf n}\rangle$. The final equality has the same form as the final form of the sum over entropy states.

Finally consider yet another orthonormal, complete, unsymmetrized, but otherwise at this stage arbitrary basis $\{\zeta_{\bf q}\}$, with $\sum_{\bf q}|\zeta_{\bf q}\rangle\,\langle\zeta_{\bf q}|=\hat{\mathrm{I}}$ One can write

[TABLE]

(One could replace here ${\bf q}\Rightarrow{\bf p}^{\prime}$, in which case the final form is just the non-diagonal sum over a single basis set, analogous to the diagonal sum over the entropy basis set.)

II.1.4 Momentum and Position States

Now particularize the analysis to momentum and position basis functions. Attard16b The momentum basis functions in the position representation ${\bf r}$ are plane waves,

[TABLE]

Here ${\bf p}$ is the configuration momentum. These form a complete orthonormal set. Periodic boundary conditions give the width of the momentum state per particle per dimension as $\Delta_{p}=2\pi\hbar/V^{1/3}$.

The position basis functions are Gaussians,

[TABLE]

Here ${\bf q}$ is the configuration position. Normalization fixes the spacing of the configuration position states as $\Delta_{q}=\sqrt{8\pi\xi^{2}}$, which gives the completeness expression, $\sum_{\bf q}|\zeta_{\bf q}({\bf r}^{\prime})\rangle\,\langle\zeta_{\bf q}({\bf r})|={e^{-({\bf r}^{\prime}-{\bf r})^{2}/8\xi^{2}}}/{(8\pi\xi^{2})^{3N/2}}\equiv\delta_{\xi}({\bf r}^{\prime}-{\bf r})$. In the limit $\xi\rightarrow 0$ the position basis functions form a complete orthonormal set.Attard16b

The transformation coefficient for the two basis sets is

[TABLE]

It is important for the treatment of permutation loops below that this is the product of individual particle factors. The product of this and its complex conjugate is obviously

[TABLE]

which product will shortly appear in the summand of the partition function.

II.1.5 Non-Commutativity Factor

The position and momentum operators do not commute, which fundamental property can be accounted for by recasting the partition function to include the so-called ‘non-commutativity factor’.

Modifying slightly an argument due to Wigner,Wigner32 one can formally commute the Maxwell-Boltzmann operator with the momentum basis function by writingAttard16b

[TABLE]

The modified energy operator induced here is

[TABLE]

The so-called non-commutativity function $w({\bf r},{\bf p})$ is extensive. Attard16b (This is one reason that the present formulation is preferable to those given by WignerWigner32 or by Kirkwood.Kirkwood33 The other reason is accounting for wave function symmetrization, which rigorously counts forbidden and multiple states correctly.) A recursion relation leading to an expansion for it will be given in §II.2 below. To leading order it vanishes, $w(\hbar=0)=w(\beta=0)=0$.

With this non-commutativity function, the expectation value of the Maxwell-Boltzmann operator that appears in the position-momentum space formulation of the partition function becomes

[TABLE]

The final equality follows because $\lim_{\xi\rightarrow 0}\zeta_{\bf q}({\bf r})=\delta({\bf q}-{\bf r})$.

II.1.6 Partition Function

The obvious merit of formulating the grand partition function in terms of the asymmetric expectation value of the Maxwell-Boltzmann operator is that the sum over entropy states has become a sum over points in classical phase space. The continuum limit of this is

[TABLE]

The volume elements are $\Delta_{p}=2\pi\hbar/V^{1/3}$ and $\Delta_{q}=\sqrt{8\pi\xi^{2}}$.Attard16b

This way of representing quantum mechanics in classical phase space is distinctly different to the way advocated by Wigner,Wigner32 Kirkwood,Kirkwood33 and followers. Groenewold46 ; Moyal49 ; Barnett03 ; Gerry05 ; Zachos05 ; Praxmeyer02 ; Dishlieva08 ; Barker10 In particular, previous authors focus upon a transform of the wave function rather than the states themselves, they do not distinguish the position representation ${\bf r}$ from the position configuration ${\bf q}$ as here, and they do not account for wave function symmetrization. Finally, and again in contrast to the present, their focus is neither the partition function nor quantum statistical mechanics.

In view of the above the partition function becomes

[TABLE]

Notice that the factors of $\xi$ cancel, which allows the limit $\xi\rightarrow 0$ to be taken. Here

[TABLE]

which is just the complex conjugate of the function $\tilde{\chi}^{\pm}({\bf\Gamma})$ defined previously.Attard16b This might be called the symmetrization factor.

II.1.7 Monomer Grand Potential

Shortly it will be shown that the symmetrization factor, $\eta^{\pm}$, breaks up into loops of particles, $l=1,2,\ldots$, which may be called monomers, dimers etc. The leading monomer term is unity, $\eta^{(1)}=1$. Hence the monomer grand partition function is

[TABLE]

This is the same for bosons as for fermions. The ratio of the full partition function to the monomer partition function is just the monomer average of the symmetrization factor,

[TABLE]

Since the leading order of the non-commutativity factor vanishes, $w(\hbar=0)=w(\beta=0)=0$ (see §II.2 below), the leading order of the monomer grand potential is just the classical grand potential,

[TABLE]

The logarithm of this gives the classical equilibrium grand potential, $\Omega_{1,0}=-k_{\mathrm{B}}T\ln\Xi_{1,0}$. Following the usual procedures of classical statistical mechanics, the difference between the quantum and the classical monomer grand potentials (ie. the quantum correction due to monomers) is just a classical equilibrium average,

[TABLE]

where $\Omega_{1}$ is the full monomer grand potential. Here and below the subscript 1,0 denotes the classical equilibrium average or classical equilibrium thermodynamic potential.

II.1.8 Permutation Loop Expansion of the Grand Potential

The permutation operator breaks up into loops

[TABLE]

The prime on the sums restrict them to unique loops, with each index being different. The first term is just the identity. The second term is a dimer loop, the third term is a trimer loop, and the fourth term is the product of two different dimers.

The symmetrization factor, $\eta^{\pm}({\bf\Gamma})$ $=\sum_{\hat{\mathrm{P}}}(\pm 1)^{p}$ $\langle\zeta_{\hat{\mathrm{P}}{\bf p}}|\zeta_{\bf q}\rangle/\langle\zeta_{{\bf p}}|\zeta_{\bf q}\rangle$, is the sum of the expectation values of these loops. Since the expectation value is the product of the individual particle factors, Eq. (2.21), the factors due to the unpermuted particles in the numerator cancel with those in the denominator.

The monomer symmetrization factor comes from just the unpermuted expectation value,

[TABLE]

The dimer overlap factor in the microstate ${\bf\Gamma}$ for particles $j$ and $k$ is

[TABLE]

The important point is the formally exact factorization of the expectation value, leaving only the permuted particles to contribute.

The symmetrization factors are localized in the sense that they are only non-zero when all the particles of the loop are close together. (More precisely, localization means that the separations between consecutive neighbors around the loop are all small.) This will be shown explicitly below, but here it can be noted that the explicit exponents give highly oscillatory and therefore canceling behavior unless the differences in configuration positions are all close to zero.

Similarly the trimer symmetrization factor for particles $j$, $k$, and $l$ is

[TABLE]

Continuing in this fashion, the symmetrization factor can be written as a series of loop products,

[TABLE]

Here the superscript is the order of the loop, and the subscripts are the atoms involved in the loop. This gives the ratio of the full to the monomer partition function as

[TABLE]

The third and following equalities write the average of the product as the product of the averages. This is valid in the thermodynamic limit, since the product of the average of two loops scales as $V^{2}$, whereas the correlated interaction of two loops scales as $V$. The combinatorial factor accounts for the number of unique loops in each term; $\eta^{\pm(l)}$ without subscripts refers to any one set of $l$ particles, since all sets give the same average.

The grand potential is essentially the logarithm of the grand partition function, $\Omega\equiv-k_{\mathrm{B}}T\ln\Xi$. Hence the difference between the full grand potential and the monomer grand potential is just the series of loop potentials,

[TABLE]

The monomer grand potential is of course $\Omega_{1}\equiv-k_{\mathrm{B}}T\ln\Xi_{1}$, with the monomer grand partition function being given by Eq. (2.29).

II.2 Expansion of the Non-Commutativity Function

Using a similar approach to Kirkwood,Kirkwood33 the temperature derivative of the defining equation (2.23), $e^{-\beta{\hat{\tilde{\cal H}}}}1=e^{-\beta{\cal H}}e^{w}$, may be rearranged to giveAttard16b

[TABLE]

Write

[TABLE]

with $w_{n}(\beta=0)=0$. (This begins at $n=1$ because the classical part must yield the classical Maxwell-Boltzmann factor, $w(\hbar=0)=0$.) This expansion leads to the recursion relation for $n>2$,

[TABLE]

It is straightforward to derive the first several coefficient functions explicitly. One has for $n=1$,

[TABLE]

for $n=2$,

[TABLE]

for $n=3$,

[TABLE]

and for $n=4$,

[TABLE]

Notice that the odd coefficient functions are pure imaginary and odd in momentum. Because ${\cal H}({\bf\Gamma})$ is an even function of momentum, the quantum weight $e^{w}$ is real and oscillatory on average. (The $\eta^{\pm(l)}$ contain terms that are either real and even, or else imaginary and odd in momentum.)

II.3 Monomer Expansion A

The quantum correction to the classical grand potential due to the monomers is just a classical average of the quantum weight due to non-commutativity, Eq. (2.32). To fourth order in $\hbar$ this is

[TABLE]

All higher order contributions have been set to zero here. It is necessary to write out the individual terms explicitly to cancel the odd imaginary terms. It is $w_{j}({\bf p},{\bf q})$ whose average is taken. This is extensive; the factor of 3 between the terms involving $w_{1}^{4}$ is a necessary part of the cancelation that accounts for the three ways of pairing the four momenta, as can be shown analytically and numerically (see, for example, §II.4.4 below). This expansion terminated at ${\cal O}(\hbar^{2})$ may be called A2, and at ${\cal O}(\hbar^{4})$ it may be called A4.

II.4 Monomer Expansion B

One can define the cumulative weight as

[TABLE]

and the $n$th approximation to the quantum correction to the monomer grand potential as

[TABLE]

One has $\lim_{n\rightarrow\infty}\Delta\Omega_{1}^{(n)}=\Omega_{1}-\Omega_{1,0}$.

II.4.1 Momentum Averages

The classical average $\left<\ldots\right>_{1,0}$ includes an average over the momenta $\left<\ldots\right>_{1,0,p}$ as well as one over the position configurations.

The classical monomer probability distribution for momentum is a Gaussian,

[TABLE]

Hence

[TABLE]

The thermal wave length is

[TABLE]

II.4.2 First order term

The momentum average of the first order term is

[TABLE]

Terms with an odd number of momenta vanish upon averaging. The combinatorial factor in the third equality arises from the number of ways of arranging the momenta in pairs: there are $(2l)!$ ways of arranging the momenta, of which interchanging the $l!$ pairs leaves the result unchanged, as does the $2^{l}$ arrangements that arise from swapping the first and second members of each of the pairs. In the third equality the thermodynamic limit has been assumed, which means that $N-1$ has been taken to equal $N$, etc. The final result is nice because it shows that the gradient of the potential is Gaussian distributed, which is to say that configurations with very large gradients are suppressed.

This gives the first order result

[TABLE]

II.4.3 Second Order Order Term

In the second order approximation to the weight, $w^{(2)}=\hbar w_{1}+\hbar^{2}w_{2}$, the terms that depend upon the momenta may be identified,

[TABLE]

This leaves the position configuration terms to be added as

[TABLE]

It is simplest to deal with the Gaussian exponent directly by completing the squares. Writing ${\bf U}^{\prime}\equiv\nabla U$ and ${\bf U}^{\prime\prime}\equiv\nabla\nabla U$, the momentum part of the exponent is

[TABLE]

where ${\bf A}\equiv\mathrm{I}-{\hbar^{2}\beta^{2}}{\bf U}^{\prime\prime}/{3m}$. One can write

[TABLE]

Since the sum of the squares of the eigenvalues is the trace of the square of the matrix, one has

[TABLE]

With these the momentum average is

[TABLE]

This exponent can be added to $w^{(2)}_{q}$ and the configuration position average taken to obtain $\Delta\Omega_{1}^{(2)}$ to ${\cal O}(\hbar^{6})$. This agrees with expansion A to ${\cal O}(\hbar^{2})$.

II.4.4 Terms ${\cal O}(\hbar^{4})$

In order to get the terms ${\cal O}(\hbar^{4})$ one needs $w_{3}$. The contribution from the dyadic triplet ${\bf p}{\bf p}{\bf p}$ can be obtained by setting one of the momenta to its most likely value,

[TABLE]

and the remaining two factors equal to their average

[TABLE]

There are three ways of doing this, and so $w_{3}$ becomes

[TABLE]

One can do a similar trick for $w_{4}$,

[TABLE]

This is multiplied by $\hbar^{4}$. The neglected contributions are of higher order in $\hbar$ (ie.  ${\cal O}(\hbar^{6})$).

In summary, the exponent that one has to average over position configurations to ${\cal O}(\hbar^{4})$ is

[TABLE]

This expansion terminated at ${\cal O}(\hbar^{2})$ may be called B2, and at ${\cal O}(\hbar^{4})$ it may be called B4.

II.5 Monomer Expansion C

Expansion A has the merit of requiring less analysis, and of requiring the classical average of extensive terms and their products. It has the disadvantage of requiring the average to be taken numerically over the momentum configurations. Also, as a fluctuation expression, it requires relatively high accuracy in the individual averages of products to get the necessary cancelation between these super-extensive terms to end up with the final extensive result for $\Omega_{1}-\Omega_{1,0}$. The larger the system size, the greater the number of configurations that need to be generated to get acceptable statistical accuracy.

Expression B has the advantage of not requiring numerical momentum averaging, which leads to higher accuracy because certain essential cancelations are enforced analytically. (In most of the results reported below, 32 momentum configurations were used for each position configuration for the average of expansion A; no momentum configurations are required for expansion B.) Also, by exponentiating the expansion, as in B, one expects faster convergence than for an expansion of the final expression, as in A. A disadvantage of expression B has is that it requires more explicit algebra. Also, more problematic, since $w$ is an extensive variable, taking the average of $e^{w}$ can lead to computational overflow problems for large systems.

One way to avoid numerical overflow in evaluating the exponent in expansion B, but to preserve the advantage of no numerical momentum average, is to use the momentum averaged $w^{(4)}$ in the expansion A. That is, write Eq. (II.4.4) as explicit powers of Planck’s constant,

[TABLE]

One can now insert this into expansion A, Eq. (II.3), with $w_{1}=w_{3}=0$,

[TABLE]

This may be called expansion C4. The averages here are classical averages over configuration positions. The neglected terms are ${\cal O}(\hbar^{6})$. Hence this can be expected to be rather similar to expansion A4, the only difference being that the momentum averages have been performed analytically to leading order before expanding the exponential.

II.6 Dimer

In general the loop potential $l\geq 2$ is given by Eq. (2.39),

[TABLE]

The zeroth order approximation for non-commutativity sets $w=0$ and is the one for which numerical results are discussed below. However the case $w=w^{(2)}$ is now analyzed.

For the dimer $l=2$, Eq. (2.35) is

[TABLE]

Combining this exponent with the part of $\hbar w_{1}({\bf\Gamma}^{N})$ that depends on ${\bf p}_{j}$ and ${\bf p}_{k}$ gives exponent

[TABLE]

One sees that the analysis of the momentum averages in §II.4.1 goes through with the only change being effective potential gradients for particles $j$ and $k$,

[TABLE]

respectively. From the zeroth order dimer result, $q_{jk}\sim\Lambda\sim\hbar$.

With these effective gradients, the momentum part of the exponent given by Eq. (2.57) becomes

[TABLE]

with the matrix ${\bf A}$ unchanged. Hence the average over the momenta, Eq. (2.60), becomes

[TABLE]

This expression is half-way between B2 and B4. Arguably, it would be most consistent to neglect here the terms ${\cal O}(\hbar^{4})$.

To this should be added the position configuration terms,

[TABLE]

The zeroth order dimer result is in the very first term, $({-\hbar^{2}\beta^{3}}/{8m})\nabla\tilde{U}\cdot\nabla\tilde{U}$, which upon setting $U=0$ becomes $-mq_{jk}^{2}/\hbar^{2}\beta=-2\pi q_{jk}^{2}/\Lambda^{2}$.

III Methodology and Results

III.1 Simulation Algorithm

Computer simulations were performed for a classical canonical equilibrium system in configuration position space.Allen87 The Metropolis algorithm was used, with the step length adjusted to give approximately a 50% acceptance rate.

A small cell, spatially based, linked list, neighbor table was used. Attard04 The cells had edge length on the order of 3/4 of the molecular diameter, which means that most cells were empty, and almost no cells were occupied by more than one particle. The total volume occupied by neighbor cells was on the order of 2–5 times the volume of the potential cut-off sphere. Tests with 250–2000 atoms confirmed that the time cost of the algorithm scaled linearly with the system size.

A cycle consists of an attempted move for every particle in turn. After every 20 cycles, the position configuration was stored for use in later averaging. For most of the results reported below, 20,000 configurations were saved. For each temperature and density, a number of equilibration cycles were run, starting from a previously equilibrated configuration at a nearby state point. After equilibration was judged complete, generally on the basis of the stability of the virial pressure, the simulation proper commenced.

For averaging, the configurations were broken into 50 blocks, and the statistical error was estimated from the fluctuations of the average of each block. The error reported below corresponds to one standard deviation, which corresponds to 68% confidence in the result.

The expansion A requires averaging over the momentum configurations in addition to averaging over the position configurations that have been saved. (The expansions B and C do not require this, but these were calculated at the same time as A.) For this purpose, for each position configuration, a number, generally 32, of independent momentum configurations were generated. Each momentum configuration corresponds to a new momentum for every atom. These were drawn randomly from a Gaussian distribution that corresponds to the classical Maxwell-Boltzmann distribution of the kinetic energy.

The position configurations do not depend upon the mass of the atoms, and the same position configurations were used for all three noble elements at that particular (dimensionless) state point. The momenta do depend upon the atomic mass, and the averaging program was run three times at each state point, once for each noble element.

The Lennard-Jones pair potential was used (see next) and cut off at $R_{\mathrm{cut}}=3.5\sigma$. Periodic boundary conditions were enforced, and the edge length of the cubic simulation cell was $L>2R_{\mathrm{cut}}$. For the case of the virial pressure and the internal energy (not reported), a tail correction was included in the reported average. For the case of the quantum correction to quadratic order in $\hbar$, a tail correction was used for earlier results.Attard16b However, the quantum corrections reported here do not include a tail correction.

III.2 Lennard-Jones Potential

The Lennard-Jones potential used here is

[TABLE]

The Lennard-Jones parameters for the noble elements used in the simulations are shown in Table 1.

For the Lennard-Jones fluid, the critical temperature, density, and pressure in dimensionless form areTDSM

[TABLE]

respectively.

III.3 Pair Potential

For the case of a general pair potential, including the present Lennard-Jones potential, the potential energy is

[TABLE]

with ${\bf q}_{jk}={\bf q}_{j}-{\bf q}_{k}$ and ${q}_{jk}=|{\bf q}_{jk}|$.

The gradient is

[TABLE]

where $\alpha=x,y,\mbox{ or }z$. Hence

[TABLE]

The dyadic second derivative is

[TABLE]

from which the Laplacian is

[TABLE]

One also has

[TABLE]

and

[TABLE]

It was found convenient to store $\nabla U$, $\nabla\nabla^{2}U$, and the diagonal elements of $\nabla\nabla U$ as vectors for use in the evaluation of $w$. This keeps the scaling of the algorithm linear with system size.

III.4 Results

The expansions for the quantum correction, A2, A4, B2, B4, C2, and C4, were implemented and their average was obtained for various temperatures, densities, and noble elements. It was confirmed numerically that the results were extensive with system size. Most results are for $N=1000$ atoms, but in some cases this was reduced to $N=500$ to reduce the statistical error in the expansion A, or to avoid overflow errors in expansion B. The expressions B and C in general had very low statistical error. The expressions A2 and A4 had a higher statistical error, particularly A4, and particularly for the larger systems. These two expressions require numerical cancelations in the fluctuation terms to get the correct extensivity. The statistical error in the A expansions was reduced by increasing the number of momentum averages per position configuration to 32. The optimum number of momentum averages was not systematically determined. (The B and C expressions do not require explicit momentum averaging.) There was good agreement between A2 and W12a and W12b,Attard16b which are the same as the Wigner and Kirkwood expressions.

Figure 1 gives results for argon for the super-critical isotherm $T^{*}=1.5$. At low densities, where the fourth order contribution is small compared to the second order, there is good agreement between A2 and B2 and between A4 and B4. This gives one high confidence in the algebra and computer programming, since these two expressions and their implementation are largely independent of each other. The programming for C is dependent on that for B.

As the density is increased $\rho\sigma^{3}\gtrsim 0.5$, the statistical error in A4 becomes almost unacceptable. At $\rho\sigma^{3}=$0.8, 32 momentum configurations per position configuration were used. Reducing the number of momentum configurations by a factor of 4 increased the error by a factor of 1.8, which suggest that this is the limiting factor. (At $\rho\sigma^{3}=0.7$, and lower, 8 momentum configurations per position configuration were used.) The data at low to moderate densities shows that the first quantum correction, A2 and B2 (not shown separately), to the pressure is negative. The second quantum correction, A4$-$A2 and B4$-$B2, is positive, but the total correction, A4 and B4, is still negative. The magnitudes of the corrections increase with increasing density.

At higher densities, $\rho\sigma^{3}\gtrsim 0.7$, there is possibly a difference between A4 and B4 beyond the statistical error. Although both are fourth order in $\hbar$, the A expression requires three successive expansions (of $w$, of $e^{w}$, and of $\ln\langle e^{w}\rangle$) before truncation, whereas the expression B performs the truncation on $w$ before evaluating $\ln\langle e^{w}\rangle$ explicitly. At higher densities the terms beyond fourth order in $\hbar$ become significant, and one might expect the results of the three expressions to differ. However, for argon at this temperature all three expansions lie within the statistical error even at the highest density.

In those cases where the quantum correction A4 agrees with B4, the higher order terms must be negligible and either gives a reliable estimate of the total quantum correction. When the two differ, the difference is probably of similar order to the neglected higher order term, and it is likely a reasonable guide to their accuracy as an estimate of the total quantum correction.

Fig. 2 gives results for neon for the isotherm $T^{*}=1.5$. At low densities, $\rho\sigma^{3}\lesssim 0.5$, there is good agreement between A4, B4, and C4. At higher densities the statistical error in A4 is noticeably larger. The data for B4 and C4 form smooth curves that indicate a monotonic increase in the quantum correction as the density is increased. (No statistical error was obtained for B4 for $\rho\sigma^{3}\geq 0.7$, presumably because of overflow error. At $\rho\sigma^{3}=0.6$ the relative statistical error in B4 was 0.3%.) The statistical error for C4 was quite small, being about 0.1% at the highest density shown. Whereas the quantum correction was negative for argon, it is positive for neon and larger in magnitude.

For neon at $T^{*}=1.5$, the quantum correction is comparable in magnitude to the classical pressure itself (Fig. 2, solid curve), whereas for argon, the quantum correction is relatively negligible. As a percentage of the classical pressure, the quantum correction B4 at $\rho\sigma^{3}=0.5$ is 1% for argon and 42% for neon. At $\rho\sigma^{3}=0.8$ it is 0.4% for argon and 25% for neon. At $\rho\sigma^{3}=0.2$ it is 0.3% for argon and 21% for neon.

Figure 3 gives the quantum correction to the classical grand potential per unit volume for helium, again at $T^{*}=1.5$. In this case only the expansions A4 and C4 are shown as computer overflow precluded results for the expansion B4. In these simulations, 500 particles, $4\times 10^{5}$ position configurations, and 32 momentum configurations per position configuration were used for the A averages. For the C averages, 1000 particles, $2\times 10^{5}$ position configurations, and 1 momentum configuration per position configuration were used.

At all densities on this isotherm for helium, there was a significance difference between A2 and B2, and between A2 and A4. (In fact, A2 is negative, whereas A4 is positive. Similarly, C2 is negative and C4 is positive and about a factor of 10 larger in magnitude.) This suggests that one should be cautious about relying on the results in Fig. 3 as an estimate of the total quantum correction. The fact that the quantum correction A4 is a couple of hundred times larger than the classical virial pressure itself also suggests that these results for helium should be used cautiously.

Figure 4 shows results for argon on the liquid branch of the sub-critical isotherm $T^{*}=1$. The relatively good agreement between A4 and B4 in this case suggests that the excluded higher order contributions are probably negligible. At the highest densities shown, there is a discrepancy between A4, B4, and C4 which suggests that the excluded higher order contributions are required to make a quantitatively accurate estimate of the full quantum correction.

Over the limited density range the quantum correction is negative with relatively little variation. This contrasts with argon at $T^{*}=1.5$ where the correction is negative and monotonic increasing in magnitude with increasing density. Except at the lowest density shown, where the classical pressure at $T^{*}=1$ is close to zero, the quantum correction for argon amounts to less than 1% of the classical pressure (see Fig. 5).

Figure 5 shows results for neon at $T^{*}=1$. In this case the quantum correction is positive and increases with increasing density. Only the results for A4 and C4 are shown, as overflow error precluded results for B4. In these simulations for A, 500 particles, $4\times 10^{5}$ position configurations, and 32 momentum configurations per position configuration were used for the averages. For C, 1000 particles, $2\times 10^{5}$ position configurations, and 1 momentum configuration per position configuration were used for the averages. The results for B2 are relatively close to those for A2, but the results for A4 are of opposite sign to those for A2, and are about a factor of 10 larger in magnitude. Similarly, the results for C4 are of opposite sign to those for C2, and are about a factor of 10 larger in magnitude. Again this suggests that one should be cautious in assuming that the excluded higher order terms are negligible. One can see that for densities $\rho\sigma^{3}\lesssim 0.8$ the quantum correction A4 is larger than the classical pressure. For larger densities the two are comparable.

Figure 6 shows A4 and C4 results for helium at $T^{*}=1$. (Overflow error precluded B4 results.) These are positive and monotonically increase with increasing density. The fact that the quantum corrections are 500 or more times greater than the corresponding classical pressure makes it questionable whether truncating the expansion at A4 yields a trustworthy estimate of the total quantum correction for helium on this sub-critical isotherm. Further doubt arises as the corrections appear to alternate in sign, A2 and C2 being negative and monotonic increasing in magnitude with increasing density in this case. (The magnitude of C4 is about 200 times that of C2 at the highest density shown.)

Figure 7 is for $k_{\mathrm{B}}T/\varepsilon=0.8$, with results for the classical virial pressure and the C4 quantum correction for argon, neon, and helium shown. In this and the following figures the quantum correction for argon has been multiplied by 100, and that for helium has been divided by 1000. Also, in these figures, 1000 atoms were used, 5,000 position configurations were saved for averaging, and one momentum configuration per position configuration was used for the A2 and A4 average (not shown). The quantum correction for argon is negative at the gas-end of the liquid branch isotherm, but then turns increasingly positive. The corrections for neon and helium are positive and monotonic increasing with increasing density. Of course, where the classical virial pressure is close to zero, the quantum corrections can be rather larger in relative terms, even though they are small in absolute terms.

Results at $k_{\mathrm{B}}T/\varepsilon=0.6$ are shown in Fig. 8. Regions where the pressure has negative slope correspond to unstable states in the thermodynamic limit. Thus one can see from the classical virial pressure, that there is a gas phase for $\rho\sigma^{3}\lesssim 0.1$, a liquid phase $0.8\lesssim\rho\sigma^{3}\lesssim 0.95$, and a solid phase $\rho\sigma^{3}\gtrsim 1$. Inclusion of the quantum correction may shift the phase boundaries. The C4 quantum correction is qualitatively similar for all three noble elements, with it being positive and mainly increasing with increasing density. The quantum correction for neon is now somewhat larger than the classical virial presure

Figure 9 shows results for $k_{\mathrm{B}}T/\varepsilon=0.5$. This is the lowest temperature for which results were obtained. Again the three phases are evident in the classical virial pressure. The C4 quantum correction is about a factor of 100 smaller in magnitude than the classical virial pressure, and that for helium about a factor of 1000 larger. However, the C4 quantum correction for neon is about the same as, or somewhat greater than the classical virial pressure. At this temperature the C4 quantum correction is positive and mainly increases as the density increases for all three noble elements.

The quantum correction due to wave function symmetrization was also simulated for dimer loops, both without and with the weight for non-commutativity. In all cases they were found to be negligible compared to the monomer non-commutativity quantum corrections. For this reason they are not included in the above figures.

IV Summary and Conclusion

In this paper a formal expression has been given for the grand potential in quantum statistical mechanics as a perturbation expansion about the classical grand partition function whose terms involve classical averages. There were two contributions to the expansion. The first arises from the non-commutativity of position and momentum operators. The second arises from the symmetrization of the wave function for fermions and for bosons.

The expansion for non-commutativity given here is similar to that given by Wigner and by Kirkwood.Wigner32 ; Kirkwood33 The main improvement here is that the phase space function that is used for the expansion is an extensive variable and the terms may be directly related to the grand potential, which is of course also extensive. To leading order in the quantum correction for non-commutativity (ie. ${\cal O}(\hbar^{2})$), the present expansions agree with those of Wigner and Kirkwood.Wigner32 ; Kirkwood33 The fourth order terms given explicitly here are new.

The loop expansion for wave function symmetrization given here was given previously by the present author.Attard16 ; Attard16b The present derivation is simpler but not materially different to the earlier ones. An important point to emerge from the present treatment is that in so far as the total entropy is the logarithm of the weighted sum of distinct allowed states, with the emphasis on distinct and allowed, the partition function of quantum statistical mechanics cannot be written as the trace of a density operator, at least not without an explicit modification of the definition of trace.

The quantum corrections proved straightforward to implement in a computer simulation. This is not so surprising as they were designed as classical statistical averages. Of the three expansions given, expansion C is arguably the best. In it the momentum averages have been carried out explicitly, which improves the statistical precision by more than an order of magnitude for more than a factor of 30 less computer time compared to expansion A. Compared to expansion B, it does not exponentiate an extensive variable, and therefore it is not susceptible to computational overflow problems. The downside is that C involves expanding the exponential and discarding some higher order terms that are retained in B.

Results were obtained for the three smallest noble elements on one supercritical and four subcritical isotherms. In general, at higher densities, the quantum correction for argon was about a factor of 100 smaller than, that for neon was of the same order as, and that for helium was about a factor of 1000 larger than the classical virial pressure. The fourth order quantum correction calculated here as an estimate of the total quantum correction appears reliable for argon, marginal for neon, and unreliable for helium.

The above discussion refers to the quantum correction for non-commutativity. The correction for wave function symmetrization was negligible for all cases studied here.

The present expansions retaining the fourth order terms for non-commutativity are rather tedious to extend to higher order terms. However the numerical results, particularly for helium, are rather slowly converging, and it appears that there is a need to go beyond the fourth order. To that end it would be nice if a more refined approach not based on a crude expansion could be found.

Finally, the quantum effects due to wave function symmetrization were found to be negligible for the noble elements in the present regime. It would be interesting to quantitatively analyze a case where symmetrization effects were dominant.