# Comment on “Critical temperature of trapped hard-sphere Bose gases”

## Abstract

In this comment, I discuss a recent path-integral Monte Carlo calculation by Pearson, Pang, and Chen ( Phys. Rev. A 58, 4796 (1998)). For bosons with a small hard-core interaction in a harmonic trap, the authors find a critical temperature which does not change with respect to the non-interacting gas. The calculation suffers from a serious discretization error of the many-particle density matrix.

In a recent article, Pearson, Pang, and Chen  studied hard-core bosons in a harmonic trap using Path-Integral Quantum Monte Carlo (PIMC) methods. The authors conclude that the critical temperature of this system for small hard-core diameters $`a`$ does not notably change with respect to the non-interacting gas. These findings are difficult to reconcile both with previous Monte Carlo work  and with mean-field calculations , which agree very nicely with each other. However, the calculation of Pearson, Pang, and Chen is flawed by a discretization error which affects the many-particle density matrix for small $`a`$.

In , the density matrix $`\rho (R,R^{},\beta )`$ at inverse temperature $`\beta `$ is reduced to an integral over high-temperature matrices $`\rho (R,R^{},\tau )`$ with $`\tau =\beta /M`$, and $`M`$ the number of ‘slices’

$$\rho (R,R^{},\beta )=𝑑R_2\mathrm{}𝑑R_M\rho (R,R_2,\tau )\mathrm{}\rho (R_M,R^{},\tau )$$
(1)   
(cf  for notation, $`R=(r_1,\mathrm{},r_N)`$ is the $`3N`$-dimensional vector of particle coordinates, with $`N`$ the number of particles; I use natural units with mass $`m=1`$, oscillator length $`a_0=1`$, $`\mathrm{}=1`$).

To evaluate $`\rho (R_i,R_{i+1},\tau )`$, the authors use the so-called “primitive approximation” (cf ) $`\rho (R_i,R_{i+1},\tau )\rho _P(R_i,R_{i+1},\tau )`$ with

$$\rho _P(R_i,R_{i+1},\tau )=c\mathrm{exp}(\frac{(R_iR_{i+1})^2}{2\tau }\tau V(R_i))$$
(2)   
For the present discussion, the harmonic potential is unessential and I will only consider the hard-core term.

The interaction $`V(r)`$ sets a small length scale, $`a`$, and the primitive approximation is justified only if, between $`R_i`$ and $`R_{i+1}`$, the paths fluctuate by much less than this scale. Under this condition we can conclude that, if paths overlap neither at slice $`i`$ nor at slice $`i+1`$, they will not collide in between. The spatial fluctuations are determined by the kinetic energy term in Eq. (2), and one arrives at the condition

$$\tau (r_ir_i^{})^2a^2$$
(3)   
In their work, Pearson, Pang and Chen consider, e. g. a value $`a=0.02`$. It follows from Eq. (3) that $`\tau `$ must satisfy $`\tau 0.0004`$. For a typical value of $`\beta =0.2`$, this implies $`M500`$. In ref. , $`M10`$ is used. The effect of the interaction is therefore largely underestimated and important deviations from the case of noninteracting bosons are not picked up.

The condition $`\tau a^2`$ renders the use of the primitive approximation unpractical for large systems of weakly interacting bosons. However, it is very instructive to check the estimate Eq. (3) for a single pair of particles in the setting of ref. . To do this, I have generated a large number of $`3`$dimensional Gaussian paths directly from the Lévy construction (cf ) for various values of $`M`$. These paths sample Eq. (1) for $`V=0`$. $`\rho (R,R^{},\beta )`$ is obtained by means of a correction factor, which differs from unity by the probability $`P`$ of pairs of paths to collide at some slice $`i`$.

As an example, I consider the diagonal density matrix for two particles separated by a distance $`|r|`$ ($`R=(0,r)`$, $`R^{}=(0,r)`$). In fig. 1, the numerically determined probability $`P`$ is plotted for $`\beta =0.2`$ as a function of $`|r|`$ for values of $`\tau =0.02`$ ($`M=10`$) to $`\tau =2\times 10^6`$ ($`M=100,000`$). As predicted in Eq. (3), we need large values of $`M`$ to obtain a converged result. In the limit $`M\mathrm{}`$, $`P`$ can also be computed directly from the eigenfunctions of the hard-sphere potential (cf ), as presented in fig. 1. This shows that the interaction of two particles need not be calculated with Monte Carlo methods. This computation is at the heart of the recent successful PIMC calculations for weakly interacting bosons, both in the trap   and in free space .