text:
# Phase diffusion of Bose-Einstein condensates in a one-dimensional optical lattice

## I Introduction

The remarkable experimental realization of Bose-Einstein condensation has generated intensive theoretical and experimental investigations for weakly interacting Bose gases. Although much progress has been made, one of the challenging problems to be solved in both theory and experiment is how to detect the phase diffusion of a Bose-Einstein condensate (BEC). The study of the phase diffusion is of crucial importance because the number of particles realized in all BEC experiments is finite. Within the formalism of the mean-field theory, it is known that the condensate can be well described by a macroscopic wave function (or called an order parameter) with a single phase . The development of the phase at low temperature is governed by the Gross-Pitaevskii (GP) equation. However, for a finite system and hence for trapped Bose gases, although a macroscopic wave function describing a BEC still exists, as is now firmly established experimentally, the phase of the macroscopic wave function will undergo a quantum diffusion process as time develops. Because the phase diffusion reflects directly the coherent nature of a condensate, the study of the phase diffusion of BECs is not only fundamental but also important for applications such as the realization of an atom laser.

For a condensate confined in a trap, the development of its phase is determined by the chemical potential of the system. Due to quantum and thermal fluctuations, the phase diffusion is described by the fluctuations of the chemical potential. Even at initial state the condensate has an exact single phase, the phase will become unpredictable due to the fluctuations of the condensate. After the realization of Bose-Einstein condensation, much theoretical attention has been paid to the phase diffusion of BECs. Because for a single condensate it is hard to observe the effect of phase diffusion based on the density distribution of a Bose-condensed gas, in the experiment conducted by the JILA group two expanding condensates are used to produce an interference pattern and hence the phase diffusion of the condensates can be detected. It was found that in this experiment there is no detectable phase diffusion on the time scale up to $`100`$ $`\mathrm{ms}`$.

Recently, the BECs in optical lattices have attracted increasing attention (see Refs  and references therein). Trapping atoms in optical lattices can provide very effective and powerful means for controlling the properties of a BEC, especially the quantum phase effects on a macroscopic scale such as the coherence of matter waves . In the present work, we investigate the effect of the phase diffusion of a BEC confined in the combined potential of a harmonically magnetic trap and a one-dimensional (1D) optical lattice. In the presence of the 1D optical lattice, there are an array of disk-shaped subcondensates. For each subcondensate, i. e. the BEC formed in the corresponding potential well induced by the optical lattice, there exist phase fluctuations which increase very slowly with the development of time. For high enough potential-well intensity, the phase fluctuations for different subcondensates are independent from each other. But when the potential-well intensity is not high, there will be a correlation for the subcondensates in different wells and hence the situation may be quite different. One expects that the phase diffusion of the BEC in the whole optical lattice can display notable new characters because of the mixing effect of the phase diffusion appearing in different subcondensates.

The paper is organized as follows. Section II is devoted to the phase diffusion of a single 2D (as well as quasi-2D) condensate for the temperature far below the critical temperature. In Sec. III, the phase diffusion of an array of disk-shaped subcondensates induced by a 1D optical lattice is investigated based on the interference pattern resulting from a superposition of expanding subcondensates in different potential wells for different holding time of the optical lattice. Sec. IV contains a discussion and summary of our results.

## II Phase diffusion of a quasi-2D condensate

We consider an array of subcondensates formed in the combined potential consisting of a harmonically magnetic trap and a periodic optical lattice. The combined potential is described by :

$$V=\frac{1}{2}m\omega _x^2x^2+\frac{1}{2}m\omega _{}^2\left(y^2+z^2\right)+sE_R\mathrm{sin}^2\left(\frac{2\pi x}{\lambda }\right),$$
(1)
where $`\omega _x`$ and $`\omega _{}`$ are respectively the axial and radial frequencies of the magnetic trap. The last term represents the potential due to the presence of the optical lattice, with $`\lambda `$ denoting the wavelength of the retroreflected laser beam, and $`s`$ being a factor characterizing the intensity of the optical potential which increases when increasing the intensity of the laser beam. In addition, $`E_R=2\pi ^2\mathrm{}^2/(m\lambda ^2)`$ is the recoil energy of an atom absorbing one photon. For the optical lattice created by the retroreflected laser beam, $`d=\lambda /2`$ is the distance between two neighboring wells induced by the optical lattice.

The optical potential can be approximated as $`V_{eff}=_k\frac{1}{2}m\stackrel{~}{\omega }_x^2(xkd)^2`$ with $`\stackrel{~}{\omega }_x=2\sqrt{s}E_R/\mathrm{}`$ . As in Ref. we assume that for the parameters appearing in the potential (1) the condition $`\omega _x<<\omega _{}<<k_BT<\stackrel{~}{\omega }_x`$ is satisfied. Thus the magnetic trap is a slowly-varying 3D potential which confines the atoms harmonically to overall cigar-shaped distribution. In the 1D optical lattice the atoms are confined to an array of 2D planes. As a result the whole condensate consists of an array of quasi-2D subcondensates like disks in shelf. Such array of quasi-2D subcondensates has been realized in the experiment.

The investigation of the phase diffusion of the above-mentioned BEC can be divided into two steps. The first step is to consider the phase diffusion of a single disk-shaped BEC. For this aim we develop a general approach to calculate the probability distribution of the phase for a condensate confined in a quasi-2D harmonic trap. For the temperature far below the critical temperature $`T_c`$, the condensate can be described by the order parameter with a phase factor $`\varphi \left(t\right)`$:

$$\mathrm{\Phi }(𝐫,t)=\mathrm{\Phi }_0\left(𝐫\right)e^{i\varphi \left(t\right)}.$$
(2)
The differential of the phase takes the form

$$d\varphi =\mu (N_0,T)dt/\mathrm{},$$
(3)
where $`\mu (N_0,T)`$ is the chemical potential of the system. We see that the particle-number fluctuations of the condensate yield fluctuations in the chemical potential, and hence lead to the phase diffusion of the condensate. Assuming that the mean ground state occupation number is $`N_0`$, in the case of the particle-number fluctuations $`\delta N_0<<N_0`$, the average phase of the condensate is then given by

$$\varphi (N_0,t)=\mu (N_0,T)t/\mathrm{}.$$
(4)
The phase diffusion of the condensate can be described by considering the phase difference $`\mathrm{\Delta }\varphi \left(t\right)=\varphi \left(t\right)\varphi (N_0,t)`$. From Eqs. (3) and (4), it is straightforward to obtain a differential equation on $`\mathrm{\Delta }\varphi \left(t\right)`$:

$$\frac{d\mathrm{\Delta }\varphi \left(t\right)}{dt}=\frac{\mu (N_0,T)}{N_0}\mathrm{\Delta }N_0\left(t\right)/\mathrm{}.$$
(5)
In the above expression, $`\mathrm{\Delta }N_0\left(t\right)`$ represents the fluctuations of the ground state occupation number around $`N_0`$ at time $`t`$. Different from the situation considered by Jin et al where a small time-dependent perturbation is used artificially to create selected collective excitations in a BEC, the collective excitations considered here are created spontaneously from the condensate due to quantum fluctuations. Thus $`\mathrm{\Delta }N_0\left(t\right)`$ changes with the development of time due to the creations and annihilations of various collective modes.

For the temperature far below the critical temperature, the collective excitations spontaneously created from the condensate play a dominant role in the fluctuations of the condensate . We therefore pay attention here to the phase diffusion contributed by the collective excitations of the system. For a pure 2D Bose-condensed gas, the collective mode is labelled by two parameters $`\{n,l\}`$. Here, the parameter $`l`$ labels the angular momentum of the excitation. In this case, one has $`\mathrm{\Delta }N_0\left(t\right)=\mathrm{\Sigma }_{nl}\mathrm{\Delta }N_{nl}\left(t\right)`$. Assuming that $`\mathrm{\Delta }\varphi _{nl}\left(t\right)`$ represents the phase diffusion due to the collective mode $`nl`$, we have

$$\frac{d\mathrm{\Delta }\varphi _{nl}\left(t\right)}{dt}=\frac{\mu (N_0,T)}{N_0}\mathrm{\Delta }N_{nl}\left(t\right)/\mathrm{}.$$
(6)
Thus the overall phase fluctuations of the condensate read $`\mathrm{\Delta }\varphi \left(t\right)=\mathrm{\Sigma }_{nl}\mathrm{\Delta }\varphi _{nl}\left(t\right)`$.

In the above equation, $`\left|\mathrm{\Delta }N_{nl}\right|`$ can be regarded as the mean number of atoms corresponding to the collective mode $`nl`$. Based on the Bogoliubov theory , $`\left|\mathrm{\Delta }N_{nl}\right|`$ is given by

$$\left|\mathrm{\Delta }N_{nl}\right|=\left(u_{nl}^2(𝐫)𝑑V+v_{nl}^2(𝐫)𝑑V\right)f_{nl}+v_{nl}^2(𝐫)𝑑V,$$
(7)
where $`f_{nl}=1/\left(\mathrm{exp}\left(\beta \epsilon _{nl}\right)1\right)`$ is the number of the collective excitations presenting at thermal equilibrium, $`\beta =1/(k_BT)`$ ($`T`$ is temperature), and $`\epsilon _{nl}=\mathrm{}\omega _{nl}`$ is the energy spectrum of the collective mode $`nl`$. In the above equation, the quantity $`v_{nl}^2(𝐫)𝑑V`$ describes the effect of quantum depletion, which does not vanish even at $`T=0`$ . In Eq. (7), $`u_{nl}(𝐫)`$ and $`v_{nl}(𝐫)`$ are determined by the following coupled equations:

$$(\frac{\mathrm{}^2}{2m}^2+V_{}\left(𝐫\right)\mu +2gn\left(𝐫\right))u_{nl}(𝐫)$$

$`+gn_0\left(𝐫\right)v_{nl}(𝐫)=\epsilon _{nl}u_{nl}(𝐫),`$

$$(\frac{\mathrm{}^2}{2m}^2+V_{}\left(𝐫\right)\mu +2gn\left(𝐫\right))v_{nl}(𝐫)$$

$$+gn_0\left(𝐫\right)u_{nl}(𝐫)=\epsilon _{nl}v_{nl}(𝐫),$$
(8)
where $`n\left(𝐫\right)`$ and $`n_0\left(𝐫\right)`$ are the density distributions of the Bose gas and condensate, respectively. $`V_{}\left(𝐫\right)=\frac{1}{2}m\omega _{}^2(y^2+z^2)`$ is the 2D harmonic potential confining the Bose gas in $`y`$ and $`z`$ directions (the second term on the right hand side of (1) ), while $`g`$ is the coupling constant describing the interatomic interaction in the condensate.

By a straightforward calculation the excitation frequency of the collective mode $`nl`$ is given by

$$\omega _{nl}=\omega _{}\left(2n^2+2n\left|l\right|+2n+\left|l\right|\right)^{1/2},$$
(9)
while $`\left|\mathrm{\Delta }N_{nl}\right|`$ reads

$$\left|\mathrm{\Delta }N_{nl}\right|=\frac{\left(N_0mg\right)^{1/2}\beta _{nl}}{2\sqrt{\pi }\mathrm{}\sqrt{2n^2+2n\left|l\right|+2n+\left|l\right|}},$$
(10)
where

$`\beta _{nl}={\displaystyle _0^1}\left(1x^2\right)x\left(H_{nl}\right)^2𝑑x/{\displaystyle _0^1}x\left(H_{nl}\right)^2𝑑x,`$ (11)
$`H_{nl}\left(x\right)=x^{|l|}{\displaystyle \underset{j=0}{\overset{n}{}}}b_jx^{2j}`$ (12)
with $`b_0=1`$ and the coefficients $`b_j`$ satisfying the recurrence relation $`b_{j+1}/b_j=\left(4j^2+4j+4j|l|4n^24n4n|l|\right)/\left(4j^2+4j|l|+8j+4|l|+4\right)`$.

As mentioned above, the collective excitations considered here are created spontaneously due to the quantum fluctuations of the condensate. They can disappear as time develops and hence have a finite longevity. The longevity can be estimated based on the time-energy uncertainty relation. For the collective mode $`nl`$, its longevity is $`\tau _{nl}=1/\omega _{}\left(2n^2+2n\left|l\right|+2n+\left|l\right|\right)^{1/2}`$. Thus at time $`t>>\tau _{nl}`$, there is a series of $`i=t/\tau _{nl}`$ successive creations and annihilations of the collective mode $`nl`$. Assuming that at the time $`t`$ the number of the collective modes created (annihilated) in the condensate is $`i_{cre}`$ ($`i_{ann}`$). Thus one has $`i=i_{cre}+i_{ann}`$. In this situation, the desired probability is given by the binomial expression

$$P_i\left(j\right)=\left(\frac{1}{2}\right)^i\frac{i!}{\left\{\frac{1}{2}(i+j)\right\}!\left\{\frac{1}{2}(ij)\right\}!},$$
(13)
where $`j=i_{cre}i_{ann}`$. For $`t>>\tau _{nl}`$, using the Stirling’s formula, the asymptotic form of the probability distribution $`P_{nl}(\mathrm{\Delta }\varphi _{nl},t)`$ of the phase diffusion due to the collective mode $`nl`$ takes the form

$$P_{nl}(\mathrm{\Delta }\varphi _{nl},t)=\frac{1}{\mathrm{\Gamma }_{nl}\sqrt{2\pi t/\tau _{nl}}}\mathrm{exp}\left(\frac{\tau _{nl}\left(\mathrm{\Delta }\varphi _{nl}\right)^2}{2\mathrm{\Gamma }_{nl}^2t}\right),$$
(14)
where the dimensionless parameter $`\mathrm{\Gamma }_{nl}`$ is given by

$$\mathrm{\Gamma }_{nl}=\frac{\mu (N_0,T)}{N_0}\frac{\left|\mathrm{\Delta }N_{nl}\right|\tau _{nl}}{\mathrm{}}.$$
(15)
From the Gaussian distribution given by Eq. (14), the phase fluctuations due to the collective mode $`nl`$ have the form:

$$\left(\delta ^2\varphi \right)_{nl}=\mathrm{\Gamma }_{nl}^2t/\tau _{nl}.$$
(16)
Assuming that the collective modes with different $`nl`$ are created and annihilated independently, one obtains that the probability distribution of the overall phase difference $`\mathrm{\Delta }\varphi \left(t\right)`$ is still a Gaussian distribution function. The probability distribution of $`\mathrm{\Delta }\varphi `$ is then given by

$$P(\mathrm{\Delta }\varphi ,t)=\sqrt{\frac{B}{\pi }}\mathrm{exp}\left(B\left(\mathrm{\Delta }\varphi \right)^2\right).$$
(17)
Based on the general probability theory, the parameter $`B`$ can be obtained through the relation $`\delta ^2\varphi =1/2B=\mathrm{\Sigma }_{nl}\left(\delta ^2\varphi \right)_{nl}`$. From the formulas (10), (15), (16) and the chemical potential $`\mu (N_0,T)=\left(N_0mg/\pi \right)^{1/2}\omega _{}`$ for 2D harmonic trap, we obtain the overall phase fluctuations $`\delta ^2\varphi `$ as

$$\delta ^2\varphi =\frac{m^2g^2\omega _{}t}{4\pi ^2\mathrm{}^4}\mathrm{{\rm Y}},$$
(18)
where the coefficient $`\mathrm{{\rm Y}}`$ takes the form

$$\mathrm{{\rm Y}}=\underset{nl}{}\frac{\beta _{nl}^2}{\left(2n^2+2n\left|l\right|+2n+\left|l\right|\right)^{3/2}}\left[f_{nl}+\frac{1}{2}\right]^2,$$
(19)
where $`\beta _{nl}`$ has been given in Eq. (11). From Eq. (18) we see that when the longevity of the collective excitations is taken into account, the phase fluctuations are proportional to the time $`t`$, rather than $`t^2`$ which is obtained when the collective excitations are regarded as quite stable.

As mentioned above, if $`\omega _x<<\omega _{}<<k_BT<\stackrel{~}{\omega }_x`$, the condensate confined in each well has a quasi-2D nature. Note that for the array of the subcondensates induced by the optical lattice with enough depth, the overlap between the subcondensates occupying different wells can be omitted. This means that the behavior of the condensate in each well has a local property. For each well the coupling constant $`g`$ is given by $`g2\sqrt{2\pi }\mathrm{}^2a_s/(m\stackrel{~}{l_x})`$ , which is fixed by a $`s`$-wave scattering length $`a_s`$ and the oscillator length $`\stackrel{~}{l_x}=\left(\mathrm{}/m\stackrel{~}{\omega }_x\right)^{1/2}`$ in the $`x`$-direction. Thus from Eq. (18) the subcondensate in each well displays the phase fluctuations

$$\delta ^2\varphi =\frac{2\mathrm{{\rm Y}}\omega _{}t}{\pi }\left(\frac{a_s}{\stackrel{~}{l_x}}\right)^2.$$
(20)
From the above result, one has $`\delta ^2\varphi \omega _{}\stackrel{~}{\omega }_x`$. We see that the confinement induced by the optical lattices has the effect of increasing the phase fluctuations. This is a natural character by considering the fact that the confinement has the effect of increasing the particle-number fluctuations .

## III Effect of phase diffusion on the interference pattern of expanding array of subcondensates

We now go to the second step, i. e. to discuss the phase diffusion of an array of subcondensates confined in the combined potential of the harmonically trapping potential and 1D optical lattice. To illustrate clearly the role of the phase diffusion in the interference pattern of the expanding subcondensates, we consider here an experiment scheme which can be realized in future experiment. Firstly, a cigar-shaped condensate is formed in a magnetic trap with the harmonic angular frequencies $`\omega _x=`$ $`2\pi \times 9`$ $`\mathrm{Hz}`$ and $`\omega _{}=2\pi \times 92`$ $`\mathrm{Hz}`$ at a temperature far below the critical temperature. The atoms are transferred into the lattice potential by increasing the power of the laser beam so that the harmonic angular frequency of the well is approximated to be $`2\pi \times 14`$ $`\mathrm{kHz}`$ in $`x`$direction. For $`\stackrel{~}{\omega }_x=2\pi \times 14`$ $`\mathrm{kHz}`$, based on the well-known Bose-Hubbard model, the tunnelling time between neighboring subcondensates is of the order of $`1`$ ms. In this situation, the array of subcondensates should be regarded as fully coherent and there is no independent random phase for the subcondensates in different wells. Although there is still a phase diffusion for the whole condensate, one can not observe this through the interference pattern. To observe the effect of the phase diffusion, one can rapidly increase the lattice potential depth to a value of $`\stackrel{~}{\omega }_x=`$ $`2\pi \times 52`$ $`\mathrm{kHz}`$ within a time (for example $`50`$ $`\mu \mathrm{s}`$) much smaller than the tunnelling time. In this situation, the phase of every subcondensate is the same after the lattice potential ramps to its final strength. For the lattice potential with $`\stackrel{~}{\omega }_x=`$ $`2\pi \times 52`$ $`\mathrm{kHz}`$, the tunnelling time between neighboring wells is of the order of $`1.5`$ $`\mathrm{s}`$. Thus, for the holding time $`t_h`$ of the optical lattice being much smaller than $`1.5`$ $`\mathrm{s}`$, the phase diffusion of the subcondensates in different wells can be regarded as independent from each other. In this situation, there is a random phase for each subcondensate with the development of the holding time and this will lead to an obvious effect on the interference pattern after the combined potential is switched off.

For the experimental parameters considered here, we have $`\delta ^2\varphi =1.60\times t_h`$. Based on the result given by Eq. (20), the probability distribution of the phase for the subcondensate in each well can be obtained through Eq. (17). For the interference pattern to wash out completely, the time scale of the holding time is obtained by using $`\delta ^2\varphi =\pi ^2`$. In the experimental parameters used here, for the subcondensate in each well $`\delta ^2\varphi `$ is still much smaller than $`\pi ^2`$ even when the holding time $`t_h`$ of the optical lattice is $`1\mathrm{s}`$. Due to the fact that there are many subcondensates induced by the optical lattice, however, we will show that there is a significant effect of the phase diffusion to the interference pattern after the optical lattice is switched off.

To observe a notable effect of the phase diffusion, we consider the case when only the optical lattice is switched off after the holding time $`t_h=0.2`$ $`\mathrm{s}`$. Assume $`t_0`$ is the time after the optical lattice is switched off. For the subcondensate confined in the $`k`$th well of the optical lattice, applying the propagator method used in Ref. it is easy to get the normalized wave function $`\phi _k(x,t_0)`$ after only the optical lattice is switched off:

$`\phi _k(x,t_0)=A_n\sqrt{{\displaystyle \frac{1}{\mathrm{sin}\omega _xt_0\left(\mathrm{ctg}\omega _xt_0+i\gamma \right)}}}\mathrm{exp}[iR_k(\delta \varphi ,t_h)]`$

$$\times \mathrm{exp}\left[\frac{\left(kd\mathrm{cos}\omega _xt_0x\right)^2}{2\sigma ^2\mathrm{sin}^2\omega _xt_0\left(\mathrm{ctg}^2\omega _xt_0+\gamma ^2\right)}\right]$$

$$\times \mathrm{exp}\left[\frac{i\left(kd\mathrm{cos}\omega _xt_0x\right)^2\mathrm{ctg}\omega _xt_0}{2\gamma \sigma ^2\mathrm{sin}^2\omega _xt_0\left(\mathrm{ctg}^2\omega _xt_0+\gamma ^2\right)}\right]$$

$$\times \mathrm{exp}\left[\frac{i\left(x^2\mathrm{cos}\omega _xt_0+k^2d^2\mathrm{cos}\omega _xt_02xkd\right)}{2\gamma \sigma ^2\mathrm{sin}\omega _xt_0}\right],$$
(21)
where $`A_n=1/\pi ^{1/4}\sigma ^{1/2}`$ is a normalization constant, the dimensionless parameter $`\gamma =\mathrm{}/m\omega _x\sigma ^2`$ with $`\sigma `$ the width of the subcondensate in each well . The factor $`R_k(\delta \varphi ,t_h)`$ represents a random phase due to the phase diffusion of the subcondensate in the $`k`$th well. In numerical calculations, $`R_k(\delta \varphi ,t_h)`$ is generated according to the probability distribution function $`P(\mathrm{\Delta }\varphi ,t_h)`$ (see (17) ). Then we obtain the density distribution in $`x`$direction as follows:

$$n(x,t_0)=\mathrm{\Xi }\left|\underset{k=k_M}{\overset{k_M}{}}\left(1\frac{k^2}{k_M^2}\right)\phi _k(x,t_0)\right|^2,$$
(22)
where $`\mathrm{\Xi }=15N_{all}k_M^3/\left(16k_M^41\right)`$. $`N_{all}`$ is the total number of atoms in the array of condensates and $`2k_M+1`$ is the total number of the subcondensates induced by the optical lattice. In the present work, by using the experimental result in , $`k_M`$ is chosen to be $`100`$.

It is obvious that the density distribution $`n(x,t_0)`$ at $`x=0`$ reaches the maximum value at time $`t_l=\left(2l1\right)\pi /2\omega _x`$, where $`l`$ is a positive integer. Displayed in Fig. 1 is $`n(x,t_0)`$ (with $`N_{all}\mathrm{\Xi }A_n^2`$ as a unit) at $`t_l`$ with the holding time $`t_h=0.2`$ $`\mathrm{s}`$. We see that the phase diffusion of the array of the subcondensates makes the central density of the interference pattern decrease significantly. Displayed in Fig. 2 is the density distribution $`n\left(x=0,t_0=t_l\right)`$ versus the holding time $`t_h`$. A exponential damping of the central density is clearly shown in the figure.

It is straightforward to calculate the interference pattern of the expanding array of BECs when the magnetic trap and optical lattice are both switched off. Shown in Fig. 3 is the interference pattern for the holding time $`t_h=0.2`$ $`\mathrm{s}`$ after the combined potential is switched off for $`30`$ $`\mathrm{ms}`$. In comparison with the case of $`t_h=0`$ $`\mathrm{s}`$ (dashed line in Fig. 4), the strong noise in Fig. 3 shows that the phase diffusion has an important effect on the interference pattern of the array of subcondensates. The solid line in Fig. 4 displays the average density distribution of twenty ensembles for the holding time $`t_h=0.2`$ $`\mathrm{s}`$. We see that the maximum density of the central and side peaks decreases significantly due to the presence of the phase diffusion. Generally speaking, there are four parameters $`\omega _{}`$, $`\stackrel{~}{\omega }_x`$, $`k_M`$, $`t_h`$ which are related closely to the effect of the phase diffusion. The present work shows that increasing these parameters will contribute to the observation of the phase diffusion for BEC in optical lattices.

## IV Discussion and summary

Based on the analysis of the collective excitations, we have investigated the phase diffusion at low temperature in an array of subcondensates formed in a combined potential consisting of harmonically magnetic trap and 1D optical lattice. We have shown that the effect of phase diffusion on the interference pattern is evident and thus can be measured from the expanding subcondensates after the optical lattice is switched off. We point out that the maximum density of the interference pattern decreases significantly due to the mixing of the phase diffusion appearing in different subcondensates formed by the optical lattice. This suggests a way to detect experimentally the notable phase diffusion of BEC. For the case of BEC in 3D optical lattice, however, it seems that the theory developed here can not be used due to the fact that there are only an average atom number of up to $`2.5`$ atoms per lattice site .

Note that a recent experiment has been conducted in Ref. where the interference pattern is measured when only the optical lattice is switched off. Clearly this type of experiment can be used to test the theoretical predictions provided in this paper. We stress that for a single subcondensate the phase fluctuations increase very slowly with the development of the time. However, the phase diffusion plays an important role in the interference pattern of the expanding subcondensates, as shown in the present work.

## Acknowledgments

We acknowledge valuable discussions with L. You. This work was supported by Natural Science Foundation of China under grant Nos. 10205011 and 10274021.


format: text

text:
# Maximum Entropy approach to a Mean Field Theory for FluidsPresented at MaxEnt 2002, the 22th International Workshop on Bayesian Inference and Maximum Entropy Methods (August 3-7, 2002, Moscow, Idaho, USA)

## 1 Introduction

The method of Maximum Entropy (ME) has been designed to solve the general problem of updating from a prior probability distribution to a posterior distribution when new information in the form of a constraint becomes available . The method reflects the deep conviction that one should not change one’s mind frivolously, that whatever was learned in the past is important. Indeed, this concern for prior information is the only guarantee that the new information itself will not be easily discarded in the future. The chosen posterior distribution should coincide with the prior as closely as possible: one should only update those aspects of one’s beliefs for which hard new evidence has been supplied -.

The purpose of this work is to explore the use of the ME method as a means to tackle problems of a very different kind. The assignment of a probability distribution that faithfully codifies our state of knowledge is only a first step towards answering our questions about a physical system. There is a second and crucial step that requires us to extract the desired answers from the probability distribution. It is not enough to write down a probability distribution, we must also be able to read it, to figure out what it means, what it implies. We must be able to calculate the expected values of whatever quantities happen to be of interest to us. This is a difficult practical problem for which no solution of universal applicability is known and, accordingly, a wide variety of different techniques, both analytical and numerical, have been developed.

The approach we explore here consists of approximating the ‘exact’ distribution by a more tractable one. Then, the question is: given a family of trial tractable distributions, which one should we select? Which distribution within the given family best approximates the exact one? When phrased in this way it is clear that this is the kind of question that can be tackled using the method of ME: one must select the trial distribution that reflects the least misrepresentation of our state of knowledge. More specifically, the chosen distribution should be that which maximizes the appropriate relative entropy.

The purpose of this work is to explore the use of ME as a method to generate approximations. As a simple illustration, in section 2, we apply the ME method to canonical Gibbs-Boltzmann distributions. The resulting variational technique, when expressed in terms of the free energies, is recognized as the Bogoliuvob variational principle . This is a well known technique for generating mean field approximations for discrete systems on a lattice. What is perhaps not as widely known is that the Bogoliuvob principle is just another application of the ME method (however, see ).

In the bulk of the paper we explore the use of the ME method to generate mean field (MF) approximations for classical fluids. Ever since the pioneering work of van der Waals many different versions of MF theories have been proposed and it is not always clear what they have in common – except perhaps for neglecting correlation effects. It appears that none of these versions have been derived as a systematic application of the ME method.

The study of classical fluids is an old and mature field. There exist extensive treatments in many excellent books and reviews . Our goal here is not so much to contribute to the study of the fluids themselves but rather to use this well-explored but still non-trivial field as a testing ground for further development of the ME method. In future work we intend to build on our current results and extend the method to generate other more sophisticated approximations involving canonical transformations to collective coordinates; only then do we expect to obtain significantly new contributions to the theory of fluids.

The study of both homogeneous and inhomogeneous fluids is currently best approached using density functional methods  which have been noted to show a remarkable formal similarity to statistical thermodynamics . In section 3 we show that the resemblance goes beyond a mere accident: density functional methods constitute a generalization of the Gibbs grand canonical ensemble to the case of a spatially dependent chemical potential. The resulting ‘density’ distribution is, however, too complicated for actual calculations and approximations are needed. In section 4 the MF approximation for the density, various thermodynamic potentials and pair distribution functions is derived. Following an idea due to Percus , a second, somewhat improved MF estimate of the two-particle distribution of a homogeneous fluid is obtained by expressing it in terms of the one-particle distribution of the fluid placed in a suitable external potential. An explicit comparison with experimental data for Argon is carried out in section 5. This serves both to determine the parameters of the Lennard-Jones potential used to model the interaction between Argon atoms and to explore the limitations of the MF approximation. Conclusions and final comments appear in section 6.

## 2 Maximum Entropy as an approximation technique

Let the microstates of a system (for example, its location in phase space or perhaps the values of spin variables) be labelled by coordinates $`q`$, and let us assume that the probability that the system is in a microstate within a particular range $`dq`$ is given by the canonical distribution

$$P(q)dq=\frac{e^{\beta H(q)}}{Z}dq\text{where}Z=e^{\beta F}=𝑑qe^{\beta H(q)}.$$
(1)
In principle knowing the distribution $`P`$ we should be able to answer any question we might care to ask but in practice complications arise. Any realistic Hamiltonian $`H(q)`$ will usually include complicated non linear interactions that make the task of calculating expected values or most other integrals over $`P`$ impossibly difficult.

A possible way around this problem consists of replacing the ‘exact’ distribution $`P`$ by another more tractable distribution $`P_0`$ that approximately represents the same state of knowledge. This is perhaps too much to ask for. What we really want is a distribution $`P_0`$ that approximately represents those aspects of the information in $`P`$ that happen to be relevant to those very few, very specific questions that we are actually asking.

The practical problem of choosing a $`P_0`$ is dealt with in two steps. First we search for a family of trial distributions that are simple enough that actual calculations are feasible and that we have some reason to suspect might codify the appropriate relevant information. This step is the difficult one because there is no known systematic procedure to carry it out; it is a matter trial and error guided by intuition. The second and easier step is to select the one member within the trial family that best resembles the exact $`P`$. This step is easier because there exists a systematic, mechanical procedure to be followed. As discussed in  (and other references therein) there is a unique selection criterion satisfying natural desiderata of consistency and objectivity: the distribution to be selected is that which maximizes the entropy of $`P_0`$ relative to $`P`$,

$$S[P_0,P]=𝑑qP_0(q)\mathrm{log}\frac{P_0(q)}{P(q)}.$$
(2)
The success of the whole method hinges crucially on the first step, the choice of the family of trial distributions. To increase the likelihood that important relevant features (*e.g.*, symmetries or dominant interactions) are not left out it is usual – but not necessary – to select trial distributions that could conceivably model an idealized system. Such trial distributions are also canonical but with a modified Hamiltonian $`H_0(q;\alpha )`$ where $`\alpha `$ are parameters labeling each distribution within the family,

$$P_0(q;\alpha )dq=\frac{e^{\beta H_0(q;\alpha )}}{Z_0}dq\text{where}Z_0(\alpha )=e^{\beta F_0(\alpha )}=𝑑qe^{\beta H_0(q;\alpha )}.$$
(3)
Substituting into eq.(2) gives,

$$S[P_0,P]=\beta \left(H_0H_0F_0+F\right),$$
(4)
where $`\mathrm{}_0`$ refers to averages over the trial $`P_0`$. The inequality $`S[P_0,P]0`$, can then be written as

$$FF_0+HH_0_0.$$
(5)
Thus, one selects $`P_0(q;\alpha )`$ by maximizing $`S`$ over all values of $`\alpha `$ or equivalently by minimizing the quantity $`F_0+HH_0_0`$. This alternate form of the variational principle and its use to generate approximations is well known. It is usually associated with the name of Bogoliuvob  and it is the main technique to generate MF approximations for discrete systems of spins on a lattice. What is perhaps not as widely known is that the Bogoliuvob variational principle is just a peculiar application of the ME method.

After this brief formal illustration we are ready to tackle the problem of fluids.

## 3 The canonical ‘density’ ensemble

Many of the most interesting questions we ask when discussing fluids can be addressed once we know the density of the fluid under various experimental conditions. Examples include the transitions between the gas and liquid phase, the structure of the gas-liquid interface, surface tension and capillarity effects and so on. We therefore would like to have a formalism where the (possibly nonuniform) density $`n(r)`$ appears as an explicit variable. It is very easy to do this, just find the probability distribution that maximizes entropy subject to constraints on the expected energy and on the expected density. The resulting distribution is a generalized canonical ensemble we will call the ‘density’ ensemble.

The logic of the previous paragraph may sound unfamiliar and require further comments. When justifying the use of the ME method to obtain the canonical Boltzmann-Gibbs distribution ($`P_qe^{\beta E_q}`$) one might say something like “we seek the minimally biased (*i.e.* maximum entropy) distribution that codifies the information we do possess (the expected energy) and nothing else”. Many authors – even those who fully appreciate the value of the concept of entropy in inductive reasoning – find this justification unsatisfactory. Indeed, they would argue, for example, that the spectrum of black body radiation is what it is independently of our knowledge. We prefer to phrase the objection differently: in most realistic situations the one thing that is not known is the expected energy, and even then one still makes correct predictions by maximizing entropy subject to a constraint on the (unknown) expected energy. But then, how can we justify the imposition of such ‘unknown’ constraints in the ME method?

We propose here that the motivation behind imposing a constraint on the expected energy should not be that this is a quantity we happen to know – because we do not – but rather that we recognize theexpected energy as the quantity we *need* to know. We recognize that in many experimental situations the expected energy, even if unknown, is the crucial quantity that codifies the relevant information and we proceed to theorize as if we knew it. All resulting predictions will contain the temperature as a free unknown parameter. Later, when it comes to actually comparing theoretical predictions with experimental data we expect that part of the data analysis will consist of adjusting the unknown temperature parameter to provide the best fit to the data. In other words, the information on expected energy or temperature is something to be obtained from the experiment itself.

Thus, the motivation behind imposing a constraint on the density is not that it is something we know, but rather that it is what we should know; it is the quantity that codifies information that is very relevant for most questions of interest and therefore it is important that it appear explicitly in the formalism.

The probability that the fluid is composed of $`N`$ particles with positions and momenta within the phase space volume $`dq_N`$ at the location $`q_N=\{p_i,r_i;i=1,\mathrm{},N\}`$ is given by

$$P(q_N;\beta ,\lambda )dq_N$$
(6)
where

$$dq_N=\frac{1}{N!h^{3N}}\underset{i=1}{\overset{N}{}}d^3p_id^3r_i,$$
(7)
and

$$P(q_N;\beta ,\lambda )=\frac{1}{Z}\mathrm{exp}\beta \left[H(q_N)+d^3r\lambda (r)\widehat{n}(r)\right],$$
(8)
where $`\lambda (r)`$ are Lagrange multipliers that enforce the constraint on the expected density $`\widehat{n}(r)`$ at each point in space and the density is

$$\widehat{n}(r)=\underset{i=1}{\overset{N}{}}\delta (rr_i).$$
(9)
$`H(q_N)`$ is the Hamiltonian,

$$H(q_N)=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2m}+\underset{i>j}{\overset{N}{}}u(r_{ij})+\underset{i=1}{\overset{N}{}}v(r_i),$$
(10)
where $`r_{ij}=r_ir_j`$, $`u(r_{ij})`$ is the two-body interparticle potential, and $`v(r)`$ is an external potential. The partition function is

$$Z(T,\lambda )=\underset{N=0}{\overset{\mathrm{}}{}}𝑑q_N\mathrm{exp}\beta \left[H(q_N)+d^3r\lambda (r)\widehat{n}(r)\right]\stackrel{\mathrm{d}ef}{=}e^{\beta \mathrm{\Omega }(T,\lambda )},$$
(11)
where we introduced the thermodynamic potential $`\mathrm{\Omega }(T,\lambda )`$, a function of the temperature $`T1/\beta `$ (units are such that $`k_B=1`$), and a functional of $`\lambda (r)`$.

The interparticle and the external potential energies can be rewritten in terms of the density,

$$\underset{i=1}{\overset{N}{}}v(r_i)=\underset{i=1}{\overset{N}{}}d^3rv(r)\delta (rr_i)=d^3rv(r)\widehat{n}(r),$$
(12)
and

$$U=\underset{i>j}{\overset{N}{}}u(r_{ij})=\frac{1}{2}\underset{ij}{\overset{N}{}}u(r_{ij})=\frac{1}{2}d^3rd^3r^{}u(rr^{})\widehat{n}^{(2)}(r,r^{}),$$
(13)
where $`\widehat{n}^{(2)}(r,r^{})`$ is the two-particle distribution,

$$\widehat{n}^{(2)}(r,r^{})=\underset{ij}{}\delta (rr_i)\delta (r^{}r_j)=\widehat{n}(r)\widehat{n}(r^{})\widehat{n}(r)\delta (rr^{}).$$
(14)
We remark that a constraint on $`\widehat{n}(r)`$ also constrains the expected number of particles

$$N=d^3r\widehat{n}(r).$$
(15)
Therefore the density ensemble is a generalization of the grand-canonical ensemble. Indeed, when $`\lambda (r)`$ is a constant the density distribution reduces to the grand-canonical distribution and $`\lambda `$ coincides with the chemical potential $`\mu `$.

Notice also that the effects of the Lagrange multiplier field $`\lambda (r)`$ are indistinguishable from the effects of the external potential $`v(r)`$. Similarly, when two-particle correlations $`\widehat{n}^{(2)}(r,r^{})`$ carry relevant information we might want to develop a formalism where they appear explicitly. The effects of the corresponding Lagrange multipliers are indistinguishable from those of the two-body potentials $`u(r,r^{})`$. Perhaps this is the ultimate reason why potentials in general are valuable in physics: they establish control over the relevant physical quantities, the density, and the two-particle correlations.

The expected value of the density is

$$\widehat{n}(r)=\frac{\delta \mathrm{\Omega }}{\delta \lambda (r)}\stackrel{\mathrm{d}ef}{=}n(r),$$
(16)
and the entropy is

$$S[P]=\underset{N}{}𝑑q_NP\mathrm{log}P=\frac{\mathrm{\Omega }}{T}$$
(17)
These two equations can be combined into

$$d\mathrm{\Omega }=SdT+d^3rn(r)\delta \lambda (r).$$
(18)
To achieve a formulation in terms of the density $`n(r)`$ rather than the multiplier $`\lambda (r)`$ we consider the Legendre transform

$$\mathrm{\Phi }[T,n]=\mathrm{\Omega }[T,\lambda ]d^3r\lambda (r)n(r),$$
(19)
so that

$$d\mathrm{\Phi }=SdTd^3r\lambda (r)\delta n(r).$$
(20)
The density functional formalism of Hohenberg and Kohn  is founded upon a theorem proving the existence of a functional $`\mathrm{\Phi }[T,n]`$ that is independent of the external potential $`v(r)`$. From the perspective afforded by the ME approach the independence of $`\mathrm{\Phi }[T,n]`$ on the potential $`v(r)`$ is a triviality achieved by construction – through the Legendre transform – and their truly nontrivial insight is the recognition of the density as the correct choice of relevant variable.

## 4 Mean field approximation

### 4.1 MF trial distributions

We want to approximate the ‘exact’ probability distribution by a more tractable trial distribution in which the interparticle interactions $`U`$ are replaced by an interaction with an external potential, the mean field. The problem consists of selecting the best trial distribution, that which best approximates the exact distribution. The drastic approximation being made is the neglect of two-particle correlations induced by the interparticle potential.

The trial MF distribution is

$$P_0(q_N;\beta ,\lambda )=\frac{1}{Z}\mathrm{exp}\beta \left[H_0(q_N)+d^3r\lambda (r)\widehat{n}(r)\right],$$
(21)
where

$$H_0(q_N)=\underset{i=1}{\overset{N}{}}\frac{p_i^2}{2m}+d^3r\left(v(r)+v_0(r)\right)\widehat{n}(r),$$
(22)
where $`v_0(r)`$ is the mean field to be determined.

First we compute several thermodynamic quantities of interest. It is convenient to absorb the external potential $`v(r)`$, the mean field $`v_0(r)`$, and the multiplier field $`\lambda (r)`$ into a single potential $`V(r)`$,

$$V(r)=v(r)+v_0(r)+\lambda (r).$$
(23)
The partition function $`Z_0`$ is

$$Z_0=\underset{N=0}{\overset{\mathrm{}}{}}\frac{1}{N!}\left[\frac{d^3p}{h^3}\mathrm{exp}\frac{\beta p^2}{2m}\right]^N\left[d^3r\mathrm{exp}\beta V(r)\right]^N\stackrel{\mathrm{d}ef}{=}e^{\beta \mathrm{\Omega }_0(\beta ,\lambda )},$$
(24)
so that

$$\mathrm{\Omega }_0[T,\lambda ]=\frac{1}{\beta \mathrm{\Lambda }^3}d^3re^{\beta V(r)}\text{where}\mathrm{\Lambda }=\left(\frac{\beta h^2}{2\pi m}\right)^{1/2}.$$
(25)
The expected density is

$$\widehat{n}(r)_0=\frac{\delta \mathrm{\Omega }_0}{\delta \lambda (r)}=\frac{\delta \mathrm{\Omega }_0}{\delta V(r)}\stackrel{\mathrm{d}ef}{=}n_0(r)$$
(26)
or

$$n_0(r)=\frac{e^{\beta V(r)}}{\mathrm{\Lambda }^3}.$$
(27)
### 4.2 ME optimization of the mean field

Within the family of trial MF distributions the one that is closest to the canonical density ensemble is that which maximizes the relative entropy

$$S[P_0|P]=\underset{N=0}{\overset{\mathrm{}}{}}𝑑q_NP_0(q_N|\beta ,\lambda ,v_0)\mathrm{log}\frac{P_0(q_N|\beta ,\lambda ,v_0)}{P(q_N|\beta ,\lambda )}$$
(28)
Substituting Eqs.(8) and (21) we obtain

$$S[P_0|P]=\beta \left[\mathrm{\Omega }\mathrm{\Omega }_0HH_0_0\right].$$
(29)
(The subscript $`0`$ in $`\mathrm{}_0`$ indicates the averages are computed over the trial distribution $`P_0`$.) Since $`S[P_0|P]0`$, we have

$`\mathrm{\Omega }[T,\lambda ]`$ $``$ $`\mathrm{\Omega }_U[T,\lambda ,v_0]\stackrel{\mathrm{d}ef}{=}\mathrm{\Omega }_0+HH_0_0`$ (30)
$`=`$ $`\mathrm{\Omega }_0+U_0{\displaystyle d^3rv_0(r)n_0(r)},`$ (31)
and maximizing $`S[P_0|P]`$ is equivalent to minimizing $`\mathrm{\Omega }_U`$ over all mean fields $`v_0.`$ The MF approximation to $`\mathrm{\Omega }`$, denoted $`\overline{\mathrm{\Omega }}`$, is defined as the best $`\mathrm{\Omega }_U`$,

$$\mathrm{\Omega }[T,\lambda ]\overline{\mathrm{\Omega }}[T,\lambda ]\stackrel{\mathrm{d}ef}{=}\underset{v_0}{\mathrm{min}}\mathrm{\Omega }_U[T,\lambda ,v_0].$$
(32)
To calculate the potential energy $`U_0`$ use Eqs.(13) and (14)

$$U_0=\frac{1}{2}d^3rd^3r^{}u(rr^{})n_0^{(2)}(r,r^{}),$$
(33)
where

$$n_0^{(2)}(r,r^{})=\widehat{n}^{(2)}(r,r^{})_0=\widehat{n}(r)\widehat{n}(r^{})_0n_0(r)\delta (rr^{}),$$
(34)
furthermore, from Eq.(21) and (25-27)

$$\widehat{n}(r)\widehat{n}(r^{})_0=n_0(r)n_0(r^{})+n_0(r)\delta (rr^{}).$$
(35)
Therefore

$$n_0^{(2)}(r,r^{})=n_0(r)n_0(r^{}),$$
(36)
and the result for $`U_0`$ is

$$U_0=\frac{1}{2}d^3rd^3r^{}u(rr^{})n_0(r)n_0(r^{}).$$
(37)
The best choice of mean field $`v_0(r)`$ is given by

$$0=\frac{\delta \mathrm{\Omega }_U}{\delta v_0(r)}=\frac{\delta \mathrm{\Omega }_U}{\delta V(r)}=\beta n_0(r)\left[v_0(r)d^3r^{}u(rr^{})n_0(r^{})\right]$$
(38)
where the multiplier $`\lambda (r)`$ and the external field $`v(r)`$ are assumed fixed and we used Eqs.(31) and (27). Therefore, the best mean field, which we will denote by $`\overline{v}(r)`$ is a function of the temperature $`T`$ and a functional of the multiplier $`\lambda (r)`$ and the external field $`v(r)`$,

$$\overline{v}=\overline{v}[T,\lambda ,v],$$
(39)
defined by the equation

$$\overline{v}(r)d^3r^{}u(rr^{})\overline{n}_0(r)=0,$$
(40)
where $`\overline{n}_0`$ is a function of $`\overline{v}`$,

$$\overline{n}_0(r)\stackrel{\mathrm{d}ef}{=}\frac{e^{\beta \overline{V}(r)}}{\mathrm{\Lambda }^3}\text{and}\overline{V}(r)\stackrel{\mathrm{d}ef}{=}v(r)+\overline{v}(r)+\lambda (r).$$
(41)
To summarize, the MF approximation to the exact thermodynamic potential is given by

$$\overline{\mathrm{\Omega }}[T,\lambda ,v]=\overline{\mathrm{\Omega }}_0+\overline{U}_0d^3r\overline{v}(r)\overline{n}_0(r)$$
(42)
where

$$\overline{\mathrm{\Omega }}_0=\frac{1}{\beta }d^3r\overline{n}_0(r),$$
(43)
$$\overline{U}_0=\frac{1}{2}d^3rd^3r^{}u(rr^{})\overline{n}_0(r)\overline{n}_0(r^{}),$$
(44)
and $`\overline{v}(r)`$ is given by Eq.(40). Throughout we will use overbars to denote quantities evaluated in the MF approximation.

Notice that in replacing the exact distribution $`P`$ by the approximate $`P_0`$ we are not quite replacing the fluid described by the Hamiltonian $`H`$ by a different fluid described by $`H_0`$: we compute the expectation of the exact $`H`$ in the approximate $`P_0`$ rather than the expectation of the approximate $`H_0`$ in the approximate $`P_0`$.

### 4.3 Density in the MF approximation

Entropy and density are obtained as derivatives of the thermodynamic potential,

$$d\overline{\mathrm{\Omega }}=\overline{S}dT+d^3r\overline{n}(r)\delta \lambda (r).$$
(45)
Using Eq.(40) and (41) the density $`n(r)`$ is given by

$$\overline{n}(r)=\frac{\delta \overline{\mathrm{\Omega }}}{\delta \lambda (r)}=\overline{n}_0(r)=\frac{e^{\beta \overline{V}(r)}}{\mathrm{\Lambda }^3}.$$
(46)
This permits us to interpret Eq.(40) as a self consistency equation: the mean field $`\overline{v}(r)`$ is generated by a molecular distribution described by the density $`\overline{n}(r)`$, while the molecules distribute themselves according to the mean field.

The Legendre transform of $`\overline{\mathrm{\Omega }}[t,\lambda ]`$,

$$\overline{\mathrm{\Phi }}[T,\overline{n}]=\overline{\mathrm{\Omega }}[T,\lambda ]d^3r\lambda (r)\overline{n}(r).$$
(47)
gives the MF approximation to the density functional,

$`\overline{\mathrm{\Phi }}[T,\overline{n}]`$ $`=`$ $`{\displaystyle d^3r\overline{n}(r)\left[T\left(\mathrm{log}\mathrm{\Lambda }^3\overline{n}(r)1\right)+v(r)\right]}`$ (48)
$`+{\displaystyle \frac{1}{2}}{\displaystyle d^3rd^3r^{}u(rr^{})\overline{n}(r)\overline{n}(r^{})}.`$
The density $`\overline{n}(r)`$ in a given external potential and multiplier fields is determined by

$$\frac{\delta \overline{\mathrm{\Phi }}}{\delta \overline{n}(r)}=\lambda (r),$$
(49)
or

$$\mathrm{log}\mathrm{\Lambda }^3\overline{n}(r)=\beta \left[v(r)+\lambda (r)+d^3r^{}u(rr^{})\overline{n}(r^{})\right],$$
(50)
which is clearly equivalent to Eqs.(40) and (41).

If the only constraint on the density is a constraint on the expected total number of particles then the multiplier field $`\lambda (r)`$ is a constant which we set equal to minus the chemical potential. Therefore,

$$T\mathrm{log}\mathrm{\Lambda }^3\overline{n}(r)+v(r)+d^3r^{}u(rr^{})\overline{n}(r^{})=\mu .$$
(51)
We can either consider this equation as determining the chemical potential given the density, or alternatively, it determines the density for a given chemical potential.

### 4.4 Correlation functions – General

First we recall some standard definitions. Density correlation functions are given by derivatives of the partition function $`Z[T,\lambda ,v]`$ with respect to $`\lambda (r_i)`$. Derivatives of the thermodynamic potential $`\mathrm{\Omega }[T,\lambda ,v]`$ yield the so called ‘connected’ correlation functions; the first is the expected density itself,

$$n(r)\stackrel{\mathrm{d}ef}{=}\widehat{n}(r)=\frac{\delta \mathrm{\Omega }}{\delta \lambda (r)}.$$
(52)
The second is the density fluctuation correlation,

$`G^{(2)}(r_1,r_2)\stackrel{\mathrm{d}ef}{=}\left(\widehat{n}(r_1)n(r_1)\right)\left(\widehat{n}(r_2)n(r_2)\right)`$ (53)
$`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\delta ^2\mathrm{\Omega }}{\delta \lambda (r_1)\delta \lambda (r_2)}}={\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\delta n(r_1)}{\delta \lambda (r_2)}}`$
Using

$$n^{(2)}(r_1,r_2)=\widehat{n}^{(2)}(r_1,r_2)=\widehat{n}(r_1)\widehat{n}(r_2)n(r_1)\delta (r_1r_2),$$
(54)
$`G^{(2)}`$ can be reexpressed as

$$G^{(2)}(r_1,r_2)=n^{(2)}(r_1,r_2)+n(r_1)\delta (r_1r_2)n(r_1)n(r_2).$$
(55)
To introduce the direct (or one-particle irreducible) correlations define $`\mathrm{\Gamma }^{(2)}`$ as the inverse of $`G^{(2)}`$, $`\mathrm{\Gamma }^{(2)}=[G^{(2)}]^1`$,

$$d^3r_3G^{(2)}(r_1,r_3)\mathrm{\Gamma }^{(2)}(r_3,r_2)=\delta (r_1r_2).$$
(56)
Therefore, using Eq.(53) and the chain rule,

$$\mathrm{\Gamma }^{(2)}(r_1,r_2)=\beta \frac{\delta \lambda (r_1)}{\delta n(r_2)}=\beta \frac{\delta ^2\mathrm{\Phi }}{\delta n(r_1)\delta n(r_2)}.$$
(57)
From Eq.(48) we see that for an ideal gas we have

$$\mathrm{\Gamma }_{\mathrm{i}d}^{(2)}(r_1,r_2)=\frac{\delta (r_1r_2)}{n(r_1)}.$$
(58)
The difference between $`\mathrm{\Gamma }^{(2)}`$ and $`\mathrm{\Gamma }_{\mathrm{i}d}^{(2)}`$ is attributed to the interactions and defines the direct correlation function $`c`$,

$$\mathrm{\Gamma }^{(2)}(r_1,r_2)=\frac{\delta (r_1r_2)}{n(r_1)}c(r_1,r_2).$$
(59)
It is convenient to introduce the pair-distribution function $`g(r_1,r_2)`$

$$n^{(2)}(r_1,r_2)=n(r_1)n(r_2)g(r_1,r_2).$$
(60)
Substituting this and Eq.(55) into Eq.(56) gives

$$g(r_1,r_2)=1+c(r_1,r_2)+d^3r_3n(r_3)c(r_3,r_2)\left[g(r_1,r_3)1\right].$$
(61)
This is called the Oernstein-Zernicke equation. Notice that there is no physical content to the OZ equation; it follows exactly from the definitions of $`\mathrm{\Gamma }^{(2)}`$, $`c`$ and $`g`$. This concludes our brief summary of definitions.

At high fluid densities there is data from x-ray and neutron diffraction experiments that gives the radial distribution $`g(r)`$ directly. At low densities the experimental data is less direct. The data that will be discussed below is derived from the virial expansion to the equation of state, in particular the second virial coefficient, and from measurements of the internal energy of the fluid.

For a homogeneous, isotropic fluid of bulk density $`N/V=n_b`$ the equation of state is given in terms of $`g(r)`$ by

$$\frac{P}{n_bkT}=1\frac{n_b}{6kT}d^3rr\frac{du\left(r\right)}{dr}g\left(r\right).$$
(62)
Using eq.(13) and (60) the internal energy is

$$E=\frac{3}{2}NkT+U$$
(63)
where, setting the external potential $`v(r)`$ equal to zero, the potential energy per particle is

$$\frac{U}{N}=\frac{n_b}{2}d^3ru(r)g\left(r\right).$$
(64)
### 4.5 MF approximation to $`g(r)`$

The MF approximation to the density-fluctuation correlation is obtained using Eqs.(53) and (46),

$$\overline{G}^{(2)}(r_1,r_2)=\overline{n}(r_1)\left[\frac{\delta \overline{v}(r_1)}{\delta \lambda (r_2)}+\delta (r_1r_2)\right].$$
(65)
The functional derivative on the right is obtained from the selfconsistency equation defining the mean field $`\overline{v}`$, Eq.(40). The expression for $`\overline{\mathrm{\Gamma }}^{(2)}`$ is much simpler and can be written in closed form,

$$\overline{\mathrm{\Gamma }}^{(2)}(r_1,r_2)=\frac{\delta (r_1r_2)}{n(r_1)}+\beta u(r_1r_2).$$
(66)
Therefore the direct correlation function is

$$\overline{c}(r_1,r_2)=\beta u(r_1r_2).$$
(67)
Using this expression for $`\overline{c}`$ in the case of a homogeneous, isotropic fluid of bulk density $`n_b`$ in the OZ equation (61), gives

$$\overline{g}(r)=1\beta u(r)\beta n_bd^3r^{}u(r^{})\left[\overline{g}(rr^{})1\right],$$
(68)
where we have set $`r_2=0`$. This is the integral equation that must be solved to obtain $`\overline{g}(r)`$ for each $`T`$ and $`n_b`$.

### 4.6 An improved MF approximation to $`g(r)`$

The approximation in which each particle moves in a MF field generated by all the others is too drastic. One can obtain a somewhat better approximation – still of the MF type – by using an elegant argument due to Percus . The strategy is to write multiparticle distribution functions for a fluid in terms of the single particle density of the same fluid placed in a suitable external potential.

The idea is simple and revolves around the product rule for probabilities. The probability $`P(r)d^3r`$ that a particle is found within $`d^3r`$ at $`r`$ is obtained from Eq.(8) by marginalizing over the remaining $`N1`$ particles,

$$P(r)=\underset{i}{}\delta (rr_i)\stackrel{\mathrm{d}ef}{=}n(r|v),$$
(69)
where the average is taken with the exact distribution, Eq.(8), and for later convenience the dependence on the external potential $`v`$ is indicated explicitly. Similarly, the probability that one particle is found at $`r`$ and a second one at $`r_s`$ is given by $`P^{(2)}(r,r_s)d^3rd^3r_s`$ with

$$P(r,r_s)=\underset{i}{}\delta (rr_i)\underset{ji}{}\delta (r_sr_j)\stackrel{\mathrm{d}ef}{=}n^{(2)}(r,r_s|v).$$
(70)
The product rule gives

$$P(r,r_s)=P(r_s)P(r|r_s),$$
(71)
but

$$P(r|r_s)=N\underset{i=1}{\overset{N1}{}}\delta (rr_i)\delta (r_sr_N)\stackrel{\mathrm{d}ef}{=}n(r|v_s)$$
(72)
is exactly the one-particle distribution of a fluid in a modified external potential, $`v_s(r)=v(r)+u(rr_s)`$. Therefore,

$$n^{(2)}(r,r_s|v)=n(r_s|v)n(r|v_s),$$
(73)
and using

$$n^{(2)}(r,r_s|v)=n(r|v)n(r_s|v)g(r,r_s|v),$$
(74)
we get

$$g(r,r_s|v)=\frac{n(r|v_s)}{n(r|v)}.$$
(75)
This is exact and can be trivially generalized to multiparticle distributions.

Now we return to the MF approximation. If there is no external potential, $`v(r)=0`$ the fluid is homogeneous and isotropic with a density $`n_b`$. We single out one molecule located at $`r_s`$ and treat it as special. Since the system is homogeneous $`n(r|v=0)=n_b`$. The remaining molecules are effectively moving in the external potential $`v_s(r)=u(rr_s)`$ due to the special molecule. The density $`\overline{n}(r|v_s)`$ at location $`r`$ is given by Eq.(51),

$$T\mathrm{log}\mathrm{\Lambda }^3\overline{n}(r|v_s)=\mu u(rr_s)d^3r^{}u(rr^{})\overline{n}(r^{}|v_s).$$
(76)
On the other hand, given $`n_b`$, the chemical potential is determined from a similar equation with the external potential set to zero,

$$T\mathrm{log}\mathrm{\Lambda }^3n_b=\mu d^3r^{}u(rr^{})n_b.$$
(77)
Subtracting these two equations gives

$$T\mathrm{log}\frac{\overline{n}(r|v_s)}{n_b}=u(rr_s)n_bd^3r^{}u(rr^{})\left[\frac{\overline{n}(r^{}|v_s)}{n_b}1\right].$$
(78)
Therefore, setting $`r_s=0`$, the MF approximation to Eq.(75) is

$$\overline{g}_s(r)=\mathrm{exp}\left(\beta u(r)\beta n_bd^3r^{}u(r^{})\left[\overline{g}_s(rr^{})1\right]\right).$$
(79)
We have described two MF approximations to $`g(r|v)`$. The first requires solving Eq.(68) for $`\overline{g}(r)`$. The second is an improved MF approximation, which does include some correlation effects, is obtained by solving Eq.(79) for $`\overline{g}_s(r)`$.

## 5 An explicit example: Argon

Realistic model interatomic potentials must include a long range weak attraction and a short range strong repulsion if they are to predict the properties of dense fluids. A popular but by no means unique choice for spherically symmetric molecules is the Lennard-Jones (LJ) potential,

$$u\left(r\right)=4\epsilon \left[\left(\frac{\sigma }{r}\right)^{12}\left(\frac{\sigma }{r}\right)^6\right],$$
(80)
where the parameter $`\epsilon `$ is interpreted as the depth of the potential well and $`\sigma `$ is the molecular or ‘collision’ diameter. The reason for choosing this model over other qualitatively similar ones is a historical one. It is not that it affords a particularly accurate representation of the potential but rather that it was convenient for analytical calculations at a time when computers were not available.

The main difficulty in solving the integral equations for $`\overline{g}(r)`$ and $`\overline{g}_s(r)`$ is immediately apparent. The integrals in Eqs.(68) and (79) are not defined if the potential $`u(r)`$ diverges as $`r0`$. As a quick remedy we modify the short range behavior and set $`u(r)`$ equal to a constant below a certain cutoff distance: $`u(r)=u(r_c)`$ for $`r<r_c`$. This arbitrary procedure is justifiable whenever the physical predictions are not sensitive to the particular value of $`r_c`$. This is a serious restriction; we find that insensitivity to the cutoff holds for very low densities, it fails for even moderately low densities, and under no circumstances can the treatment be extended to liquid densities. This is not totally unexpected: behavior at high densities depends crucially on the short range correlations neglected in MF approximations.

First we discuss the iterative solution of Eqs.(68) and (79) for $`\overline{g}(r)`$ and $`\overline{g}_s(r)`$. We assume the dilute gas approximation,

$$g_d(r)=e^{\beta u(r)},$$
(81)
as the initial guess that is substituted into the integrals in the right hand side of Eqs.(68) and (79). We find that for sufficiently low densities the iterative process converges well for both the $`\overline{g}(r)`$ and $`\overline{g}_s(r)`$ and the two agree. For increasing densities the solutions develop oscillations of increasingly large amplitude until the iterations fail to converge.

In Fig.1 we show the dilute gas approximation (DGA) $`g_d`$, the mean field approximation (MFA) $`\overline{g}`$, and the improved mean field approximation (IMFA) $`\overline{g}_s`$ for Argon at $`n_b=0.00125`$Å$`^{3\text{ }}`$and $`T=150.5`$K. The values used for the LJ parameters are $`\epsilon =0.0103`$eV and $`\sigma =3.4`$Å (the choice of $`\epsilon `$ and $`\sigma `$ will be discussed below) and the cutoff is set at $`r_c=2.9`$Å (the calculation is insensitive to changes of $`r_c`$ within a range from about $`2.5`$ to $`3.1`$Å). The temperature has been chosen sufficiently low and the density sufficiently high that one is already approaching the regime where the iterative solutions fail to converge. We find that only within the restricted region $`n_b/kT0.16/`$Å<sup>3</sup>eV the solutions converge satisfactorily.

The somewhat higher peak of the IMFA can be attributed to the fact that it partially takes some correlation effects into account. At longer distances we find that MFA and IMFA are in close agreement.

Next we briefly address the two closely related questions of the determination of the LJ parameters $`\epsilon `$ and $`\sigma `$ from experimental data and the overall agreement or not of the IMFA predictions with experiment. We will use experimental data given in  for the equation of state of Argon and for its internal energy. The equation of state data are summarized through the coefficients $`B\left(T\right)`$ and $`C\left(T\right)`$ in the virial expansion

$$\frac{P}{n_bkT}=1+B\left(T\right)n_b+C\left(T\right)n_b^2+\mathrm{}$$
(82)
and data on the internal energy is given in terms of the potential energy per particle $`U/N`$. These data will be compared with the DGA and IMFA predictions derived from Eqs.(62) and (64) and with the expression for the second virial coefficient given by

$$B\left(T\right)=\frac{2\pi }{3kT}_0^{\mathrm{}}𝑑rr^3e^{\beta u\left(r\right)}\frac{du\left(r\right)}{dr}.$$
(83)
The generally accepted values of the LJ parameters are those given in  by fitting $`B(T)`$ at $`T=273`$K,

$$\epsilon _1=0.01034\text{eV}\text{and }\sigma _1=3.405\text{Å\hspace{0.17em}.}$$
(84)
In Fig.2 the DGA and IMFA equations of state for the parameters $`\epsilon _1`$ and $`\sigma _1`$ are compared with experimental data. The IMFA curve is in sharp disagreement with the data. Also shown are the DGA and IMFA predictions for a second set of parameters,

$$\epsilon _2=0.01032\text{eV}\text{and}\sigma _2=3.151\text{Å}.$$
(85)
The DGA curve fits the experimental curve for low density equally well for either set of parameters, while the IMFA curve for the second set is clearly the best match with the data. The difference between the $`\epsilon `$ parameters is not significant.

On the other hand, as shown in Fig.3 when the DGA and IMFA predictions for the potential energy per particle are compared with experiment it is not the second set of LJ parameters but the first that yields the better match.

To study this perhaps surprising result from a different perspective we plot in Fig.4 the second virial coefficient $`B`$ at $`T=273`$K for both the DGA and the IMFA as a function of the collision diameter $`\sigma `$ for $`\epsilon =0.01034`$eV. Our calculation shows that there are two solutions for $`\sigma `$ that match the experimental value, and these two values are precisely the $`\sigma _1`$ and $`\sigma _2`$ in Eqs.(84) and (85).

While it is clear that the IMFA is a definite improvement over the DGA a number of interesting questions remain unanswered. First, there is the question of sensitivity of the MF approximations to the short range part of the potential and the issue of convergence of the iterative solutions. Both deserve to be studied further. Then, there is the fact that different sets of LJ parameters are needed to fit different data. Indeed, we found that different values of $`\epsilon `$ and $`\sigma `$ are needed to fit data at different temperatures. This can be interpreted as due to intrinsic inaccuracies in the approximations (DGA, MFA or IMFA) employed or, perhaps more likely, to the fact that the Lennard-Jones potential is not a good representation of the true interatomic potential. Finally there are the inevitable uncertainties in the experimental data which can be considerable .

## 6 Conclusions

The method of maximum entropy is traditionally used to assign appropriate probability distributions. Here we have explored a different use of the method of maximum entropy, as a technique to generate approximations to probability distributions that may be intractably difficult to calculate with. We find that the resulting variational technique includes the well-known Bogoliuvob variational principle as a special case. The extension to other generalized canonical ensembles is straightforward and we explicitly show its application to the density functional formalism and the derivation of mean field approximations for classical fluids. A side result of some interest is a simple proof that the density functional formalism is itself an application of the method of maximum entropy.

Numerical calculations for Argon gas were compared with experimental data. Just as with other mean field approximations, the particular MF versions studied here represent improvements over the dilute gas approximation but remain restricted to very small densities where particle correlations are small. However, as with so many other variational principles the power of the method hinges on the right choice of family of trial distributions. Indeed, following the ME-based variational approach still has the advantage that the same method can easily be used to generate other approximations without these shortcomings. If one’s interest lies in the physics of denser fluids the next natural step consists in choosing a family of trial distributions that provides a better representation of correlation effects.