Quantum Monte Carlo in Classical Phase Space. Mean Field and Exact Results for a One Dimensional Harmonic Crystal

TL;DR

This paper applies Monte Carlo simulations in classical phase space to a one-dimensional quantum harmonic crystal, analyzing symmetrization effects and comparing mean field approximations with exact results across temperature regimes.

Contribution

It introduces a phase space Monte Carlo method for quantum crystals and evaluates mean field approximations against exact solutions, including a novel cluster mean field approach.

Findings

High-temperature singlet mean field is very accurate.

Pair mean field improves accuracy at lower temperatures.

Cluster mean field accounts for non-commutativity and is extendable to 3D.

Abstract

Monte Carlo simulations are performed in classical phase space for a one-dimensional quantum harmonic crystal. Symmetrization effects for spinless bosons and fermions are quantified. The algorithm 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, and that the pair mean field approximation gives a systematic improvement in the intermediate and low temperature regime. The latter is derived from a cluster mean field approximation that accounts for the non-commutativity of position and momentum, and that can be applied in three dimensions.

Equations

$${⟨\hat{\mathcal{H}}⟩}_{N,V,T}^{\pm }$$

$$Z^{\pm }(N,V,T)$$

$$|\mathbf{q}⟩=\delta (\mathbf{r}-\mathbf{q}),\text{mbox}and|\mathbf{p}⟩=\frac{e^{-\mathbf{p}\cdot \mathbf{r}/i\hslash }}{V^{N/2}}.$$

$${\eta }_q^{\pm }(\mathbf{p},\mathbf{q})\equiv \frac{1}{⟨\mathbf{p}|\mathbf{q}⟩}\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p⟨\hat{\mathrm{P}}\mathbf{p}|\mathbf{q}⟩.$$

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

$${\eta }_{q;1\dots l}^{\pm (l)}=(\pm 1{)}^{l-1}{e}^{{\mathbf{q}}_{1l}\cdot {\mathbf{p}}_l/i\hslash }\prod _{j=1}^{l-1}{e}^{{\mathbf{q}}_{j+1,j}\cdot {\mathbf{p}}_j/i\hslash },$$

$$e^{-\beta \mathcal{H}(\mathbf{p},\mathbf{q})}{W}_p(\mathbf{p},\mathbf{q})=\frac{⟨\mathbf{q}|e^{-\beta \hat{\mathcal{H}}}|\mathbf{p}⟩}{⟨\mathbf{q}|\mathbf{p}⟩}.$$

$$e^{-\beta \mathcal{H}(\mathbf{p},\mathbf{q})}{W}_p(\mathbf{p},\mathbf{q})$$

$$U(\mathbf{q})$$

$$U_j({\mathbf{q}}_j;\mathbf{q})$$

$$U_j({\mathbf{q}}_j;\mathbf{q})={\overline{U}}_j(\mathbf{q})+\frac{1}{2}[{\mathbf{q}}_j-{\overline{\mathbf{q}}}_j][{\mathbf{q}}_j-{\overline{\mathbf{q}}}_j]:{\overline{\underline{\underline{U}}}}_j^{\prime \prime },$$

$${\omega }_{j\alpha }(\mathbf{q})=\sqrt{{\lambda }_{j\alpha }(\mathbf{q})/m},\alpha =x,y,\dots ,d.$$

$$U(\mathbf{q})$$

$${\mathcal{H}}^{\mathrm{SHO}}(\mathbf{p},\mathbf{q}-\overline{\mathbf{q}})=\sum _{j=1}^N{\overline{U}}_j+\frac{1}{2}\sum _{j,\alpha }\hslash {\omega }_{j\alpha }\left[P_{j\alpha }^2+Q_{j\alpha }^2\right].$$

$$W_p^{\mathrm{mf}}(\mathbf{\Gamma })$$

$$W_{p,j\alpha }^{\mathrm{SHO}}(P_{j\alpha },Q_{j\alpha })$$

$$U(\mathbf{q})$$

$$U_{\alpha }({\mathbf{q}}_{\alpha };\mathbf{q})$$

$$U_{\alpha }({\mathbf{q}}_{\alpha };\mathbf{q})={\overline{U}}_{\alpha }(\mathbf{q})+\frac{1}{2}{\overline{\underline{\underline{U}}}}_{\alpha }^{\prime \prime }:[{\mathbf{q}}_{\alpha }-{\overline{\mathbf{q}}}_{\alpha }][{\mathbf{q}}_{\alpha }-{\overline{\mathbf{q}}}_{\alpha }].$$

$$W^{\mathrm{mf}}(\mathbf{\Gamma })=\prod _{\alpha }{W}_{\alpha }^{\mathrm{SHO}}({\mathbf{P}}_{\alpha },{\mathbf{Q}}_{\alpha })=\prod _{\alpha ,b}{W}_{\alpha ,b}^{\mathrm{SHO}}(P_{\alpha ,b},Q_{\alpha ,b}).$$

$$U(\mathbf{q})$$

$$U_j(q_j;\mathbf{q})$$

$${\overline{d}}_j(\mathbf{q})$$

$${\overline{U}}_j^{\prime \prime }\equiv {\nabla }_j{\nabla }_j{U}_j(q_j;\mathbf{q})=\kappa +\lambda +\frac{\lambda }{2}[{\delta }_{j1}+{\delta }_{jN}],$$

$$U_j(q_j;\mathbf{q})$$

$$U(\mathbf{q})$$

$$U_{\alpha }({\mathbf{q}}_{\alpha };\mathbf{q})$$

$${\overline{\underline{\underline{U}}}}_{\alpha }^{\prime \prime }$$

$$\equiv$$

$$\left(\begin{array}[]{c}\overline{d}_{2\alpha-1}\\ \overline{d}_{2\alpha}\end{array}\right)=\frac{-1/2}{K_{\alpha}^{\prime}K_{\alpha}^{\prime\prime}-1}\left(\begin{array}[]{c}K_{\alpha}^{\prime\prime}d_{2\alpha-2}-d_{2\alpha+1}\\ -d_{2\alpha-2}+K_{\alpha}^{\prime}d_{2\alpha+1}\end{array}\right).$$

$${\mu }_{\alpha }^{\pm }$$

Full text

I Introduction

The impediment to applying quantum statistical mechanics to condensed matter systems at terrestrial temperatures and densities is the horrendous scaling with system size entailed by conventional approaches. Parr94 ; Morton05 ; Bloch08 ; McMahon12 ; Pollet12 ; Hernando13 Although the field is not exactly moribund, it would be fair to say that the rate of progress has been disappointing, and that there appears no obvious way of overcoming this fundamental limitation of the existing methods.

The pressing need to approach the problem from a different direction has motivated the author to develop a new methodology that is based on a formally exact transformation that expresses the quantum partition function as an integral over classical phase space.QSM ; STD2 ; Attard18a The validity of the formulation has been verified analytically and numerically for the quantum ideal gas,QSM ; STD2 and for non-interacting quantum harmonic oscillators.Attard18b

Application of the algorithm to interacting systems indicate that, for a given statistical error, the algorithm scales sub-linearly with system size.Attard16 ; Attard17 ; Attard18c This places it in a class of its own for the treatment of condensed matter quantum systems.

Validation of the algorithm for an interacting system has been provided by quantitative tests for interacting Lennard-Jones particles against benchmarks obtained with more conventional methods by Hernando and Vaníček.Hernando13 The latter results were exact, apart from the fact that they were numerical and ‘only’ 50 energy levels were used for the statistical averages. One should not understate the difficulty of obtaining numerically 50 energy levels for the Lennard-Jones system; it is undoubtedly a greater computational challenge than obtaining analytically the 20,000 energy levels used below for the present one-dimensional harmonic crystal. If nothing else the figure of 50 levels does underscore the insurmountable intractability of conventional quantum approaches.

In fact the original motivation for the present paper was to sort out a small discrepancy between the putative exact results of Hernando and VaníčekHernando13 and the author’s mean field, classical phase space results at the highest temperature studied.Attard18c It was unclear whether the differences 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 the exact results. Accordingly the author has undertaken to establish his own exact results for use as benchmarks, and it was found that 20,000 energy levels were necessary for reliable results at temperatures high enough to give a departure from the ground state.

The results of these tests are reported below. It is found that the mean field approximation as originally formulated is essentially exact at the highest temperatures, but the error can be on the order of 5–10% at intermediate and low temperatures, depending upon model parameters. Accordingly, this papers develops a systematic general improvement to the phase space algorithm that might be called the cluster mean field approximation. This is here implemented at the singlet level, which is the original version, and also at the pair level, which is new. It is found that the pair mean field approximation gives essentially exact results at intermediate temperatures. At low temperatures such that the system is predominantly in the ground state, the singlet and pair mean field approximations are found to be in error by on the order of 5%.

This paper relies upon earlier work, which will not be re-derived here. The reader is referred Ref. STD2, for the derivation of quantum statistical mechanics, to Ref. Attard18a, for the derivation of the phase space formulation (see also Ref. Attard19, for some formal details concerning symmetrization of multiparticle states), and to Ref. Attard19, for the exact phonon analysis of the present one-dimensional harmonic crystal.

II Analysis and Model

II.1 Phase Space Formulation

The details of the classical phase space formulation of quantum statistical mechanics can vary with the particular quantity being averaged.Attard18a For the average energy, which represents a certain class of operators, namely those that are a linear combination of functions purely of the momentum operator or of the position operator, the canonical equilibrium average can be writtenAttard18a

[TABLE]

where the partition function is

[TABLE]

In these $N$ is the number of particles, assumed identical and spinless, $V$ is the volume, $T$ is the temperature, $\beta=1/k_{\mathrm{B}}T$ is the inverse temperature, $h$ is Planck’s constant, and $k_{\mathrm{B}}$ is Boltzmann’s constant. Also, ${\bf\Gamma}=\{{\bf p},{\bf q}\}$ is a point in phase space, with the vector of particles’ momenta being ${\bf p}=\{{\bf p}_{1},{\bf p}_{2},\ldots,{\bf p}_{N}\}$, and that of the particles’ position being ${\bf q}=\{{\bf q}_{1},{\bf q}_{2},\ldots,{\bf q}_{N}\}$. Finally, $\hat{\cal H}={\cal H}(\hat{\bf p},\hat{\bf q})$ is the energy or Hamiltonian operator, and ${\cal H}({\bf p},{\bf q})$ is the classical Hamiltonian function.

The plus sign is for bosons and the minus sign is for fermions. Actually, since this has been formulated for spinless particles, these refer to the fully symmetrized and fully anti-symmetrized spatial part of the wave function (see Appendix C of Ref. Attard19, ). Throughout the words ‘boson’ and ‘fermion’ should be understood in this sense.

The unsymmetrized position and momentum eigenfunctions in the position representation are respectivelyMessiah61

[TABLE]

The symmetrization function is defined asAttard18a

[TABLE]

Here the sum is over the $N!$ permutation operators $\hat{\mathrm{P}}$, whose parity is $p$. The imaginary part of these is odd in momentum, $\eta^{\pm}_{q}({\bf p},{\bf q})^{*}=\eta^{\pm}_{q}(-{\bf p},{\bf q})$.

The symmetrization function 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. The prime signifies that the sum is over unique loops (ie. each configuration of particles in loops occurs once only) with each index different (ie. no particle may belong to more than one loop). In general the $l$-loop symmetrization factor is

[TABLE]

where ${\bf q}_{jk}\equiv{\bf q}_{j}-{\bf q}_{k}$. This corrects a typographical error in Eq. (3.4) of Ref. Attard18a,

The commutation function, which is essentially the same as the function introduced by WignerWigner32 and analyzed by Kirkwood,Kirkwood33 is defined by

[TABLE]

Again one has $W_{p}({\bf p},{\bf q})^{*}=W_{p}(-{\bf p},{\bf q})$. High temperature expansions for the commutation function have been given.Wigner32 ; Kirkwood33 ; STD2 ; Attard18b

The commutation function in phase space can also be written as a series of energy eigenfunctions and eigenvalues, $\hat{\cal H}|{\bf n}\rangle=E_{n}|{\bf n}\rangle$. Using the completeness properties of these one obtains

[TABLE]

This exact expression forms the basis of the mean field approximation to the commutation function.

II.2 Cluster Mean Field Approximation

In the classical phase space formulation of quantum statistical mechanics, the symmetrization function is relatively trivial to obtain and implement. The commutation function is more of a challenge, with the most successful approach using a mean field approximation that exploits the analytic form of the commutation function in the case of independent simple harmonic oscillators.Attard18b This has previously been tested for the simulation of a Lennard-Jones system. Attard18c

This section begins with a summary of the singlet mean field approximation to the commutation function.Attard18b Then the cluster mean field approximation is given.

II.2.1 Singlet Mean Field Approximation

In general, the particles of the sub-system interact via the potential energy, which is the sum of one-body, two-body, three-body terms, etc.,

[TABLE]

Distributing the energy equally, the energy of particle $j$ can be defined as

[TABLE]

with $U({\bf q})=\sum_{j=1}^{N}U_{j}({\bf q}_{j};{\bf q})$. The argument $({\bf q}_{j};{\bf q})$ means that ${\bf q}_{j}$ is here separated out from ${\bf q}$.

The potential energy of particle $j$ in configuration ${\bf q}$ may be expanded to second order about its local minimum at $\overline{\bf q}_{j}({\bf q})$,

[TABLE]

where the minimum value of the potential is $\overline{U}_{j}({\bf q})\equiv U_{j}(\overline{\bf q}_{j};{\bf q})$. The $d\times d$ second derivative matrix for particle $j$ at the minimum, $\overline{\underline{\underline{U}}}_{j}^{\prime\prime}=\left.\nabla_{j}\nabla_{j}U_{j}({\bf q}_{j};{\bf q})\right|_{{\bf q}_{j}=\overline{\bf q}_{j}}$, is assumed positive definite.

For configurations ${\bf q}$ that have no local minimum in the potential, or that have too large a displacement $|{\bf q}_{j}-\overline{\bf q}_{j}|$, the corresponding single particle commutation function can be set to unity, $W_{j}({\bf\Gamma})=1$, (or, in the multi-dimensional case, the commutation function of the corresponding mode). This is justified by analytic results for the simple harmonic oscillator.Attard18b

The positive definite second derivative matrix has $d$ eigenvalues $\lambda_{j\alpha}({\bf q})>0$, and orthonormal eigenvectors, $\overline{\underline{\underline{U}}}_{j}^{\prime\prime}{\bf X}_{j\alpha}=\lambda_{j\alpha}{\bf X}_{j\alpha}$, $\alpha=x,y,\ldots,d$. For molecule $j$ in configuration ${\bf q}$ the eigenvalues define the frequencies

[TABLE]

With this the potential energy is

[TABLE]

Here $Q_{j\alpha}\equiv\sqrt{{m\omega_{j\alpha}}/{\hbar}}\,Q_{j\alpha}^{\prime}$, and ${\bf Q}_{j}^{\prime}=\underline{\underline{X}}_{j}^{\mathrm{T}}[{\bf q}_{j}-\overline{\bf q}_{j}]$. (This corrects a typographical error in Eqs (2.7) and (2.8) in Ref. Attard18c, .)

The mean field approximation combined with the second order expansion about the local minima maps each configuration ${\bf\Gamma}$ to a system of $dN$ independent harmonic oscillators with frequencies $\omega_{j\alpha}$ displacements $Q_{j\alpha}$, and momenta $P_{j\alpha}=\left\{\underline{\underline{X}}_{j}^{\mathrm{T}}{\bf p}_{j}\right\}_{\alpha}/\sqrt{m\hbar\omega_{j\alpha}}$. (This corrects a typographical error in Ref. Attard19, .)

With this harmonic approximation for the potential energy, the effective Hamiltonian in a particular configuration can be written

[TABLE]

The commutation function for the interacting system for a particular configuration can be approximated as the product of commutation functions for effective non-interacting harmonic oscillators which have the local displacement as their argument. With this the mean field commutation function is

[TABLE]

The harmonic oscillator commutation function for a single mode isAttard18b

[TABLE]

The prefactor $e^{-iP_{j\alpha}Q_{j\alpha}}$ corrects the prefactor $e^{-ip_{j\alpha}q_{j\alpha}/\hbar}$ given in Eq. (5.10) of Ref. [Attard18b, ]. Here $\mathrm{H}_{n}(z)$ is the Hermite polynomial of degree $n$. The imaginary terms here are odd in momentum. As justified by analytic results for the simple harmonic oscillator,Attard18b for configurations such that $\overline{\underline{\underline{U}}}_{j}^{\prime\prime}({\bf q})$ is not positive definite (ie. a particular eigenvalue is not positive, $\lambda_{j\alpha}\leq 0$), or that the displacement $Q_{j\alpha}$ exceeds a predetermined cut-off, the corresponding commutation function can be set to unity, $W_{p,{j\alpha}}^{\mathrm{SHO}}=1$.

For the averages, the momentum integrals can be performed analytically, both here and in combination with the symmetrization function. This considerably reduces computer time and substantially increases accuracy.

II.2.2 Cluster Mean Field Approximation

Any configuration ${\bf q}$ can be decomposed into disjoint clusters labeled $\alpha=1,2,\ldots$. Of the different criteria that can be used to define a cluster, perhaps the simplest is that two particles belong to the same cluster if, and only if, they are connected by at least one path of bonds. Two particles are bonded if their separation is less than a nominated length. Some clusters, perhaps the great majority, will consist of a single particle.

An even simpler definition can be made in one dimension. In this case define the pair cluster $\alpha$, $\alpha=1,2,\ldots,N/2$, as the nearest neighbors $\{2\alpha-1,2\alpha\}$, irrespective of their actual separation. This criteria is used in the results presented below.

Using a separation-based criterion for the definition of a cluster is useful not only for the mean field approximation to the commutation function, but also for the calculation of the symmetrization function. Depending on the chosen bond length, only permutations of particles within the same cluster need to be considered. (This idea is not used in the results presented below.)

The cluster energy is the internal energy plus the relevant proportion of the interaction energy with other clusters: half for pair interactions, one third for triplet interactions, etc. The total potential energy is

[TABLE]

where the energy of cluster $\alpha$ is

[TABLE]

There are $N_{\alpha}$ particles in cluster $\alpha$, with positions ${\bf q}_{\alpha}=\{{\bf q}_{\alpha,1},{\bf q}_{\alpha,2},\ldots,{\bf q}_{\alpha,N_{\alpha}}\}$, where ${\bf q}_{\alpha,j}={\bf q}_{k}$ for one of the $k\in\alpha$. This is a $(dN_{\alpha})$-dimensional vector.

The second order expansion about the minimum energy cluster configuration, $\overline{\bf q}_{\alpha}$, is

[TABLE]

The second derivative matrix is $\overline{\underline{\underline{U}}}_{\alpha}^{\prime\prime}\equiv\left.\nabla_{\alpha}\nabla_{\alpha}U_{\alpha}({\bf q}_{\alpha};{\bf q})\right|_{\overline{\bf q}_{\alpha}}$, which is $(dN_{\alpha})\times(dN_{\alpha})$. One has to find the eigenvalues, assumed positive, and eigenvectors of this. These give the cluster phonon mode frequencies $\omega_{\alpha,b}$, $b=1,2,\ldots,dN_{\alpha}$, and mode amplitudes, ${\bf Q}_{\alpha}=\underline{\underline{\omega}}_{\alpha}\underline{\underline{X}}_{\alpha}^{\mathrm{T}}[{\bf q}_{\alpha}-\overline{\bf q}_{\alpha}]$. The frequency matrix is diagonal with elements $\omega_{\alpha;bc}=\sqrt{m\omega_{\alpha,b}/\hbar}\,\delta_{bc}$. As before, the momentum of mode $b$ in cluster $\alpha$ is $P_{\alpha,b}=\left\{\underline{\underline{X}}_{\alpha}^{\mathrm{T}}{\bf p}_{\alpha}\right\}_{b}/\sqrt{m\hbar\omega_{\alpha,b}}$.

This formulation is essentially the same as the singlet mean field theory, and one may similarly define the cluster mean field commutation function as the product of simple harmonic oscillator commutation functions, one for each phonon mode of each cluster. The cluster mean field commutation function of the configuration is

[TABLE]

From the computational point of view, a felicitous aspect of the cluster mean field approximation is that there are exactly as many modes as in the singlet mean field approximation. This means that all of the sub-routines called to obtain the various statistical averages in the singlet approximation can be called without change in the cluster mean field approximation.

II.3 Harmonic Crystal Model

II.3.1 Potential Energy

Following earlier work,TDSM ; Attard19 consider a one-dimensional harmonic crystal in which the particles are attached by linear springs to each other and to lattice sites. Let the coordinate of the $j$th particle be $q_{j}$, and let its lattice position (ie. in the lowest energy state) be $\tilde{\overline{q}}_{j}=j\Delta_{q}$. The tilde signifies that these are ordered. The lattice spacing is also the relaxed inter-particle spring length. There are fixed ‘wall’ particles at $q_{0}=0$ and $q_{N+1}=(N+1)\Delta_{q}$. Let $d_{j}\equiv q_{j}-\tilde{\overline{q}}_{j}$ be the displacement from the lattice position; for the wall particles, $d_{0}=d_{N+1}=0$. The system has over-all number density $\rho=\Delta_{q}^{-1}$.

There is an external harmonic potential of spring constant $\kappa$ acting on each particle centered at its lattice site. The inter-particle spring has strength $\lambda$ and relaxed length $\Delta_{q}$. With these the potential energy is

[TABLE]

The energy eigenfunctions and eigenvectors can be obtained explicitly for this model by expressing it in terms of phonon modes.Attard19

This model potential is not invariant with respect to the permutation of the positions of the particles, and it therefore violates a fundamental axiom of quantum mechanics. Messiah61 This was discussed in detail in Appendix A of Ref. Attard19, , where the origin, interpretation, and justification for the potential was given. To those remarks may be added the fact that the formulation of quantum statistical mechanics in classical phase space is unchanged by a non-symmetric potential (unpublished).

II.3.2 Singlet Mean Field

The energy of particle $j$ in configuration ${\bf q}$ is

[TABLE]

The total potential energy is just $U({\bf q})=\sum_{j=1}^{N}U_{j}({\bf q}_{j};{\bf q})$, and $d_{j}\equiv q_{j}-\tilde{\overline{q}}_{j}$.

The gradient vanishes when

[TABLE]

The second derivative is

[TABLE]

which is independent of the configuration ${\bf q}$.

The potential energy of particle $j$ in configuration ${\bf q}$ may be expanded to second order about its local minimum at $\overline{q}_{j}({\bf q})$,

[TABLE]

where $\overline{U}_{j}({\bf q})\equiv U_{j}(\overline{q}_{j}({\bf q});{\bf q})$. This second order expansion for the potential is exact for the present harmonic crystal. Note that the most likely position of particle $j$ for the current configuration, $\overline{q}_{j}({\bf q})$, is not the same as its lattice position $\tilde{\overline{q}}_{j}$.

For each molecule define the frequency $\omega_{j}=\sqrt{\overline{U}^{\prime\prime}_{j}/m}$. This is the same for all configurations ${\bf q}$. With this the potential energy is

[TABLE]

Here $Q_{j}\equiv\sqrt{{m\omega_{j}}/{\hbar}}\,[{q}_{j}-\overline{q}_{j}]$. As above, the momenta are ${P}_{j}={p}_{j}/\sqrt{m\hbar\omega_{j}}$.

This is the approximation in the singlet mean field approach: each frequency mode is a displacement of a single particle. One may now apply Eqs (2.14), (2.15), and (II.2.1) for the singlet mean field commutation function.

II.3.3 Pair Mean Field Approximation

For the one dimensional harmonic crystal, define a cluster pair $\alpha$ as the nearest neighbors $\{2\alpha-1,2\alpha\}$, $\alpha=1,2,\ldots,\lfloor N/2\rfloor$. For simplicity, assume $N$ even. (Contrariwise, include particle $N$ as a singlet cluster.)

The energy of cluster $\alpha$ is

[TABLE]

where the displacement is $d_{j}=q_{j}-\tilde{\overline{q}}_{j}$. Note that the interaction with the wall particles, when present, has to be counted fully. Note also that $d_{0}=d_{N+1}=0$.

The second derivative matrix is

[TABLE]

Using it gives the optimum cluster displacement as

[TABLE]

The eigenvalues of the second derivative matrix (without $-\lambda$) are

[TABLE]

Since the $K_{\alpha}$ are negative, so are the $\mu_{\alpha}^{\pm}$.

Writing the eigenvectors as ${\bf u}_{\alpha}^{\pm}=c_{\alpha}^{\pm}(1,a_{\alpha}^{\pm})$, $c_{\alpha}^{\pm}\equiv 1/\sqrt{1+a_{\alpha}^{\pm 2}}$, from the eigenvalue equation one obtains $a_{\alpha}^{\pm}=\mu_{\alpha}^{\pm}-K_{\alpha}^{\prime}$. This is presumably equivalent to $a_{\alpha}^{\mp}=1/[\mu_{\alpha}^{\mp}-K_{\alpha}^{\prime\prime}]$.

The mode frequencies are $\omega_{\alpha}^{\pm}=\sqrt{-\lambda\mu_{\alpha}^{\pm}/m}$, which are real. The mode amplitude is

[TABLE]

where $\Delta d_{2\alpha-1}\equiv d_{2\alpha-1}-\overline{d}_{2\alpha-1}=q_{2\alpha-1}-\overline{q}_{2\alpha-1}$, and $\Delta d_{2\alpha}\equiv d_{2\alpha}-\overline{d}_{2\alpha}=q_{2\alpha}-\overline{q}_{2\alpha}$.

The contribution to the total energy from cluster $\alpha$ in configuration ${\bf\Gamma}$ is

[TABLE]

In the computational implementation of the algorithm, these can be used to replace directly the singlet mean field terms, $P_{\alpha}^{+}\Rightarrow P_{2\alpha-1}$, $P_{\alpha}^{-}\Rightarrow P_{2\alpha}$, and similarly for $Q_{\alpha}^{\pm}$ and $\omega_{\alpha}^{\pm}$.

II.3.4 Symmetrization Function

Symmetrization consists of a sum over all particle permutations. However because of the highly oscillatory Fourier contributions to the loop symmetrization function, Eq. (2.6), only permutations of closely separated particles actually contribute to the statistical average.

In view of this one can define the permutation length, Attard19

[TABLE]

One can see that $d_{\mathrm{m}}=0$ corresponds to the identity permutation, $d_{\mathrm{m}}=2$ corresponds to a single nearest neighbor transposition (dimer), and $d_{\mathrm{m}}=4$ corresponds to either two distinct nearest neighbor transpositions (double dimer), or else a single cyclic permutation of three consecutive particles (trimer). One expects that the contributions to the symmetrized wave function will decrease with increasing permutation length.

Hence one can set an upper limit on the length of the permutations that are included in the symmetrization function. The numerical results below show that by far the greatest contribution comes from the identity permutation alone, $d_{\mathrm{m}}=0$. In some cases a measurable change occurs by including also nearest neighbor permutations, $d_{\mathrm{m}}=2$. Measurable but smaller change occurs upon also including permutations of length $d_{\mathrm{m}}=4$ (not shown below). For $N$ particles, the number of permutation terms that contribute to the symmetrization function is 1 for $d_{\mathrm{m}}=0$, $1+(N-1)$ for $d_{\mathrm{m}}\leq 2$, and $N+(N-2)(N-3)/2+2(N-2)$ for $d_{\mathrm{m}}\leq 4$.

II.4 Simulation Algorithm

The simulation algorithm was as previously described.Attard18c Briefly the Metropolis algorithm in position space was used with the usual classical Maxwell-Boltzmann weight. The various momentum integrals were performed analytically. Averages were evaluated by umbrella sampling using the commutation function and symmetrization function as weight. Since three versions of the commutation function (unity (ie. classical), singlet mean field, pair mean field), and three versions of the symmetrization function ($d_{\mathrm{m}}=0$ (ie. classical), and $d_{\mathrm{m}}\leq 2$, for bosons and for fermions) were tested, 9 different averages for each quantity were obtained simultaneously.

Typically enough configurations were generated to make the statistical error less than 1%, sometimes much less. In the simulations the time depends on how many Hermite polynomials are used for the commutation function (eight in the results reported below; tests with six and twelve showed no great effect), and the cut-off for the mode amplitude beyond which the commutation function was set to unity ($Q^{\mathrm{cut}}=1$ in the results below; tests with $Q^{\mathrm{cut}}=2$ showed no great effect).

Interestingly enough, for an accuracy of about 1%, the Monte Carlo algorithm was a factor of about $2,000$ times more efficient (in terms of total computer time) than the quasi-analytic exact phonon method.Attard19 The main bottleneck in the latter was the crude numerical quadrature method that was used to evaluate the symmetrization function and density profile, and this was exacerbated by the large number of energy levels that were required for accurate results at higher temperatures. For this particular comparison at $\beta\hbar\omega_{\mathrm{LJ}}=2$, $l^{\mathrm{max}}=20,000$ energy levels were necessary; grid parameters $\Delta_{Q}=0.15$ and $Q^{\mathrm{max}}=6$ gave a quadrature error of about .8%. The phonon method is more efficient in one respect, namely that it requires negligible computer time for each additional temperature point; the simulations give results for only one temperature at a time.

The Lennard-Jones frequency used to scale the results below is $\omega_{\mathrm{LJ}}=3.28\times 10^{12}\,$Hz, the mass is $m=3.35\times 10^{-26}\,$kg, the well-depth is $\varepsilon=4.93\times 10^{-22}\,$J, and the equilibrium separation is $r_{\mathrm{e}}=3.13\times 10^{-10}\,$m. These parameters are appropriate for neon.Sciver12

III Results

Figure 1 shows the average energy as a function of inverse temperature using exact calculations and classical phase space simulations. At low temperatures the phonon modes are in the ground state, $E_{1}=\sum_{n=1}^{N}\hbar\omega_{n}/2$, and at high temperatures the system approaches the classical result, $\langle{\cal H}\rangle_{\mathrm{cl}}=N/\beta$. It can be seen in the main part of the figure that both the singlet and pair mean field approximations are in relatively good agreement with these limiting results and with the exact calculations over the whole temperature regime shown. The commutation function provides a significant correction to the classical results at high temperatures. In the regime of Fig. 1, symmetrization effects are negligible and only the $d_{\mathrm{m}}=0$ calculations are shown.

There are two approximations in the exact calculations: the number of energy levels used and the domain and spacing of the grid used for the numerical quadrature. (I persist in calling these results ‘exact’ because they use explicit analytic expressions for the energy eigenvalues and eigenfunctions.) The quadrature affects only the average density profile without symmetrization, and also the average energy and the average density profile with symmetrization. Hence the exact calculations of the average energy in the absence of symmetrization effects, $d_{\mathrm{m}}=0$, as in Fig. 1, are approximate only as regards to the number of energy levels that are retained.

The exact calculations in Fig. 1 use 10,000 energy levels. These are adequate for low and intermediate temperatures, $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 0.24$, judged in part by comparison with results using 5,000 energy levels. The exact data begins to underestimate the classical result for higher temperatures than this, and are not shown in Fig. 1. One might speculate that the exact quantum result for the average energy of the harmonic crystal should approach the classical limit from above. Using fewer energy levels reduces the domain of inverse temperatures in which the exact calculations are reliable.

Both the singlet and pair mean field simulations are practically exact at higher temperatures, $\beta\hbar\omega_{\mathrm{LJ}}\lesssim 0.5$ in this case. As the temperature is decreased, the singlet mean field energy lies between the exact energy and the classical energy. The pair mean field result lies between the singlet mean field energy and the exact energy. It can be seen that at low temperatures the classical energy is substantially less than the exact ground state energy, but the pair mean field energy is only slightly less than the exact ground state energy. It may be concluded that the mean field approach is better than a high temperature expansion in that it yields the dominant quantum correction to the classical result over the entire range of temperatures.

The inset of Fig. 1 scales the average energy by the inverse temperature and focusses on the low temperature regime. For inverse temperatures $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 3$, the exact results with 10,000 energy levels are practically indistinguishable from the ground state energy. At low temperatures, both mean field classical phase space approximations give a lower energy than the ground state. For example, in the case of Fig. 1, the ground state energy is $E_{1}/\hbar\omega_{\mathrm{LJ}}=3.385$. At $\beta\hbar\omega_{\mathrm{LJ}}=10$, the singlet mean field theory gives $\langle\hat{\cal H}\rangle/\hbar\omega_{\mathrm{LJ}}=3.116\pm.002$, and the pair mean field theory gives $3.305\pm.001$. At this temperature the classical result given by the simulation was $0.4002\pm.0002$, which is rather close to the exact classical result of 0.4. Here and throughout, the statistical error quoted for the simulations is twice the standard error on the mean, which is the 96% confidence level.

One can conclude from the data that at low temperatures such that the system is close to the ground state, the mean field approximations remain viable. The pair mean field approach substantially reduces the error in the singlet mean field approach. In the absence of exact data, the difference between the pair and singlet mean field results would give a guide to the quantitative accuracy of the former.

It is worth mentioning that at each temperature the classical phase space simulations took about five minutes on a desktop personal computer to obtain the quoted accuracy. In comparison, it took about 2 days to obtain the exact results with these energy levels and quadrature grid, the latter being the time limiting part of the computations. (This is independent of how many temperature points are saved.)

Figure 2 shows the average energy for a weakly coupled crystal. There is no singlet potential, $\kappa=0$, and the nearest neighbor spring constant has been decreased by a factor of 50, $\lambda/m\omega_{\mathrm{LJ}}^{2}=0.02$. The lattice spacing and spring length is unchanged. The ground state energy is $E_{1}=0.3757$, which is approximately one tenth that for the parameters of Fig. 1.

In the main body of Fig. 2 the results of both mean field approximations appear indistinguishable from the exact results. At $\beta\hbar\omega_{\mathrm{LJ}}=3$, the exact energy is 1.3707 (for $l_{\mathrm{max}}=20,000$; 1.3627 for 10,000), the singlet mean field gives $1.3322\pm.0009$ and the pair mean field gives $1.3443\pm.0008$. The exact classical result here is 4/3, compared to $1.3329\pm.0007$ given by the simulations. These results are for $d=0$, so no symmetrization effects are included. For inverse temperatures $\beta\hbar\omega_{\mathrm{LJ}}\lesssim 1.8$, the mean field algorithms appear more reliable than the exact results with $l_{\mathrm{max}}=20,000$.

The inset of Fig. 2 shows the average energy for bosons less that for fermions, $[\langle\hat{\cal H}\rangle^{+}-\langle\hat{\cal H}\rangle^{-}]/\hbar\omega_{\mathrm{LJ}}$, $d_{\mathrm{m}}\leq 2$. Recall that for these spinless particles, ‘bosons’ means the fully symmetrized spatial energy eigenfunctions, and ‘fermions’ means the fully anti-symmetrized spatial energy eigenfunctions. At lower temperatures, $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 1.5$, the difference is positive, which means that the energy for bosons is greater than that for fermions. It can be seen that the peak difference, which occurs at $\beta\hbar\omega_{\mathrm{LJ}}=1.9$, is about 3% of the actual energy (exact results). (The error in the numerical quadrature used for the exact result is 2% at this temperature for the energy with $d_{\mathrm{m}}=0$. Hopefully for $d_{\mathrm{m}}\leq 2$ this error is the same for bosons as for fermions and therefore cancels.) The classical phase space results may be described as qualitatively correct and perhaps semi-quantitative in accuracy. For $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 2$, the singlet mean field approximation overestimates the energy difference, whereas the pair mean field approximation perhaps halves the error. The classical results, with commutation function $W=1$, performs surprisingly well in this low temperature regime. At higher temperatures than this all four approaches indicate that the energy difference turns negative. The extent to which the singlet and pair mean field predictions agree with each other gives an indication of their reliability in this regime. Although the exact results are terminated at the estimated limit of reliability of the energy, one should note that the results in the inset of Fig. 2 represent the difference between two relatively large terms, and the effects of any errors or approximation are accordingly magnified.

In Ref. Attard19, , the non-monotonic behavior of the energy difference was attributed to two competing effects: on the one hand the thermal wavelength increases with decreasing temperature, and on the other the particles become more confined to their lattice positions as the temperature decreases, which reduces the amount of overlapping wave function and non-zero symmetrization exchange.

Figure 3 shows the density profile for this weak coupling case at $\beta\hbar\omega_{\mathrm{LJ}}=2$. The density peaks are rather broad, with the central two particles merging into a single peak. Interestingly enough, the density profile spills over beyond the wall particles at $q_{0}=0$ and $q_{5}=5r_{\mathrm{e}}$. There is good agreement between all four methods, with the classical phase space simulations being closer to the exact results at the shoulders of the density profile.

The inset of the figure shows the difference between the density of bosons and that for fermions with symmetrization effects accounted for by only nearest neighbor transpositions (dimers). The phase space simulations may again be described as quantitatively correct. Evidently the bulk of the symmetrization effects are captured using the classical commutation function, $W=1$.

Figure 4 is for a high density case with lattice spacing $\Delta_{q}=r_{\mathrm{e}}/10$. It can be seen that the exact energy approaches the classical energy from above as the temperature is increased. Without symmetrization effects (main part of figure) the mean field simulations lie between the exact and the classical results, with the pair mean field results lying closer to the exact results than the singlet mean field results across the temperature range shown.

Both mean field results lie below the ground state energy at low temperatures. For example, in this case the ground state energy is $E_{1}/\hbar\omega_{\mathrm{LJ}}=2.6569$, and at $\beta\hbar\omega_{\mathrm{LJ}}=10$, the singlet mean field gives $\langle\hat{\cal H}\rangle/\hbar\omega_{\mathrm{LJ}}=2.336\pm.001$, and the pair mean field gives $2.575\pm.001$. Both are substantially more accurate than the classical simulation result of $0.4001\pm 0.0002$.

The inset to Fig. 4 shows the energy for bosons less that for fermions, calculated by including only nearest neighbor transpositions, $d_{\mathrm{m}}=2$. It can be seen that there is a pole in the exact results for fermions at $\beta\hbar\omega_{\mathrm{LJ}}\approx 0.72$. (This pole disappears when two nearest neighbor dimer transpositions or a cyclic permutation of three consecutive particles, $d_{\mathrm{m}}=4$, are included.)Attard19 It can be seen that the classical, $W=1$, singlet mean field, and pair mean field commutation functions all give this pole for $d_{\mathrm{m}}\leq 2$ at about the same location. There is little to choose between the three at high temperatures; the apparent agreement of the mean field approximations with the exact results for $\beta\hbar\omega_{\mathrm{LJ}}\leq 0.3$ should not be taken seriously because this is about the limit of reliability of the exact results with $l_{\mathrm{max}}=20,000$ in this case. At intermediate and low temperatures, $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 1$, the classical results lie closer to the exact results than do the mean field results. It is difficult to obtain the results for fermions accurately when the denominator passes through zero. In such a regime the neglected higher order terms in the expansion contribute significantly

Figure 5 shows the density profile in this high density case at $\beta\hbar\omega_{\mathrm{LJ}}=1$. It can be seen that there is essentially a single density peak in the center of the system, and that the density spills beyond the wall particles at $q_{0}=0$ and $q_{5}=0.5$. There is little to choose between the three simulation algorithms, at least in the case of monomers (ie. no symmetrization effects, $d_{\mathrm{m}}=0$). Compared to the exact calculations all three simulation algorithms are slightly broader at the peak.

For the case of bosons, Fig. 5B, including nearest neighbor transpositions makes the profile slightly narrower and more sharply peaked. There is again good agreement between the three simulation algorithms, with the pair mean field approach slightly underestimating the height of the density peak. For the case of fermions, Fig. 5C, the exact calculations give a single peak, that is narrower and higher than for bosons. There is no dip in the center of this peak, as one might have expected from Fermi repulsion. (At the lower temperature of $\beta\hbar\omega_{\mathrm{LJ}}=2$, the pair mean field profile shows a bifurcated peak, but the other methods show a single peak similar to here). The classical and the singlet mean field approaches underestimate the height of the peak, with the latter being closer to the exact results than the former. The pair mean field approach overestimates the height of the peak. All three simulation algorithms miss the broad base to the density profile at the walls that is given by the exact calculations. Arguably for reliable results for fermions at this density and temperature, one should go beyond single dimer transpositions.

IV Conclusion

This paper has been concerned with ascertaining the accuracy of a quantum Monte Carlo algorithm for an interacting system. The algorithm is based on a formulation of quantum statistical mechanics in classical phase space and it uses a mean field approximation for the commutation function. For the tests reliable exact results were required as benchmarks, and these were obtained for a one-dimensional harmonic crystal for which the energy eigenvalues and eigenstates can be expressed analytically in closed form. The analytic nature of the exact results for this model allowed up to 20,000 energy levels to be used to establish the benchmark results for the tests. Previous tests of the mean field classical phase space formulation of quantum statistical mechanicsAttard18c used benchmarks established for a one-dimensional Lennard-Jones model with 50 energy eigenvalues obtained numerically. Hernando13 Here it was found that the number of energy levels has a significant effect on the statistical average at higher temperatures, and so the present benchmarks can be relied upon in this regime.

Two versions of the mean field approximation were tested: the singlet version, which has previously been published, Attard18b ; Attard18c and a pair version, which is new here. It was found that the singlet mean field approach was qualitatively correct over the whole temperature (and density and coupling) regime studied. It appeared to be exact in the high temperature limit, and it was generally better than 10% accurate in the low temperature regime in which the system was predominately in the quantum ground state. The pair mean field algorithm significantly improved the accuracy of the singlet algorithm in the intermediate and low temperature regime. In the absence of benchmark results, the difference between the singlet and pair mean field predictions can be used as a guide to the quantitative accuracy of the latter.

There appear to be at least two advantages to the present mean field treatment of the commutation function compared to evaluating it from high temperature expansions. Wigner32 ; Kirkwood33 ; STD2 ; Attard18b First, the mean field expressions remain accurate across the entire temperature regime, including the ground state. Second, algebraically the mean field expressions are relatively simple, and computationally they are easy to implement and efficient to evaluate. In contrast the high temperature expansions rapidly become algebraically complex as higher order terms are included,STD2 ; Attard18b and there are corresponding challenges in their computational implementation and numerical evaluation.

The present paper also explored wave function symmetrization effects. This was at the dimer level, which means the transposition of nearest neighbor particles. It was found, somewhat surprisingly, that combining classical Monte Carlo in classical phase space with the symmetrization function (ie. neglecting the commutation function) in some, but not all, cases gave as good results as those obtained retaining the mean field commutation function. At the highest density studied, some features of the symmetrized system were not captured entirely by the present phase space simulations, particularly in the case of fermions. This suggests that retaining further terms in the symmetrization loop expansion (eg. double dimer and trimer) may be necessary in some cases. It may also be worth reflecting on the underlying philosophy of the mean field approach in the presence of wave function symmetrization.

Finally, a rather interesting conclusion from the present and earlier resultsAttard18c is that the classical component dominates, not just in the high temperature limit but even in the quantum ground state (for structure, not energy). Of course quantum effects are not entirely absent, and when present these are captured by the present mean field commutation function and also the symmetrization function, but it is clear that these truly are a perturbation on the classical prediction. This underscores the utility of treating real world condensed matter systems via quantum statistical mechanics formulated as an integral over classical phase space, rather than formulating it as a sum over quantum states, or by parameterizing the wave function.