Table Of ContentQuantum Statistical Mechanics Results for Argon, Neon, and Helium Using Classical
Monte Carlo
Phil Attard
(Dated: 1 Feb., 2017, [email protected])
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
virialpressureonallisothermsatliquid-likedensities,thequantumcorrectionisfoundtobeO(10−2)
for argon, O(100) for neon, and O(103) for helium. Thefirst order correction dueto wavefunction
symmetrization is also obtained, butthis is relatively negligible.
7
1
I. INTRODUCTION model the noble elements. One advance on the earlier
0
work is to include the second non-zero quantum correc-
2
The Lennard-Jones potential is arguably the simplest tion (ie. the terms that are fourth order in Planck’s con-
b
stant). Besides the numerical results, the analytic ex-
e inter-particle potential that contains the two features of
pansionbywhichthis is accomplishedis ofsome interest
F realistic atoms and molecules that are necessary for the
properties of condensed matter. These are a short-range in its own right. Kirkwood’s expansion10 is such that
1
the higher order terms beyond the second order one are
repulsion that reflects the particle size, and a long-range
non-extensive,whichposesdifficultiessincethegrandpo-
] attractionthatis responsibleforthe cohesionofthe par-
h tential, which the expansion is meant to yield, must be
ticles and for the existence of the liquid phase, as well
p extensive. This means that Kirkwood’s expansion can-
as for the nature of the solid phase. That the attraction
- notbeusedbeyondthesecondorderwithoutsomesortof
t decays with the sixth power of separation is due to the
n re-summation. Three related but different expansions to
correlated electronic fluctuations that induce temporary
a
fourthorderinPlanck’sconstantaregiveninthepresent
u dipolesaccordingtotheirpolarizability,andagainthisis
paper. Unlike Kirkwood’s expansion, the present expan-
q a feature of all real atoms and molecules.
[ Although extensive tabulations of Lennard-Jones pa- sions are based on a quantity that is strictly extensive.
1 rameters for various molecules exist,1 the potential is A second contribution of the present paper is to in-
cludefullytheeffectsofwavefunctionsymmetrization.11
v mostsuitedforsphericalatomsandmoleculessinceitisa
Kirkwood10 gave the primary correction in this series,
6 function ofthe separationalone. The noble gasatomsin
9 particularhavebeenfittedtolow-densityscatteringdata, but his correction differs from the first term used in the
0 amongst other things, to find the appropriate Lennard- present computations by a factor of two due to double
0 Jones parameters.2 counting.12 Although the effects of wave function sym-
0 metrizationaresmallforthenobleelementsintheregime
Because of its simplicity and realism, the Lennard-
.
2 explored here, it is nevertheless of interest to establish
Jones potential has a long history in statistical mechan-
0 that fact. It is also of interest to see how the effects of
ics. In particular, since the early days of computer sim-
7 symmetrization and non-commutativity are combined.
ulation it has been used to predict the physical proper-
1
: ties of condensed matter systems, to establish the vi-
v
ability of various algorithms, to test the performance
i
X of different approximations and expansions, and to ob-
II. FORMALISM
tain benchmark results to which more approximate ap-
r
a proaches could be compared.3–8
A. Partition Function
One of the early concerns of these computer simula-
tions was quantum effects, and how these modify the
classical equation of state of the noble elements.3–8 The In statistical mechanics, the partition function is the
primary quantum correction that was explored was the total number of states of the total system, and its loga-
onegivenbyWignerandKirkwood,whicharisesfromthe rithmis in essencethe totalentropy. Some careneeds to
non-commutativityofthe positionandmomentumoper- be taken with the counting of states, so that forbidden
ators,andwhichissecondorderinPlanck’sconstant.9,10 statesarenotincluded, andsothatthesamestateisnot
Asisoftenthecase,thefirstterminanexpansionisvery counted more than once.
muchsimplertodealwiththanthesubsequentterms. Al- The symmetrization of the wave function plays a rˆole
thoughKirkwoodgavearecursionformulaforthehigher in this, as states that differ only by the permutation of
orderterms,10 the mainfocus ofthis earlierworkwas on the particle labels are not distinct. Further, states with
the primary quantum correction. repeated label values may be allowedfor bosons but for-
The present paper also explores the quantum correc- bidden for fermions. The way to properly formulate the
tions to the classical equation of state using classical partitionfunctionto accountforthese effects isto define
Monte Carlo simulations of Lennard-Jones fluids that a characteristic function for the microstates.
2
1. Characteristic Function For fermions, the 11 and 22 states are excluded. For
bosonsthe11and22statesareeachcountedwithweight
Abasiswavefunctionformedfromtheorthonormalset 1. The12=21stateiscountedonceforbothbosonsand
of unsymmetrized wave functions {φ } has symmetrized fermions.
n
form
1
φ± = (±1)pφ . (2.1)
n N!χ± Pˆn
n XPˆ
HerePˆ isthepermutpationoperator,pisitsparity,theup- 3. Partition Function
per sign is for bosons, and the lower sign is for fermions.
The characteristic function of the state, χ±n, is inversely Letφn be anentropy(energy)eigenfunction, Hˆ|φni=
proportional to the number of non-zero distinct permu- H |φ i. Here Hˆ is the Hamiltonian or energy opera-
n n
tations of the wave function. Specifically, normalization,
tor, and the {φ } form a complete orthonormal unsym-
hφ±|φ±i=1, gives n
n n metrized set, and n labels entropy (energy) microstates.
χ± = (±1)phφ |φ i. (2.2) As just mentioned, the partition function is the total
n n Pˆn
XPˆ number of allowed, distinct states of the total system,
which is the total number of reservoir-weightedallowed,
Withthecharacteristicfunction,thesumoverdistinct,
distinctstatesofthesub-system. Theemphasisisonthe
allowed states of some symmetric function f =f can
Pˆn n wordsallowedanddistinct,becauseitwouldbewrongto
be written as a sum over all states,
count forbidden states, or to count the same state more
1 than once. Hence the grand partition function is11,12
′f = χ±f . (2.3)
n N! n n
n n
X X
∞
Ξ± = zN ′e−βHn
2. Example: Two Particles in Two States
N=0 n
X X
∞ zN
As an example, consider a two particle system, each = N! χ±ne−βHn
of which can exist in one of two one-particle states. One N=0 n
X X
has ∞ zN
= (±1)phφ |φ ie−βHn
χ±11 = hφ11|φ11i±hφ11|φ11i NX=0 N! Xn XPˆ Pˆn n
= (cid:26) 20 = ∞ zNN! (±1)phφPˆn|e−βHˆ|φni. (2.7)
= χ±. (2.4) NX=0 Xn XPˆ
22
Since two fermions cannot be in the same state, the 11
HereN isthenumberofparticles,thefugacityisz ≡eβµ,
and 22 states are forbidden. Also
where µ is the chemical potential, and β = 1/k T is
B
χ± = hφ |φ i±hφ |φ i sometimes calledthe inversetemperature,with k being
12 12 12 12 21 B
= 1 Boltzmann’s constant and T the temperature.
= χ±. (2.5) Thestandardexpressionforthepartitionfunctionisas
21
the trace of the Maxwell-Boltzmann operator or density
The state 12 is the same as the state 21. matrix.13–15 It can be derived from the collapse of the
Hence the sum over distinct, allowed states of some
sub-system wave function into entropy microstates due
symmetric function fPˆn =fn can be written as the sum to entanglement with the reservoir.16,17 Unlike the stan-
over all states by using the characteristic function,
dardexpression,the presentexpressionexplicitly assigns
the correct weight to forbidden states and to duplicate
f +f +f
′fn = f11 12 22 states.11 The present expression is dependent upon us-
12
n (cid:26) ing entropy microstates and their characteristic function
X
f + 1[f +f ]+f χ±(φ ). It does not appear possible to write the parti-
= 11 2 12 21 22 n
1[f +f ] tionfunctionasthetraceofanoperatorwithoutaltering
(cid:26) 2 12 21
the definition of trace.
1
= χ±f +χ±f +χ±f +χ±f
2 11 11 12 12 21 21 22 22 Consider another orthonormal, complete, unsym-
1(cid:8) (cid:9) metrizedbutotherwiseatthisstagearbitrarybasis{ζ }.
= χ±f . (2.6) p
N! n n Since this is complete, |ζ ihζ | =ˆI, the identity op-
n p p p
X
P
3
erator in this form can be inserted to yield as ∆ = 8πξ2, which gives the completeness expres-
q
sion, |ζ (r′)ihζ (r)| = e−(r′−r)2/8ξ2/(8πξ2)3N/2 ≡
Ξ± = ∞ zNN! (±1)phφPˆn|e−βHˆ|ζpihζp|φni δtiξo(nrs′ −fPorrqm).paqIncotmhpelelqitmeiotrtξho→no0rmtahlespeto.s1i2tion basis func-
NX=0 Xn Xp XPˆ
The transformationcoefficient for the two basis sets is
∞ zN
= (±1)phφ |e−βHˆ|ζ ihζ |φ i
N! n Pˆp p n (8πξ2)3N/4
NX=0 Xn Xp XPˆ hζp|ζqi ≡ VN/2 e−ξ2p2/h¯2eq·p/ih¯
∞ zN
= N! (±1)phζp|e−βHˆ|ζPˆpi. (2.8) = N (8πξ2)3/4e−ξ2p2j/h¯2eqj·pj/ih¯. (2.12)
NX=0 Xp XPˆ V/2
j=1
Y
The second equality follows from the replacements n ⇒
Pˆn and p⇒Pˆp, and the fact that hζ |φ i=hζ |φ i. It is important for the treatment of permutation loops
Pˆp Pˆn p n
The final equality has the same form as the final formof below that this is the product of individual particle fac-
tors. The product of this and its complex conjugate is
the sum over entropy states.
obviously
Finally consider yet another orthonormal, complete,
unsymmetrized,butotherwiseatthisstagearbitraryba-
sis {ζq}, with q|ζqihζq|=ˆI One can write hζp|ζqihζq|ζpi≡ (8πξV2N)3N/2e−2ξ2p2/h¯2, (2.13)
∞ PzN
Ξ± = (±1)phζ |ζ ihζ |e−βHˆ|ζ i whichproductwillshortlyappearinthesummandofthe
N! p q q Pˆp partition function.
NX=0 Xq,pXPˆ
∞ zN
= (±1)phζ |ζ ihζ |e−βHˆ|ζ i.
N! Pˆp q q p
NX=0 Xq,pXPˆ 5. Non-Commutativity Factor
(2.9)
The position and momentum operators do not com-
(One could replace here q ⇒ p′, in which case the final mute, which fundamental property can be accounted for
formis just the non-diagonalsum overa single basisset, by recasting the partition function to include the so-
analogoustothediagonalsumovertheentropybasisset.) called ‘non-commutativity factor’.
Modifying slightly an argument due to Wigner,9 one
can formally commute the Maxwell-Boltzmann operator
4. Momentum and Position States with the momentum basis function by writing12
Now particularize the analysis to momentum and po- e−βHˆζ (r) = ζ (r)e−βHˆ˜1
sition basis functions.12 The momentum basis functions p p
in the position representation r are plane waves, ≡ ζp(r)e−βH(r,p)ew(r,p). (2.14)
ζ (r) = 1 e−p·r/ih¯ The modified energy operator induced here is
p
VN/2
N e−pj·rj/ih¯ Hˆ˜ ≡ ep·r/ih¯Hˆe−p·r/ih¯
= j=1 V1/2 . (2.10) = ep·r/ih¯ −¯h2∇2+U(r) e−p·r/ih¯
Y 2m
(cid:20) (cid:21)
Here p is the configuration momentum. These form a p2 i¯h ¯h2
complete orthonormalset. Periodicboundaryconditions = +U(r)− p·∇− ∇2. (2.15)
2m m 2m
give the width of the momentum state per particle per
dimension as ∆p =2π¯h/V1/3. The so-called non-commutativity function w(r,p) is
The position basis functions are Gaussians, extensive.12 (This is one reason that the present for-
mulation is preferable to those given by Wigner9 or by
e−(r−q)2/4ξ2 Kirkwood.10 The other reason is accounting for wave
ζ (r) =
q (2πξ2)3N/4 functionsymmetrization,whichrigorouslycountsforbid-
den and multiple states correctly.) A recursion relation
= N e−(rj−qj)2/4ξ2. (2.11) leadingtoanexpansionforitwillbegivenin§IIBbelow.
(2πξ2)3/4 To leading order it vanishes, w(h¯ =0)=w(β =0)=0.
j=1
Y
With this non-commutativity function, the expecta-
Here q is the configuration position. Normalization tion value of the Maxwell-Boltzmann operator that ap-
fixes the spacing of the configuration position states pearsintheposition-momentumspaceformulationofthe
4
partition function becomes 7. Monomer Grand Potential
hζ |e−βHˆ|ζ i = hζ |e−βHˆζ i
q p q p
Shortly it will be shown that the symmetrization fac-
= ζ (r)|ζ (r)e−βH(r,p)ew(r,p) tor, η±, breaks up into loops of particles, l = 1,2,...,
q p
which may be called monomers,dimers etc. The leading
= hDζq|ζpie−βH(q,p)ew(q,p). E(2.16) monomer term is unity, η(1) = 1. Hence the monomer
The final equality follows because lim ζ (r) = δ(q− grand partition function is
ξ→0 q
r).
∞ zN
Ξ = dΓe−βH(Γ)ew(Γ). (2.20)
1 N!h3N
N=0 Z
6. Partition Function X
This is the same for bosons as for fermions. The ratio of
The obvious merit of formulating the grand partition thefullpartitionfunctiontothemonomerpartitionfunc-
function in terms of the asymmetric expectation value tion is just the monomer average of the symmetrization
of the Maxwell-Boltzmannoperatoris that the sum over factor,
entropy states has become a sum over points in classical
Ξ± 1 ∞ zN
phase space. The continuum limit of this is = dΓe−βH(Γ)ew(Γ)η±(Γ)
Ξ Ξ N!h3N
1 1 1 N=0 Z
⇒ dΓ. (2.17) X
(∆ ∆ )3N = η± . (2.21)
q,p p q Z 1
X
The volume elements are ∆ = 2π¯h/V1/3 and ∆ = Since the(cid:10)lead(cid:11)ing order of the non-commutativity fac-
p q
8πξ2.12 torvanishes, w(h¯ =0)=w(β =0)=0(see§IIBbelow),
the leading order of the monomer grandpotential is just
This way of representing quantum mechanics in clas-
p the classical grand potential,
sical phase space is distinctly different to the way ad-
vocated by Wigner,9 Kirkwood,10 and followers.18–25 In ∞ zN
particular, previous authors focus upon a transform of Ξ ≡Ξ (w =0)= dΓe−βH(Γ). (2.22)
1,0 1 N!h3N
thewavefunctionratherthanthestatesthemselves,they
N=0 Z
X
do notdistinguishthe positionrepresentationr fromthe
positionconfigurationqashere,andtheydonotaccount The logarithm of this gives the classical equilib-
rium grand potential, Ω = −k T lnΞ . Following
for wave function symmetrization. Finally, and again in 1,0 B 1,0
the usual procedures of classical statistical mechanics,
contrasttothepresent,theirfocusisneitherthepartition
the difference between the quantum and the classical
function nor quantum statistical mechanics.
monomer grand potentials (ie. the quantum correction
In view of the above the partition function becomes
due to monomers)is just a classicalequilibrium average,
∞ zN
Ξ± = N! (±1)phζPˆp|φqihζq|e−βHˆ|ζpi e−β[Ω1−Ω1,0] = Ξ1 =hewi , (2.23)
NX=0 Xq,pXPˆ Ξ1,0 1,0
∞ zN
= dΓ where Ω is the full monomer grandpotential. Here and
N!(∆ ∆ )3N 1
N=0 p q Z below the subscript 1,0 denotes the classical equilibrium
X
× (±1)phζ |ζ ihζ |ζ ie−βH(q,p)ew(q,p) average or classical equilibrium thermodynamic poten-
Pˆp q q p tial.
XPˆ
∞ zN
= dΓ
N!(∆ ∆ )3N 8. Permutation Loop Expansion of the Grand Potential
N=0 p q Z
X
×η±(q,p)hζ |ζ ihζ |ζ ie−βH(q,p)ew(q,p)
p q q p The permutation operator breaks up into loops
∞ zN
= N!h3N dΓe−βH(Γ)ew(Γ)η±(Γ). (2.18) (±1)p Pˆ = ˆI± ′ Pˆij + ′ PˆijPˆjk
N=0 Z
X XPˆ Xi,j iX,j,k
Noticethatthefactorsofξ cancel,whichallowsthelimit
+ ′ Pˆ Pˆ ±... (2.24)
ξ →0 to be taken. Here ij kl
i,j,k,l
1 X
η±(Γ) ≡ (±1)phζ |ζ i, (2.19)
hζ |ζ i Pˆp q The prime on the sums restrict them to unique loops,
p q
XPˆ with each index being different. The first term is just
which is just the complex conjugate of the function the identity. The second term is a dimer loop, the third
χ˜±(Γ) defined previously.12 This might be called the termis atrimerloop,andthe fourthtermisthe product
symmetrization factor. of two different dimers.
5
The symmetrization factor, η±(Γ) = (±1)p theratioofthefulltothemonomerpartitionfunctionas
Pˆ
hζ |ζ i/hζ |ζ i, is the sum of the expectation values of
Pˆp q p q P
these loops. Since the expectation value is the product Ξ±
= η±
of the individual particle factors, Eq. (2.12), the factors Ξ 1
1
due to the unpermutedparticlesinthe numeratorcancel
(cid:10) (cid:11)
with those in the denominator. = 1+ ′η±(2) + ′η±(3)
ij ijk
The monomer symmetrization factor comes from just * + * +
Xij 1 Xijk 1
the unpermuted expectation value,
+ ′η±(2)η±(2) +...
ij kl
hζ |ζ i * +
η(1)(Γ)≡ p q =1. (2.25) Xijkl 1
hζp|ζqi N N!
= 1+ η±(2) + η±(3)
(N −2)!2 (N −3)!3
(cid:28) (cid:29)1 (cid:28) (cid:29)1
The dimer overlap factor in the microstate Γ for par- 2
1 N!
ticles j and k is + η±(2) +...
2 (N −2)!2
(cid:28) (cid:29)1
±hζ |ζ i 1 ∞ N! ml
ηj±k(2)(Γ) = hζPˆpj|kζpqiq = {Xml}ml!Yl=2(cid:28)(N −l)!lη±(l)(cid:29)1
= ±hζpk|ζqjihζpj|ζqki = ∞ ∞ 1 N! η±(l) ml
= ±heζ(pqjk|−ζqqjji)·hpζjp/kih¯|ζeq(kqij−qk)·pk/ih¯. (2.26) Yl∞=2mXl=0ml!(cid:28)(N −l)!l (cid:29)1
N!
= exp η±(l) . (2.29)
(N −l)!l
The important point is the formally exact factorization l=2 (cid:28) (cid:29)1
Y
of the expectation value, leaving only the permuted par-
ticles to contribute. The third and following equalities write the average of
the product as the product of the averages. This is valid
The symmetrization factors are localized in the sense
inthe thermodynamic limit, since the product ofthe av-
that they are only non-zero when all the particles of
erage of two loops scales as V2, whereas the correlated
the loop are close together. (More precisely, localization
interaction of two loops scales as V. The combinatorial
means that the separations between consecutive neigh-
factor accounts for the number of unique loops in each
bors around the loop are all small.) This will be shown
term; η±(l) without subscripts refers to any one set of l
explicitlybelow,buthereitcanbenotedthattheexplicit
particles, since all sets give the same average.
exponents give highly oscillatoryand thereforecanceling
behaviorunless the differences in configurationpositions The grandpotential is essentially the logarithmof the
are all close to zero. grand partition function, Ω≡−kBT lnΞ. Hence the dif-
ferencebetweenthefullgrandpotentialandthemonomer
Similarlythetrimersymmetrizationfactorforparticles
grand potential is just the series of loop potentials,
j, k, and l is
Ξ±
ηj±k(l3)(Γ) = hζPˆhjζkPˆ|kζlpi|ζqi (2.27) −β[Ω±−Ω1] = ln Ξ1
p q ∞
N!
hζ |ζ ihζ |ζ ihζ |ζ i = η±(l)
= pk qj pj ql pl qk (N −l)!l
hζ |ζ ihζ |ζ ihζ |ζ i l=2(cid:28) (cid:29)1
pj qj pk qk pl ql X ∞
= e(qj−qk)·pk/ih¯e(qk−ql)·pl/ih¯e(ql−qj)·pj/ih¯. ≡ −β Ω±. (2.30)
l
l=2
X
Continuing in this fashion, the symmetrization factor
can be written as a series of loop products, The monomer grand potential is of course Ω1 ≡
−k T lnΞ , with the monomer grand partition function
B 1
being given by Eq. (2.20).
η±(Γ) = 1+ ′η±(2)(Γ)+ ′η±(3)(Γ)
ij ijk
ij ijk
X X
+ ′η±(2)(Γ)η±(2)(Γ)+... (2.28)
ij kl B. Expansion of the Non-Commutativity Function
ijkl
X
Usinga similarapproachto Kirkwood,10 the tempera-
Here the superscript is the order of the loop, and the
subscripts arethe atoms involvedin the loop. This gives ture derivative of the defining equation (2.14), e−βHˆ˜1 =
6
e−βHew, may be rearrangedto give12 C. Monomer Expansion A
∂w i¯h ¯h2
= eβU−wp·∇ ew−βU + eβU−w∇2 ew−βU The quantum correction to the classical grand poten-
∂β m 2m
tialduetothemonomersisjustaclassicalaverageofthe
i¯h (cid:8) ¯h2(cid:9) (cid:8) (cid:9) quantum weight due to non-commutativity, Eq. (2.23).
= p·∇(w−βU)+ ∇(w−βU)·∇(w−βU)
m 2m To fourth order in ¯h this is
n
+∇2(w−βU) . (2.31)
−β[Ω −Ω ] (2.38)
1 1,0
Write o = lnhewi1,0
∞ 1 1 1
w ≡ wn¯hn, (2.32) = ln 1+w+ 2!w2+ 3!w3+ 4!w4
(cid:28) (cid:29)1,0
n=1
X 1 1 1
with w (β = 0) = 0. (This begins at n = 1 because the = w+ w2+ w3+ w4
n 2 3! 4!
classicalpartmustyieldtheclassicalMaxwell-Boltzmann (cid:28) (cid:29)1,0
factor, w(h¯ = 0) = 0.) This expansion leads to the re- 1 1 2
− w+ w2
cursion relation for n>2, 2 2
(cid:28) (cid:29)1,0
∂∂wβn = mi p·∇wn−1+ 21mn−2∇wn−2−j ·∇wj = ¯h2hw2i1,0+ ¯h22 w12 1,0
− β ∇wn−2·∇U +Xj=01 ∇2wn−2. (2.33) + ¯h24 w22−hw(cid:10)2i21,(cid:11)0 1,0+ ¯h44! w14 1,0− ¯h84 w12 21,0
m 2m +¯h4hDw i +¯h4hwEw i (cid:10) (cid:11) (cid:10) (cid:11)
4 1,0 1 3 1,0
It is straightforward to derive the first several coeffi-
cient functions explicitly. One has for n=1, + ¯h4 w2w − ¯h4 w2 hw i .
2 1 2 1,0 2 1 1,0 2 1,0
−iβ2
w = p·∇U, (2.34) (cid:10) (cid:11) (cid:10) (cid:11)
1 Allhigherordercontributionshavebeensettozerohere.
2m
It is necessary to write out the individual terms explic-
for n=2,
itly to cancel the odd imaginary terms. It is w (p,q)
j
β3 1 β3 β2 whose average is taken. This is extensive; the factor of
w2 = 6m2pp:∇∇U + 2m 3 ∇U ·∇U − 2 ∇2U , 3 between the terms involving w14 is a necessary part of
(cid:26) (cid:27) the cancelation that accounts for the three ways of pair-
(2.35)
ing the four momenta, as can be shown analytically and
for n=3, numerically (see, for example, §IID4 below). This ex-
pansion terminated at O(h¯2) may be called A2, and at
iβ4 .. 5iβ4 O(h¯4) it may be called A4.
w = ppp.∇∇∇U + p(∇U):∇∇U
3 24m3 24m2
iβ3
− p·∇∇2U. (2.36) D. Monomer Expansion B
6m2
and for n=4,
One can define the cumulative weight as
w4 = −5!im4β45(p·∇)4U − 30βm53(∇U)pp...∇∇∇U w(n) =h¯w1+¯h2w2+...+¯hnwn, (2.39)
β5 β4
− (∇U)(∇U):∇∇U + ∇U ·∇∇2U andthenthapproximationtothequantumcorrectionto
15m2 16m2 the monomer grand potential as
β4 β4
+ pp:∇∇∇2U + ∇2(∇U ·∇U)
48m3 48m2 ∆Ω(n) =−k T ln ew(n) . (2.40)
1 B
β3 β5 1,0
− ∇2∇2U − (p·∇∇U)·(p·∇∇U). D E
24m2 40m3 One has lim ∆Ω(n) =Ω −Ω .
(2.37) n→∞ 1 1 1,0
Noticethattheoddcoefficientfunctionsarepureimag-
1. Momentum Averages
inary and odd in momentum. Because H(Γ) is an even
function of momentum, the quantum weight ew is real
andoscillatoryonaverage. (Theη±(l) containtermsthat The classical average h...i includes an average over
1,0
areeitherrealandeven,orelseimaginaryandoddinmo- the momenta h...i as well as one over the position
1,0,p
mentum.) configurations.
7
Theclassicalmonomerprobabilitydistributionformo- This leaves the position configurationterms to be added
mentum is a Gaussian, as
℘(p)=[2πmkBT]−3N/2e−βp2/2m. (2.41) w(2) = ¯h2β3∇U ·∇U − ¯h2β2∇2U. (2.47)
q 6m 4m
Hence
It is simplest to deal with the Gaussian exponent di-
hppi =mk T I. (2.42) rectly by completing the squares. Writing U′ ≡∇U and
1,0,p B U′′ ≡∇∇U, the momentum part of the exponent is
The thermal wave length is
−β
p2+w(2)
Λ=[2π¯h2/mk T]1/2. (2.43) 2m p
B
−β i¯hβ2 ¯h2β3
= p2− p·U′+ pp:U′′
2m 2m 6m2
2. First order term
−β 2mi¯hβ2 2
= A: p+ A−1U′
2m 2β 2m
The momentum averageof the first order term is (cid:20) (cid:21)
¯h2β3
ew(1) = eh¯w1 − 8m A−1 :U′U′ (2.48)
1,0,p
1,0,p
D E (cid:10)∞ 1(cid:11) −i¯hβ2 n where A≡I−¯h2β2U′′/3m. One can write
= h[p·∇U]ni
n! 2m 1,0,p
n=0 (cid:18) (cid:19) ¯h2β2
X A−1 =I+ U′′+O(h¯4), (2.49)
∞ 1 −i¯hβ2 2l (2l)! 3m
=
(2l)! 2m 2ll! Since the sum of the squares of the eigenvalues is the
l=0 (cid:18) (cid:19)
X trace of the square of the matrix, one has
×(mk T)l[(∇U)·(∇U)]l
B
= ∞ (−1)l ¯h2β l β2(∇U)·(∇U) l ln|A|1/2 = 12 ln[1+λj] (2.50)
l! 8m j
l=0 (cid:18) (cid:19) X
X (cid:2) (cid:3) 1 1
−β2Λ2 = [λ − λ2+...]
= exp (∇U)·(∇U) . (2.44) 2 j 2 j
16π j
(cid:26) (cid:27) X
−¯h2β2 1¯h4β4
Terms with an odd number of momenta vanish upon av- = ∇2U − (∇∇U):(∇∇U).
6m 4 9m2
eraging. The combinatorial factor in the third equality
arisesfromthenumberofwaysofarrangingthemomenta With these the momentum averageis
in pairs: there are (2l)! ways of arranging the momenta,
of which interchanging the l! pairs leaves the result un- ewp(2)
changed, as does the 2l arrangements that arise from 1,0,p
D E
swapping the first and second members of each of the 2πmk TA−1 1/2 −¯h2β3
pairs. Inthe thirdequalitythe thermodynamiclimit has = B exp A−1 :U′U′
[2πmk T]3N/2 8m
been assumed, which means that N −1 has been taken (cid:12) B (cid:12) (cid:26) (cid:27)
(cid:12) (cid:12)
to equal N, etc. The final result is nice because it shows −¯h2β3
= exp ∇U ·∇U
thatthegradientofthepotentialisGaussiandistributed, 8m
(cid:26)
whichisto saythatconfigurationswithverylargegradi- ¯h2β2 ¯h4β4
ents are suppressed. + ∇2U + (∇∇U):(∇∇U)
6m 36m2
This gives the first order result
¯h4β5
− ∇U∇U :∇∇U +O(h¯6) . (2.51)
∆Ω(1) =−k T ln e−β2Λ2(∇U)·(∇U)/16π . (2.45) 24m2
1 B (cid:27)
1,0
D E Thisexponentcanbeaddedtow(2)andtheconfiguration
q
position average taken to obtain ∆Ω(2) to O(h¯6). This
3. Second Order Order Term 1
agrees with expansion A to O(h¯2).
In the second order approximation to the weight,
w(2) = ¯hw1 + h¯2w2, the terms that depend upon the 4. Terms O(¯h4)
momenta may be identified,
wp(2) = −2i¯hmβ2p·∇U + ¯h62mβ23pp:∇∇U. (2.46) conIntriobrudteirontofrogmettthheedtyeardmicstOri(ph¯le4t)popnpe cnaenedbsewob3.taiTnhede
8
by setting one of the momenta to its most likely value, position configurations to O(h¯4) is
p= −i¯hβA−1U′ = −i¯hβ∇U +O(h¯3), (2.52) w(4) (2.56)
2 2 −¯h2β3
= ∇U ·∇U
and the remaining two factors equal to their average 8m
¯h2β2 ¯h4β4
+ ∇2U + (∇∇U):(∇∇U)
hppi =mk TA−1 =mk TI+O(h¯2). (2.53) 6m 36m2
1,0,p B B
¯h4β5
− ∇U∇U :∇∇U
There are three ways of doing this, and so w3 becomes 24m2
¯h2β3 ¯h2β2
¯h3w = i¯h3β4ppp...∇∇∇U + 5i¯h3β4p(∇U):∇∇U + 6m ∇U ·∇U − 4m ∇2U
3 24m3 24m2 ¯h4β4 5h¯4β5
i¯h3β3 − (∇U)·∇∇2U + (∇U)(∇U):∇∇U
− p·∇∇2U 48m2 48m2
6m2 11¯h4β3 17¯h4β4
i¯h3β4−i¯hβ − ∇2∇2U + (∇U)·∇∇2U
= 3 mk T(∇U)·∇∇2U 240m2 240m2
24m3 2 B ¯h4β5 ¯h4β4
5i¯h3β4−i¯hβ − (∇U)(∇U):∇∇U + (∇∇U):(∇∇U)
+ (∇U)(∇U):∇∇U 15m2 60m2
24m2 2 ¯h2β3 ¯h2β2
i¯h3β3−i¯hβ = ∇U ·∇U − ∇2U
− (∇U)·∇∇2U 24m 12m
6m2 2 2h¯4β4 ¯h4β5
−¯h4β4 + (∇∇U):(∇∇U)− ∇U∇U :∇∇U
= (∇U)·∇∇2U 45m2 240m2
48m2 ¯h4β4 11¯h4β3
5h¯4β5 + (∇U)·∇∇2U − ∇2∇2U +O(h¯6).
+ (∇U)(∇U):∇∇U +O(h¯6). (2.54) 20m2 240m2
48m2
This expansion terminated at O(h¯2) may be called B2,
One can do a similar trick for w4, and at O(h¯4) it may be called B4.
w = −i4β5(p·∇)4U − β5 (∇U)pp...∇∇∇U
4 5!m4 30m3 E. Monomer Expansion C
β5 β4
− (∇U)(∇U):∇∇U + ∇U ·∇∇2U
15m2 16m2 Expansion A has the merit of requiring less analysis,
β4 β4 and of requiring the classical average of extensive terms
+ pp:∇∇∇2U + ∇2(∇U ·∇U)
48m3 48m2 and their products. It has the disadvantage of requiring
β3 β5 the average to be taken numerically over the momen-
− ∇2∇2U − (p·∇∇U)·(p·∇∇U) tum configurations. Also, as a fluctuation expression, it
24m2 40m3
requires relatively high accuracy in the individual aver-
−3β3 β4
= ∇2∇2U − (∇U)·∇∇2U agesofproductstogetthenecessarycancelationbetween
5!m2 30m2
these super-extensive terms to end up with the final ex-
β5 β4
− (∇U)(∇U):∇∇U + ∇U ·∇∇2U tensive result for Ω1−Ω1,0. The larger the system size,
15m2 16m2 the greaterthe number ofconfigurationsthat needto be
β3 β4 generated to get acceptable statistical accuracy.
+ ∇2∇2U + ∇2(∇U ·∇U)
48m3 48m2 Expression B has the advantage of not requiring nu-
β3 β4 merical momentum averaging, which leads to higher ac-
− ∇2∇2U − (∇∇U):(∇∇U)
24m2 40m2 curacybecausecertainessentialcancelationsareenforced
−11β3 17β4 analytically. (In most of the results reported below, 32
= ∇2∇2U + (∇U)·∇∇2U momentum configurations were used for each position
240m2 240m2
configurationforthe averageofexpansionA;nomomen-
β5
− (∇U)(∇U):∇∇U tum configurations are required for expansion B.) Also,
15m2 by exponentiating the expansion, as in B, one expects
β4 faster convergence than for an expansion of the final ex-
+ (∇∇U):(∇∇U). (2.55)
60m2 pression, as in A. A disadvantage of expression B has is
that it requires more explicit algebra. Also, more prob-
Thisismultipliedby¯h4. Theneglectedcontributionsare lematic, since w is an extensive variable, taking the av-
of higher order in ¯h (ie. O(h¯6)). erageof ew canleadto computationaloverflowproblems
Insummary,theexponentthatonehastoaverageover for large systems.
9
One wayto avoidnumericaloverflowin evaluatingthe One sees that the analysis of the momentum averagesin
exponent in expansion B, but to preserve the advantage §IID1 goesthroughwiththe only changebeing effective
of no numerical momentum average, is to use the mo- potential gradients for particles j and k,
mentumaveragedw(4) intheexpansionA.Thatis,write
2m
Eq. (2.56) as explicit powers of Planck’s constant, ∇ U˜ ≡ ∇ U − q ,
j j ¯h2β2 jk
w(4) ≡¯h2w˜2+¯h4w˜4. (2.57) and ∇ U˜ ≡ ∇ U − 2m q , (2.62)
k k ¯h2β2 kj
One can now insert this into expansion A, Eq. (2.38),
with w1 =w3 =0, respectively. From the zeroth order dimer result, qjk ∼
Λ∼¯h.
−β[Ω −Ω ] (2.58) With these effective gradients, the momentum part of
1 1,0
the exponent given by Eq. (2.48) becomes
= ln ew(4)
1,0 −β
= ¯h2Dhw˜ i E+ ¯h4 w˜2−hw˜ i2 +¯h4hw˜ i . 2mp2+w˜p(2)
2 1,0 2 2 2 1,0 1,0 4 1,0 −β 2mi¯hβ2 2
D E = A: p+ A−1∇U˜
2m 2β 2m
This may be called expansion C4. The averages here (cid:20) (cid:21)
are classical averages over configuration positions. The ¯h2β3
neglected terms are O(h¯6). Hence this can be expected − 8m A−1 :∇U˜∇U˜, (2.63)
to be rather similar to expansion A4, the only difference
with the matrix A unchanged. Hence the average over
beingthatthemomentumaverageshavebeenperformed
the momenta, Eq. (2.51), becomes
analytically to leading order before expanding the expo-
nential.
ew˜p(2)
1,0,p
D E
−¯h2β3
F. Dimer = exp ∇U˜ ·∇U˜
8m
(cid:26)
¯h2β2 ¯h4β4
In general the loop potential l ≥ 2 is given by + ∇2U + (∇∇U):(∇∇U)
Eq. (2.30), 6m 36m2
¯h4β5
− ∇U˜∇U˜ :∇∇U . (2.64)
N! 24m2
−βΩ± ≡ η±(l) (cid:27)
l (N −l)!l
(cid:28) (cid:29)1 This expression is half-way between B2 and B4. Ar-
−1 N!
= ew η±(l)ew .(2.59) guably, it would be most consistent to neglect here the
D E1,0 (cid:28)(N −l)!l (cid:29)1,0 terms O(h¯4).
To this should be added the position configuration
The zeroth order approximation for non-commutativity terms,
sets w = 0 and is the one for which numerical results
are discussed below. However the case w = w(2) is now ¯h2β3 ¯h2β2
w(2) = ∇U ·∇U − ∇2U. (2.65)
analyzed. q 6m 4m
For the dimer l =2, Eq. (2.26) is
The zerothorder dimer resultis in the very firstterm,
(−¯h2β3/8m)∇U˜·∇U˜,whichuponsettingU =0becomes
N!
η±(l)(ΓN) = ± e−qjk·pjk/ih¯. (2.60) −mq2 /¯h2β =−2πq2 /Λ2.
(N −2)!2 jk jk
j<k
X
Combining this exponent with the partof ¯hw (ΓN) that III. METHODOLOGY AND RESULTS
1
depends on p and p gives exponent
j k
A. Simulation Algorithm
¯hw˜ (p ,p ;qN)
1 j k
−1 i¯hβ2 Computer simulations were performed for a classical
= pjk·qjk− {pj ·∇jU +pk·∇kU} canonical equilibrium system in configuration position
i¯h 2m
space.8TheMetropolisalgorithmwasused,withthestep
−i¯hβ2 2m
= p · ∇ U − q length adjusted to give approximately a 50% acceptance
2m j j ¯h2β2 jk
(cid:26) (cid:20) (cid:21) rate.
2m Asmallcell,spatiallybased,linkedlist,neighbortable
+p · ∇ U − q . (2.61)
k k ¯h2β2 kj wasused.26 Thecellshadedgelengthontheorderof3/4
(cid:20) (cid:21)(cid:27)
10
of the molecular diameter, which means that most cells
were empty, and almost no cells were occupied by more TABLE I: Lennard-Jones Parameters.2
than one particle. The total volume occupied by neigh-
bor cellswasonthe orderof2–5times the volumeofthe σ rmin ǫ/kB m
(nm) (nm) (K) (amu)
potentialcut-offsphere. Testswith250–2000atomscon-
firmedthatthetime costofthe algorithmscaledlinearly
He 0.2556 0.2869 10.22 4.002602
with the system size.
Ne 0.2789 0.3131 35.7 20.1797
A cycle consists of an attempted move for every par- Ar 0.3418 0.3837 124.0 39.948
ticle in turn. After every 20 cycles, the position config-
uration was stored for use in later averaging. For most
of the results reportedbelow, 20,000configurationswere
saved. For each temperature and density, a number of The Lennard-Jones parameters for the noble elements
equilibration cycles were run, starting from a previously used in the simulations are shown in Table I.
equilibrated configuration at a nearby state point. After For the Lennard-Jones fluid, the critical temperature,
equilibration was judged complete, generally on the ba- density, and pressure in dimensionless form are?
sis of the stability of the virial pressure, the simulation
k T p σ3
proper commenced. B c =1.35, ρ σ3 =0.35, and c =0.142, (3.2)
c
ǫ k T
For averaging, the configurations were broken into 50 B c
blocks, and the statistical error was estimated from the respectively.
fluctuations of the average of each block. The error
reported below corresponds to one standard deviation,
which corresponds to 68% confidence in the result. C. Pair Potential
The expansion A requires averaging over the momen-
tum configurations in addition to averagingover the po- For the case of a general pair potential, including the
sition configurations that have been saved. (The expan- present Lennard-Jones potential, the potential energy is
sions B and C do not require this, but these were calcu-
lated at the same time as A.) For this purpose, for each U(q)= u(q ), (3.3)
jk
position configuration, a number, generally 32, of inde-
j<k
X
pendentmomentumconfigurationsweregenerated. Each
momentum configuration corresponds to a new momen- with qjk =qj −qk and qjk =|qjk|.
tum for every atom. These were drawn randomly from The gradient is
a Gaussian distribution that corresponds to the classical
N
q
Maxwell-Boltzmann distribution of the kinetic energy. {∇U} = (k6=j)u′(q ) jk,α, (3.4)
j,α jk
The position configurations do not depend upon the qjk
k=1
mass of the atoms, and the same position configurations X
were used for all three noble elements at that particular where α=x,y, or z. Hence
(dimensionless) state point. The momenta do depend p ·q
upon the atomic mass, and the averaging program was p·∇U = (k6=j)u′(qjk) j jk
q
jk
run three times at each state point, once for each noble j k
XX
p ·q
element. = u′(q ) jk jk. (3.5)
jk
The Lennard-Jones pair potential was used (see next) qjk
j<k
andcutoffatR =3.5σ. Periodicboundaryconditions X
cut
wereenforced,andtheedgelengthofthecubicsimulation The dyadic second derivative is
cell was L > 2R . For the case of the virial pressure
cut ∂ q
and the internal energy (not reported), a tail correction {∇∇U} = (l6=j)u′(q ) jl;α
jα;kγ jl
∂q q
was included in the reported average. For the case of kγ jl
l
X
the quantum correction to quadratic order in ¯h, a tail
q q
correction was used for earlier results.12 However, the = δ (l6=j) u′′(q ) jl;α jl;γ
jk jl q2
quantum corrections reported here do not include a tail l ( jl
X
correction.
δ q q
+u′(q ) αγ −u′(q ) jl;α jl;γ
jl q jl q3
jl jl )
B. Lennard-Jones Potential
q q
+[1−δ ] u′′(q ) jk;α kj;γ
jk jk
q q
The Lennard-Jones potential used here is (cid:26) jk jk
δ q q
σ 12 σ 6 −u′(q ) αγ −u′(q ) jk;α kj;γ ,
4ǫ − , r <R jk q jk q3
u(r)= r r cut (3.1) jk jk )
0, (cid:20)(cid:16) (cid:17) (cid:16) (cid:17) (cid:21) r >R . (3.6)
cut
