Quantum Monte Carlo in Classical Phase Space. Mean Field and Exact Results for a One Dimensional Harmonic Crystal
aa r X i v : . [ qu a n t - ph ] A p r Quantum Monte Carlo in Classical Phase Space. Mean Field and Exact Results for aOne Dimensional Harmonic Crystal.
Phil Attard [email protected]/QSM19/mfxtl/mfxtl.tex
Begun 25-Mar-2019
Monte Carlo simulations are performed in classical phase space for a one-dimensional quantumharmonic crystal. Symmetrization effects for spinless bosons and fermions are quantified. Thealgorithm is tested for a range of parameters against exact results that use 20,000 energy levels.It is shown that the singlet mean field approximation is very accurate at high temperatures, andthat the pair mean field approximation gives a systematic improvement in the intermediate and lowtemperature regime. The latter is derived from a cluster mean field approximation that accountsfor the non-commutativity of position and momentum, and that can be applied in three dimensions.
I. INTRODUCTION
The impediment to applying quantum statistical me-chanics to condensed matter systems at terrestrial tem-peratures and densities is the horrendous scaling withsystem size entailed by conventional approaches.
Al-though the field is not exactly moribund, it would be fairto say that the rate of progress has been disappointing,and that there appears no obvious way of overcoming thisfundamental limitation of the existing methods.The pressing need to approach the problem from a dif-ferent direction has motivated the author to develop anew methodology that is based on a formally exact trans-formation that expresses the quantum partition functionas an integral over classical phase space.
The valid-ity of the formulation has been verified analytically andnumerically for the quantum ideal gas, and for non-interacting quantum harmonic oscillators. Application of the algorithm to interacting systems in-dicate that, for a given statistical error, the algorithmscales sub-linearly with system size.
This places itin a class of its own for the treatment of condensed mat-ter quantum systems.Validation of the algorithm for an interacting sys-tem has been provided by quantitative tests for inter-acting Lennard-Jones particles against benchmarks ob-tained with more conventional methods by Hernando andVan´ıˇcek. The latter results were exact, apart from thefact that they were numerical and ‘only’ 50 energy levelswere used for the statistical averages. One should not un-derstate the difficulty of obtaining numerically 50 energylevels for the Lennard-Jones system; it is undoubtedly agreater computational challenge than obtaining analyti-cally the 20,000 energy levels used below for the presentone-dimensional harmonic crystal. If nothing else thefigure of 50 levels does underscore the insurmountableintractability of conventional quantum approaches.In fact the original motivation for the present paperwas to sort out a small discrepancy between the putativeexact results of Hernando and Van´ıˇcek and the author’smean field, classical phase space results at the highesttemperature studied. It was unclear whether the differ-ences in the original comparison were due to the mean field approximation used in the phase space simulation,or else to the limited number of energy levels used in theexact results. Accordingly the author has undertaken toestablish his own exact results for use as benchmarks,and it was found that 20,000 energy levels were neces-sary for reliable results at temperatures high enough togive a departure from the ground state.The results of these tests are reported below. It isfound that the mean field approximation as originallyformulated is essentially exact at the highest tempera-tures, but the error can be on the order of 5–10% at in-termediate and low temperatures, depending upon modelparameters. Accordingly, this papers develops a system-atic general improvement to the phase space algorithmthat might be called the cluster mean field approxima-tion. This is here implemented at the singlet level, whichis the original version, and also at the pair level, whichis new. It is found that the pair mean field approxima-tion gives essentially exact results at intermediate tem-peratures. At low temperatures such that the system ispredominantly in the ground state, the singlet and pairmean field approximations are found to be in error by onthe order of 5%.This paper relies upon earlier work, which will not bere-derived here. The reader is referred Ref. 8 for thederivation of quantum statistical mechanics, to Ref. 9 forthe derivation of the phase space formulation (see alsoRef. 14 for some formal details concerning symmetriza-tion of multiparticle states), and to Ref. 14 for the exactphonon analysis of the present one-dimensional harmoniccrystal.
II. ANALYSIS AND MODELA. Phase Space Formulation
The details of the classical phase space formulation ofquantum statistical mechanics can vary with the partic-ular quantity being averaged. For the average energy,which represents a certain class of operators, namelythose that are a linear combination of functions purely ofthe momentum operator or of the position operator, thecanonical equilibrium average can be written D ˆ H E ± N,V,T (2.1)= 1 h dN N ! Z ± Z d Γ e − β H ( Γ ) W p ( Γ ) η ± q ( Γ ) H ( Γ ) , where the partition function is Z ± ( N, V, T ) (2.2)= 1 h dN N ! Z d Γ e − β H ( Γ ) W p ( Γ ) η ± q ( Γ ) . In these N is the number of particles, assumed identi-cal and spinless, V is the volume, T is the temperature, β = 1 /k B T is the inverse temperature, h is Planck’s con-stant, and k B is Boltzmann’s constant. Also, Γ = { p , q } is a point in phase space, with the vector of particles’momenta being p = { p , p , . . . , p N } , and that of theparticles’ position being q = { q , q , . . . , q N } . Finally,ˆ H = H (ˆ p , ˆ q ) is the energy or Hamiltonian operator, and H ( p , q ) is the classical Hamiltonian function.The plus sign is for bosons and the minus sign is forfermions. Actually, since this has been formulated forspinless particles, these refer to the fully symmetrizedand fully anti-symmetrized spatial part of the wave func-tion (see Appendix C of Ref. 14). Throughout the words‘boson’ and ‘fermion’ should be understood in this sense.The unsymmetrized position and momentum eigen-functions in the position representation are respectively | q i = δ ( r − q ) , and | p i = e − p · r /i ¯ h V N/ . (2.3)The symmetrization function is defined as η ± q ( p , q ) ≡ h p | q i X ˆP ( ± p h ˆP p | q i . (2.4)Here the sum is over the N ! permutation operators ˆP,whose parity is p . The imaginary part of these is odd inmomentum, η ± q ( p , q ) ∗ = η ± q ( − p , q ).The symmetrization function can be written as a seriesof loop products, η ± q ( Γ ) = 1 + X ij ′ η ± (2) q ; ij + X ijk ′ η ± (3) q ; ijk + X ijkl ′ η ± (2) q ; ij η ± (2) q ; kl + . . . (2.5)Here the superscript is the order of the loop, and thesubscripts are the atoms involved in the loop. The primesignifies that the sum is over unique loops (ie. each config-uration of particles in loops occurs once only) with eachindex different (ie. no particle may belong to more thanone loop). In general the l -loop symmetrization factor is η ± ( l ) q ;1 ...l = ( ± l − e q l · p l /i ¯ h l − Y j =1 e q j +1 ,j · p j /i ¯ h , (2.6) where q jk ≡ q j − q k . This corrects a typographical errorin Eq. (3.4) of Ref. 9The commutation function, which is essentially thesame as the function introduced by Wigner and an-alyzed by Kirkwood, is defined by e − β H ( p , q ) W p ( p , q ) = h q | e − β ˆ H | p ih q | p i . (2.7)Again one has W p ( p , q ) ∗ = W p ( − p , q ). High tempera-ture expansions for the commutation function have beengiven. The commutation function in phase space can also bewritten as a series of energy eigenfunctions and eigenval-ues, ˆ H| n i = E n | n i . Using the completeness properties ofthese one obtains e − β H ( p , q ) W p ( p , q ) = h q | e − β ˆ H | p ih q | p i (2.8)= 1 h q | p i X n e − βE n h q | n i h n | p i . This exact expression forms the basis of the mean fieldapproximation to the commutation function.
B. Cluster Mean Field Approximation
In the classical phase space formulation of quantumstatistical mechanics, the symmetrization function is rel-atively trivial to obtain and implement. The commuta-tion function is more of a challenge, with the most suc-cessful approach using a mean field approximation thatexploits the analytic form of the commutation functionin the case of independent simple harmonic oscillators. This has previously been tested for the simulation of aLennard-Jones system. This section begins with a summary of the singlet meanfield approximation to the commutation function. Thenthe cluster mean field approximation is given.
1. Singlet Mean Field Approximation
In general, the particles of the sub-system interact viathe potential energy, which is the sum of one-body, two-body, three-body terms, etc., U ( q ) = N X j =1 u (1) ( q j ) + N X j
0, and orthonormal eigenvectors, U ′′ j X jα = λ jα X jα , α = x, y, . . . , d . For molecule j inconfiguration q the eigenvalues define the frequencies ω jα ( q ) = q λ jα ( q ) /m, α = x, y, . . . , d. (2.12)With this the potential energy is U ( q ) = N X j =1 U j + 12 N X j =1 ( q j − q j )( q j − q j ) : U ′′ j = N X j =1 U j + 12 X j,α ¯ hω jα Q jα . (2.13)Here Q jα ≡ p mω jα / ¯ h Q ′ jα , and Q ′ j = X T j [ q j − q j ].(This corrects a typographical error in Eqs (2.7) and (2.8)in Ref. 13.)The mean field approximation combined with the sec-ond order expansion about the local minima maps eachconfiguration Γ to a system of dN independent harmonicoscillators with frequencies ω jα displacements Q jα , andmomenta P jα = n X T j p j o α / p m ¯ hω jα . (This corrects atypographical error in Ref. 14.)With this harmonic approximation for the potentialenergy, the effective Hamiltonian in a particular configu-ration can be written H SHO ( p , q − q ) = N X j =1 U j + 12 X j,α ¯ hω jα (cid:2) P jα + Q jα (cid:3) . (2.14) The commutation function for the interacting systemfor a particular configuration can be approximated asthe product of commutation functions for effective non-interacting harmonic oscillators which have the local dis-placement as their argument. With this the mean fieldcommutation function is W mf p ( Γ ) ≈ W SHO p ( p , q − q )= e β H SHO ( p , q − q ) h q − q | e − β ˆ H SHO | p ih q − q | p i = Y j,α W SHO p,jα ( P jα , Q jα ) . (2.15)The harmonic oscillator commutation function for a sin-gle mode is W SHO p,jα ( P jα , Q jα ) (2.16)= √ e − iP jα Q jα e β ¯ hω jα [ P jα + Q jα ] / e − [ P jα + Q jα ] / × ∞ X n jα =0 i n jα e − β ¯ hω jα ( n jα +1 / n jα n jα ! H n jα ( P jα )H n jα ( Q jα ) . The prefactor e − iP jα Q jα corrects the prefactor e − ip jα q jα / ¯ h given in Eq. (5.10) of Ref. [10]. HereH n ( z ) is the Hermite polynomial of degree n . The imag-inary terms here are odd in momentum. As justified byanalytic results for the simple harmonic oscillator, forconfigurations such that U ′′ j ( q ) is not positive definite(ie. a particular eigenvalue is not positive, λ jα ≤ Q jα exceeds a predeterminedcut-off, the corresponding commutation function can beset to unity, W SHO p,jα = 1.For the averages, the momentum integrals can be per-formed analytically, both here and in combination withthe symmetrization function. This considerably reducescomputer time and substantially increases accuracy.
2. Cluster Mean Field Approximation
Any configuration q can be decomposed into disjointclusters labeled α = 1 , , . . . . Of the different criteriathat can be used to define a cluster, perhaps the simplestis that two particles belong to the same cluster if, andonly if, they are connected by at least one path of bonds.Two particles are bonded if their separation is less thana nominated length. Some clusters, perhaps the greatmajority, will consist of a single particle.An even simpler definition can be made in one di-mension. In this case define the pair cluster α , α =1 , , . . . , N/
2, as the nearest neighbors { α − , α } , irre-spective of their actual separation. This criteria is usedin the results presented below.Using a separation-based criterion for the definition ofa cluster is useful not only for the mean field approx-imation to the commutation function, but also for thecalculation of the symmetrization function. Dependingon the chosen bond length, only permutations of parti-cles within the same cluster need to be considered. (Thisidea is not used in the results presented below.)The cluster energy is the internal energy plus the rele-vant proportion of the interaction energy with other clus-ters: half for pair interactions, one third for triplet inter-actions, etc. The total potential energy is U ( q ) = N X j =1 u (1) ( q j ) + X j 1. Potential Energy Following earlier work, consider a one-dimensionalharmonic crystal in which the particles are attached bylinear springs to each other and to lattice sites. Let thecoordinate of the j th particle be q j , and let its latticeposition (ie. in the lowest energy state) be ˜ q j = j ∆ q . Thetilde signifies that these are ordered. The lattice spacingis also the relaxed inter-particle spring length. There arefixed ‘wall’ particles at q = 0 and q N +1 = ( N + 1)∆ q .Let d j ≡ q j − ˜ q j be the displacement from the latticeposition; for the wall particles, d = d N +1 = 0. Thesystem has over-all number density ρ = ∆ − q .There is an external harmonic potential of spring con-stant κ acting on each particle centered at its latticesite. The inter-particle spring has strength λ and relaxedlength ∆ q . With these the potential energy is U ( q ) = κ N X j =1 [ q j − ˜ q j ] + λ N X j =0 [ q j +1 − q j − ∆ q ] = κ N X j =1 d j + λ N X j =0 [ d j +1 − d j ] . (2.21)The energy eigenfunctions and eigenvectors can be ob-tained explicitly for this model by expressing it in termsof phonon modes. This model potential is not invariant with respect tothe permutation of the positions of the particles, andit therefore violates a fundamental axiom of quantummechanics. This was discussed in detail in Appendix Aof Ref. 14, where the origin, interpretation, and justifica-tion for the potential was given. To those remarks maybe added the fact that the formulation of quantum sta-tistical mechanics in classical phase space is unchangedby a non-symmetric potential (unpublished). 2. Singlet Mean Field The energy of particle j in configuration q is U j ( q j ; q )= κ q j − ˜ q j ] + λ q j − q j +1 + ∆ q ] + λ q j − q j − − ∆ q ] + λ q − ∆ q ] δ j + λ q N − N ∆ q ] δ jN = κ d j + λ (cid:8) [ d j − d j +1 ] + [ d j − d j − ] (cid:9) + λ d j [ δ j + δ jN ] . (2.22)The total potential energy is just U ( q ) = P Nj =1 U j ( q j ; q ),and d j ≡ q j − ˜ q j .The gradient vanishes when d j ( q ) = λ U ′′ j [ d j +1 + d j − ] . (2.23)The second derivative is U ′′ j ≡ ∇ j ∇ j U j ( q j ; q ) = κ + λ + λ δ j + δ jN ] , (2.24)which is independent of the configuration q .The potential energy of particle j in configuration q may be expanded to second order about its local mini-mum at q j ( q ), U j ( q j ; q ) = U j ( q ) + U ′′ j q j − q j ( q )] . (2.25)where U j ( q ) ≡ U j ( q j ( q ); q ). This second order expan-sion for the potential is exact for the present harmoniccrystal. Note that the most likely position of particle j for the current configuration, q j ( q ), is not the same asits lattice position ˜ q j .For each molecule define the frequency ω j = q U ′′ j /m .This is the same for all configurations q . With this thepotential energy is U ( q ) = N X j =1 U j + 12 N X j =1 U ′′ j [ q j − q j ( q )] = N X j =1 U j + 12 X j ¯ hω j Q j . (2.26)Here Q j ≡ p mω j / ¯ h [ q j − q j ]. As above, the momentaare P j = p j / p m ¯ hω j .This is the approximation in the singlet mean field ap-proach: each frequency mode is a displacement of a sin-gle particle. One may now apply Eqs (2.14), (2.15), and(2.16) for the singlet mean field commutation function. 3. Pair Mean Field Approximation For the one dimensional harmonic crystal, define acluster pair α as the nearest neighbors { α − , α } , α = 1 , , . . . , ⌊ N/ ⌋ . For simplicity, assume N even.(Contrariwise, include particle N as a singlet cluster.)The energy of cluster α is U α ( q α ; q ) (2.27)= κ (cid:8) d α − + d α (cid:9) + λ n [ d α − d α − ] + 12 [ d α − − d α − ] + 12 [ d α +1 − d α ] (cid:27) + λ d α − − d α − ] δ α − , + λ d α +1 − d α ] δ α,N , where the displacement is d j = q j − ˜ q j . Note that theinteraction with the wall particles, when present, has tobe counted fully. Note also that d = d N +1 = 0.The second derivative matrix is U ′′ α = (cid:18) κ + λ + λ δ α − , − λ − λ κ + λ + λ δ α,N (cid:19) ≡ − λ (cid:18) K ′ α K ′′ α (cid:19) . (2.28)Using it gives the optimum cluster displacement as (cid:18) d α − d α (cid:19) = − / K ′ α K ′′ α − (cid:18) K ′′ α d α − − d α +1 − d α − + K ′ α d α +1 (cid:19) . (2.29)The eigenvalues of the second derivative matrix (with-out − λ ) are µ ± α = 12 ( K ′ α + K ′′ α ) ± p ( K ′ α − K ′′ α ) + 4 . (2.30)Since the K α are negative, so are the µ ± α .Writing the eigenvectors as u ± α = c ± α (1 , a ± α ), c ± α ≡ / p a ± α , from the eigenvalue equation one obtains a ± α = µ ± α − K ′ α . This is presumably equivalent to a ∓ α = 1 / [ µ ∓ α − K ′′ α ].The mode frequencies are ω ± α = q − λµ ± α /m , which arereal. The mode amplitude is Q α = ω α X T α ∆ d α (2.31)= q mω + α / ¯ h { c + α ∆ d α − + c + α a + α ∆ d α } q mω − α / ¯ h { c − α ∆ d α − + c − α a − α ∆ d α } , where ∆ d α − ≡ d α − − d α − = q α − − q α − , and∆ d α ≡ d α − d α = q α − q α .The contribution to the total energy from cluster α inconfiguration Γ is H α ( Γ ) = 12 m (cid:2) p α − + p α (cid:3) + U α + 12 U ′′ α : [ q α − q α ] [ q α − q α ]= U α + ¯ hω + α (cid:8) P α, + + Q α, + (cid:9) + ¯ hω − α (cid:8) P α, − + Q α, − (cid:9) . (2.32)In the computational implementation of the algorithm,these can be used to replace directly the singlet meanfield terms, P + α ⇒ P α − , P − α ⇒ P α , and similarly for Q ± α and ω ± α . 4. Symmetrization Function Symmetrization consists of a sum over all particle per-mutations. However because of the highly oscillatoryFourier contributions to the loop symmetrization func-tion, Eq. (2.6), only permutations of closely separatedparticles actually contribute to the statistical average.In view of this one can define the permutation length, d m (ˆP) ≡ N X j =1 | j − j ′ | , j ′ ≡ { ˆP j } j . (2.33)One can see that d m = 0 corresponds to the identity per-mutation, d m = 2 corresponds to a single nearest neigh-bor transposition (dimer), and d m = 4 corresponds toeither two distinct nearest neighbor transpositions (dou-ble dimer), or else a single cyclic permutation of threeconsecutive particles (trimer). One expects that the con-tributions to the symmetrized wave function will decreasewith increasing permutation length.Hence one can set an upper limit on the length of thepermutations that are included in the symmetrizationfunction. The numerical results below show that by farthe greatest contribution comes from the identity per-mutation alone, d m = 0. In some cases a measurablechange occurs by including also nearest neighbor per-mutations, d m = 2. Measurable but smaller change oc-curs upon also including permutations of length d m = 4(not shown below). For N particles, the number of per-mutation terms that contribute to the symmetrizationfunction is 1 for d m = 0, 1 + ( N − 1) for d m ≤ 2, and N + ( N − N − / N − 2) for d m ≤ D. Simulation Algorithm The simulation algorithm was as previouslydescribed. Briefly the Metropolis algorithm inposition space was used with the usual classicalMaxwell-Boltzmann weight. The various momentumintegrals were performed analytically. Averages wereevaluated by umbrella sampling using the commutationfunction and symmetrization function as weight. Sincethree versions of the commutation function (unity (ie.classical), singlet mean field, pair mean field), and threeversions of the symmetrization function ( d m = 0 (ie.classical), and d m ≤ 2, for bosons and for fermions)were tested, 9 different averages for each quantity wereobtained simultaneously.Typically enough configurations were generated tomake the statistical error less than 1%, sometimes muchless. In the simulations the time depends on how manyHermite polynomials are used for the commutation func-tion (eight in the results reported below; tests with sixand twelve showed no great effect), and the cut-off for themode amplitude beyond which the commutation functionwas set to unity ( Q cut = 1 in the results below; tests with Q cut = 2 showed no great effect).Interestingly enough, for an accuracy of about 1%, theMonte Carlo algorithm was a factor of about 2 , 000 timesmore efficient (in terms of total computer time) than the A v e r a g e E n e r gy , < H > / ħ w L J b < H > 05 0 1 2 3 4 5 A v e r a g e E n e r gy , Inverse Temperature, b ħ w LJ b ħ w LJ FIG. 1: Average energy versus inverse temperature ( N = 4,∆ q = r e , λ = κ = mω and d m = 0). The solid curve is theexact result using l max = 10 , 000 energy levels. The trianglesare the singlet and the circles are the pair mean field MonteCarlo simulations. The statistical error is less than the size ofthe symbols. The dotted curve is the classical result, h E i cl = N/β . The dashed line is the ground state, E = P Nn =1 ¯ hω n / Inset. Average energy divided by temperature, β h ˆ Hi , versusinverse temperature. quasi-analytic exact phonon method. The main bottle-neck in the latter was the crude numerical quadraturemethod that was used to evaluate the symmetrizationfunction and density profile, and this was exacerbatedby the large number of energy levels that were requiredfor accurate results at higher temperatures. For this par-ticular comparison at β ¯ hω LJ = 2, l max = 20 , 000 energylevels were necessary; grid parameters ∆ Q = 0 . 15 and Q max = 6 gave a quadrature error of about .8%. Thephonon method is more efficient in one respect, namelythat it requires negligible computer time for each addi-tional temperature point; the simulations give results foronly one temperature at a time.The Lennard-Jones frequency used to scale the resultsbelow is ω LJ = 3 . × Hz, the mass is m = 3 . × − kg, the well-depth is ε = 4 . × − J, and theequilibrium separation is r e = 3 . × − m. Theseparameters are appropriate for neon. III. RESULTS Figure 1 shows the average energy as a function ofinverse temperature using exact calculations and clas-sical phase space simulations. At low temperaturesthe phonon modes are in the ground state, E = P Nn =1 ¯ hω n / 2, and at high temperatures the system ap-proaches the classical result, hHi cl = N/β . It can beseen in the main part of the figure that both the singletand pair mean field approximations are in relatively goodagreement with these limiting results and with the exactcalculations over the whole temperature regime shown.The commutation function provides a significant correc-tion to the classical results at high temperatures. In theregime of Fig. 1, symmetrization effects are negligible andonly the d m = 0 calculations are shown.There are two approximations in the exact calcula-tions: the number of energy levels used and the do-main and spacing of the grid used for the numericalquadrature. (I persist in calling these results ‘exact’ be-cause they use explicit analytic expressions for the energyeigenvalues and eigenfunctions.) The quadrature affectsonly the average density profile without symmetrization,and also the average energy and the average density pro-file with symmetrization. Hence the exact calculationsof the average energy in the absence of symmetrizationeffects, d m = 0, as in Fig. 1, are approximate only asregards to the number of energy levels that are retained.The exact calculations in Fig. 1 use 10,000 energy lev-els. These are adequate for low and intermediate temper-atures, β ¯ hω LJ > ∼ . 24, judged in part by comparison withresults using 5,000 energy levels. The exact data beginsto underestimate the classical result for higher tempera-tures than this, and are not shown in Fig. 1. One mightspeculate that the exact quantum result for the averageenergy of the harmonic crystal should approach the clas-sical limit from above. Using fewer energy levels reducesthe domain of inverse temperatures in which the exactcalculations are reliable.Both the singlet and pair mean field simulations arepractically exact at higher temperatures, β ¯ hω LJ < ∼ . β ¯ hω LJ > ∼ 3, theexact results with 10,000 energy levels are practically in-distinguishable from the ground state energy. At lowtemperatures, both mean field classical phase space ap-proximations give a lower energy than the ground state.For example, in the case of Fig. 1, the ground state en-ergy is E / ¯ hω LJ = 3 . β ¯ hω LJ = 10, the singletmean field theory gives h ˆ Hi / ¯ hω LJ = 3 . ± . . ± . . ± . A v e r a g e E n e r gy , < H > / ħ w L J -0.100.10.2 H > + - < H > - ] / ħ w L J A v e r a g e E n e r gy , Inverse Temperature, b ħ w LJ -0.2-0.1 0 1 2 3 4 5 [< H Inverse Temperature, b ħ w LJ FIG. 2: Average energy versus inverse temperature ( N = 4,∆ q = r e , κ = 0, λ/mω = 0 . d m = 0). The solid curveis the exact result using l max = 20 , 000 energy levels. Thetriangles are the singlet and the circles are the pair meanfield Monte Carlo simulations, respectively. The dotted curveis the classical result, h E i cl = N/β , and the dashed line isthe ground state energy. Inset. The boson energy minus thefermion energy using d m ≤ 2. The crosses are Monte Carlosimulations with the classical commutation function, W = 1. mean field approach substantially reduces the error in thesinglet mean field approach. In the absence of exact data,the difference between the pair and singlet mean field re-sults would give a guide to the quantitative accuracy ofthe former.It is worth mentioning that at each temperature theclassical phase space simulations took about five minuteson a desktop personal computer to obtain the quoted ac-curacy. In comparison, it took about 2 days to obtainthe exact results with these energy levels and quadraturegrid, the latter being the time limiting part of the compu-tations. (This is independent of how many temperaturepoints are saved.)Figure 2 shows the average energy for a weakly coupledcrystal. There is no singlet potential, κ = 0, and thenearest neighbor spring constant has been decreased bya factor of 50, λ/mω = 0 . 02. The lattice spacing andspring length is unchanged. The ground state energy is E = 0 . β ¯ hω LJ = 3, the exact energy is 1.3707 (for l max = 20 , . ± . . ± . . ± . d = 0, so no symmetrization effects are included.For inverse temperatures β ¯ hω LJ < ∼ . 8, the mean fieldalgorithms appear more reliable than the exact resultswith l max = 20 , h ˆ Hi + − h ˆ Hi − ] / ¯ hω LJ , d m ≤ D e n s it y , r r e + -r - ] r e D e n s it y , Position, z/r e -0.04-0.02 -1 0 1 2 3 4 5 6 [ r + - FIG. 3: Density profile at β ¯ hω LJ = 2 ( N = 4, ∆ q = r e , κ = 0, λ/mω = 0 . d m = 0). The solid curve is the exactresult using l max = 20 , 000 energy levels, and the symbols arethe classical phase space Monte Carlo simulations, with thecrosses being the classical result, W = 1, the triangles usingthe singlet and the circles using the pair mean field commu-tation function. Inset. The boson density minus fermiondensity using d m ≤ the fully symmetrized spatial energy eigenfunctions, and‘fermions’ means the fully anti-symmetrized spatial en-ergy eigenfunctions. At lower temperatures, β ¯ hω LJ > ∼ . 5, the difference is positive, which means that the en-ergy for bosons is greater than that for fermions. Itcan be seen that the peak difference, which occurs at β ¯ hω LJ = 1 . 9, is about 3% of the actual energy (exactresults). (The error in the numerical quadrature usedfor the exact result is 2% at this temperature for the en-ergy with d m = 0. Hopefully for d m ≤ β ¯ hω LJ > ∼ 2, the singlet mean field ap-proximation overestimates the energy difference, whereasthe pair mean field approximation perhaps halves the er-ror. The classical results, with commutation function W = 1, performs surprisingly well in this low tempera-ture regime. At higher temperatures than this all fourapproaches indicate that the energy difference turns neg-ative. The extent to which the singlet and pair mean fieldpredictions agree with each other gives an indication oftheir reliability in this regime. Although the exact re-sults are terminated at the estimated limit of reliabilityof the energy, one should note that the results in the insetof Fig. 2 represent the difference between two relativelylarge terms, and the effects of any errors or approxima-tion are accordingly magnified.In Ref. 14, the non-monotonic behavior of the energydifference was attributed to two competing effects: on theone hand the thermal wavelength increases with decreas-ing temperature, and on the other the particles becomemore confined to their lattice positions as the tempera-ture decreases, which reduces the amount of overlappingwave function and non-zero symmetrization exchange. A v e r a g e E n e r gy , < H > / ħ w L J -505 H > + -< H > - ] / ħ w L J 246 0 1 2 3 4 5 A v e r a g e E n e r gy , Inverse Temperature, b ħ w LJ -10 0 1 2 3 [< H Inverse Temperature, b ħ w LJ FIG. 4: Average energy versus inverse temperature ( N = 4,∆ q = r e / κ = 0, λ = mω , d m = 0). The solid curveis the exact result using l max = 20 , 000 energy levels. Thetriangles are the singlet and the circles are the pair meanfield Monte Carlo simulations, respectively. The dotted curveis the classical result, h E i cl = N/β , and the dashed line isthe ground state energy. Inset. The boson energy minus thefermion energy using d m ≤ 2. The crosses are simulationswith the classical commutation function, W = 1. Figure 3 shows the density profile for this weak cou-pling case at β ¯ hω LJ = 2. The density peaks are ratherbroad, with the central two particles merging into a sin-gle peak. Interestingly enough, the density profile spillsover beyond the wall particles at q = 0 and q = 5 r e .There is good agreement between all four methods, withthe classical phase space simulations being closer to theexact results at the shoulders of the density profile.The inset of the figure shows the difference betweenthe density of bosons and that for fermions with sym-metrization effects accounted for by only nearest neigh-bor transpositions (dimers). The phase space simula-tions may again be described as quantitatively correct.Evidently the bulk of the symmetrization effects are cap-tured using the classical commutation function, W = 1.Figure 4 is for a high density case with lattice spacing∆ q = r e / 10. It can be seen that the exact energy ap-proaches the classical energy from above as the tempera-ture is increased. Without symmetrization effects (mainpart of figure) the mean field simulations lie between theexact and the classical results, with the pair mean fieldresults lying closer to the exact results than the singletmean field results across the temperature range shown.Both mean field results lie below the ground state en-ergy at low temperatures. For example, in this casethe ground state energy is E / ¯ hω LJ = 2 . β ¯ hω LJ = 10, the singlet mean field gives h ˆ Hi / ¯ hω LJ =2 . ± . . ± . . ± . d m = 2. It can be seen that there isa pole in the exact results for fermions at β ¯ hω LJ ≈ . d m = 4, are included.) It can be seen thatthe classical, W = 1, singlet mean field, and pair meanfield commutation functions all give this pole for d m ≤ β ¯ hω LJ ≤ . l max = 20 , 000 in this case. At intermediateand low temperatures, β ¯ hω LJ > ∼ 1, the classical results liecloser to the exact results than do the mean field results.It is difficult to obtain the results for fermions accuratelywhen the denominator passes through zero. In such aregime the neglected higher order terms in the expansioncontribute significantlyFigure 5 shows the density profile in this high densitycase at β ¯ hω LJ = 1. It can be seen that there is essentiallya single density peak in the center of the system, and thatthe density spills beyond the wall particles at q = 0 and q = 0 . 5. There is little to choose between the threesimulation algorithms, at least in the case of monomers(ie. no symmetrization effects, d m = 0). Compared tothe exact calculations all three simulation algorithms areslightly broader at the peak.For the case of bosons, Fig. 5B, including nearestneighbor transpositions makes the profile slightly nar-rower and more sharply peaked. There is again goodagreement between the three simulation algorithms, withthe pair mean field approach slightly underestimating theheight of the density peak. For the case of fermions,Fig. 5C, the exact calculations give a single peak, thatis narrower and higher than for bosons. There is nodip in the center of this peak, as one might have ex-pected from Fermi repulsion. (At the lower temperatureof β ¯ hω LJ = 2, the pair mean field profile shows a bifur-cated peak, but the other methods show a single peaksimilar to here). The classical and the singlet mean fieldapproaches underestimate the height of the peak, withthe latter being closer to the exact results than the for-mer. The pair mean field approach overestimates theheight of the peak. All three simulation algorithms missthe broad base to the density profile at the walls thatis given by the exact calculations. Arguably for reliableresults for fermions at this density and temperature, oneshould go beyond single dimer transpositions. IV. CONCLUSION This paper has been concerned with ascertaining theaccuracy of a quantum Monte Carlo algorithm for aninteracting system. The algorithm is based on a formu-lation of quantum statistical mechanics in classical phasespace and it uses a mean field approximation for the com-mutation function. For the tests reliable exact resultswere required as benchmarks, and these were obtained D e n s it y , r ( z ) r e (A)01234-0.25 0 0.25 0.5 0.75 D e n s it y , Position, z/r e D e n s it y , r ( z ) r e (B)01234-0.25 0 0.25 0.5 0.75 D e n s it y , Position, z/r e D e n s it y , r ( z ) r e (C)02468-0.25 0 0.25 0.5 0.75 D e n s it y , Position, z/r e FIG. 5: Density profile for β ¯ hω LJ = 1 ( N = 4, ∆ q = r e / κ = 0, λ = mω ). The solid curve is the exact result using l max = 20 , 000 energy levels. The crosses are the classical, W = 1, the triangles are the singlet, and the circles are thepair mean field Monte Carlo simulations, respectively. (A) d m = 0. (B) Bosons, d m ≤ (C) Fermions, d m ≤ for a one-dimensional harmonic crystal for which the en-ergy eigenvalues and eigenstates can be expressed analyt-ically in closed form. The analytic nature of the exact re-sults for this model allowed up to 20,000 energy levels tobe used to establish the benchmark results for the tests.Previous tests of the mean field classical phase space for-mulation of quantum statistical mechanics used bench-marks established for a one-dimensional Lennard-Jonesmodel with 50 energy eigenvalues obtained numerically. Here it was found that the number of energy levels has asignificant effect on the statistical average at higher tem-0peratures, and so the present benchmarks can be reliedupon in this regime.Two versions of the mean field approximation weretested: the singlet version, which has previously beenpublished, and a pair version, which is new here. Itwas found that the singlet mean field approach was quali-tatively correct over the whole temperature (and densityand coupling) regime studied. It appeared to be exactin the high temperature limit, and it was generally bet-ter than 10% accurate in the low temperature regimein which the system was predominately in the quantumground state. The pair mean field algorithm significantlyimproved the accuracy of the singlet algorithm in the in-termediate and low temperature regime. In the absenceof benchmark results, the difference between the singletand pair mean field predictions can be used as a guide tothe quantitative accuracy of the latter.There appear to be at least two advantages to thepresent mean field treatment of the commutation func-tion compared to evaluating it from high temperatureexpansions. First, the mean field expressions re-main accurate across the entire temperature regime, in-cluding the ground state. Second, algebraically the meanfield expressions are relatively simple, and computation-ally they are easy to implement and efficient to evalu-ate. In contrast the high temperature expansions rapidlybecome algebraically complex as higher order terms areincluded, and there are corresponding challenges intheir computational implementation and numerical eval-uation.The present paper also explored wave function sym- metrization effects. This was at the dimer level, whichmeans the transposition of nearest neighbor particles.It was found, somewhat surprisingly, that combiningclassical Monte Carlo in classical phase space with thesymmetrization function (ie. neglecting the commutationfunction) in some, but not all, cases gave as good resultsas those obtained retaining the mean field commutationfunction. At the highest density studied, some featuresof the symmetrized system were not captured entirely bythe present phase space simulations, particularly in thecase of fermions. This suggests that retaining furtherterms in the symmetrization loop expansion (eg. doubledimer and trimer) may be necessary in some cases. Itmay also be worth reflecting on the underlying philoso-phy of the mean field approach in the presence of wavefunction symmetrization.Finally, a rather interesting conclusion from thepresent and earlier results is that the classical com-ponent dominates, not just in the high temperature limitbut even in the quantum ground state (for structure,not energy). Of course quantum effects are not en-tirely absent, and when present these are captured bythe present mean field commutation function and alsothe symmetrization function, but it is clear that thesetruly are a perturbation on the classical prediction. Thisunderscores the utility of treating real world condensedmatter systems via quantum statistical mechanics for-mulated as an integral over classical phase space, ratherthan formulating it as a sum over quantum states, or byparameterizing the wave function. R. G. Parr and W. Yang, Density-Functional Theory ofAtoms and Molecules , (Oxford University Press, 2nd ed.1994). K. Morton and D. Mayers, Numerical Solution of PartialDifferential Equations, An Introduction , (Cambridge Uni-versity Press, 2nd ed. 2005). I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. , 885 (2008). J. M. McMahon, M. A. Morales, C. Pierleoni, and D. M.Ceperley, Rev. Mod. Phys. , 1607 (2012). L. Pollet, Rep. Prog. Phys. , 094501 (2012). A. Hernando and J. Van´ıˇcek, Phys. Rev. A , 062107(2013). arXiv:1304.8015v2 [quant-ph] (2013). P. Attard, Quantum Statistical Mechanics: Equilibriumand Non-Equilibrium Theory from First Principles , (IOPPublishing, Bristol, 2015). P. Attard, Entropy Beyond the Second Law. Thermody-namics and Statistical Mechanics for Equilibrium, Non-Equilibrium, Classical, and Quantum Systems , (IOP Pub-lishing, Bristol, 2018). P. Attard, “Quantum Statistical Mechanics in Classi-cal Phase Space. Expressions for the Multi-Particle Den-sity, the Average Energy, and the Virial Pressure”,arXiv:1811.00730 [quant-ph] (2018). P. Attard, “Quantum Statistical Mechanics in Classical Phase Space. Test Results for Quantum Harmonic Oscilla-tors”, arXiv:1811.02032 (2018). P. Attard, “Quantum Statistical Mechanics as an ExactClassical Expansion with Results for Lennard-Jones He-lium”, arXiv:1609.08178v3 [quant-ph] (2016). P. Attard, “Quantum Statistical Mechanics Results forArgon, Neon, and Helium Using Classical Monte Carlo”,arXiv:1702.00096 (2017). P. Attard, “Quantum Statistical Mechanics in ClassicalPhase Space. III. Mean Field Approximation Benchmarkedfor Interacting Lennard-Jones Particles”, arXiv:1812.03635[quant-ph] (2018). P. Attard, “Fermionic Phonons: Exact Analytic Resultsand Quantum Statistical Mechanics for a One DimensionalHarmonic Crystal”, arXiv:1903.06866 [quant-ph] (2019). A. Messiah, Quantum Mechanics , (North-Holland, Ams-terdam, Vols I and II, 1961). E. Wigner, Phys. Rev. , 749, (1932). J. G. Kirkwood, Phys. Rev. , 31, (1933). P. Attard, Thermodynamics and Statistical Mechanics:Equilibrium by Entropy Maximisation (Academic Press,London, 2002). S. W. van Sciver,