text:
# Coherence properties of a continuous atom laser

## 1 Evaporative cooling of an atomic beam

### 1.1 A continuous injection of atoms in a magnetic guide

In the experimental scenario considered at the École normale supérieure the continuous source of atoms is provided by standard laser cooling and trapping techniques taking place inside a cell. Atoms in the cell are captured, cooled and trapped in the $`xy`$ plane using a two-dimensional magneto-optical trap, that is with laser beams in the $`xy`$ plane and a dipolar magnetic field $`\stackrel{}{B}x\stackrel{}{e}_xy\stackrel{}{e}_y`$, where $`\stackrel{}{e}_{x,y}`$ are unit vectors along $`x,y`$ axes. The atomic motion along $`z`$ is controlled with a standard moving molasses technique: the two counter-propagating laser beams along $`z`$ have different frequencies so that the atoms are cooled around a non-zero mean velocity $`\overline{v}_0`$.

The already cold atomic beam emerging from this set-up is sent towards a magnetic guide of axis $`z`$. The magnetic guide is produced by the superposition of a uniform magnetic field along $`z`$ and a dipolar magnetic field $`x\stackrel{}{e}_xy\stackrel{}{e}_y`$. This provides a transverse confinement of the atoms being in the right Zeeman sublevels, the atoms in the wrong Zeeman sublevels being not trapped or even expelled. The magnetic guide provides a trapping potential in the $`xy`$ plane with a harmonic bottom that we write as

$$U(x,y)=\frac{1}{2}m\omega _{}^2(x^2+y^2)$$
(1)
where $`m`$ is the atomic mass and $`\omega _{}`$ the transverse oscillation frequency of the atoms.

The following parameters of the injected atoms are expected to be realistic for 87Rb atoms . The initial velocity dispersion of the atoms is $`\mathrm{\Delta }v_0=20`$ cm/s corresponding to an initial temperature of 400 $`\mu `$K; such a high temperature is obtained after spatial compression of the cloud usually performed to increase the collision rate in preparation of evaporative cooling. We take an injection velocity $`\overline{v}_0=3\mathrm{\Delta }v_0`$ larger than the velocity spread so that the incoming atoms form a beam. Assuming an injected flux of atoms of $`3\times 10^9`$ atoms/s and an oscillation frequency $`\omega _{}=2\pi \times 1`$ kHz we find an initial on-axis thermal density of $`8\times 10^{11}`$ atoms/cm3. The initial phase space density is $`7\times 10^71`$. The initial collision rate of atoms is related to the on-axis density $`n_0`$ and to the collisional cross section $`\sigma `$ by $`\gamma _{\mathrm{coll}}^{(0)}=(2/\sqrt{\pi })n_0\sigma \mathrm{\Delta }v_0`$. For the $`s`$-wave collisional cross-section of rubidium ($`\sigma =7.6\times 10^{16}`$ m2) we get $`\gamma _{\mathrm{coll}}^{(0)}`$ 100 s-1 which is much smaller than $`\omega _{}`$:

$$\frac{\gamma _{\mathrm{coll}}^{(0)}}{\omega _{}}=0.021.$$
(2)
The atoms have therefore the time to perform a full transverse oscillation in the trapping potential before experiencing a collision.

### 1.2 Modeling of evaporative cooling

We assume that the atoms are subject to evaporative cooling in the magnetic guide, e.g. by application of a $`z`$-dependent radio-frequency flipping the atoms to untrapped or expelled Zeeman sublevels when they are too far from $`z`$ axis.

The dynamics of evaporative cooling can be described by the classical Boltzmann equation on the phase space atomic density $`f(\stackrel{}{r},\stackrel{}{p})`$ as long as the phase space density remains small. The collision terms in the Boltzmann equation are simplified by the assumption (justified for rubidium) that atomic interactions take place in the $`s`$-wave only and have a constant (momentum independent) total cross-section $`\sigma =8\pi a_{3d}^2`$ ($`a_{3d}`$ is the scattering length). The effect of evaporative cooling is modeled by setting $`f(\stackrel{}{r},\stackrel{}{p})`$ to zero in the domain $`x^2+y^2>\mathrm{\Lambda }^2(z)`$ for two-dimensional evaporation, or in the domain $`x^2>\mathrm{\Lambda }^2(z)`$ for one-dimensional evaporation, where $`\mathrm{\Lambda }(z)`$ is an adjustable cut in position space.

We have first looked for an approximate analytical solution of the Boltzmann equation, adapting to our geometry the truncated Gaussian ansatz put forward in . As argued in  it is then more rigorous to restrict to a one-dimensional evaporation scheme. The ansatz then takes the following form:

$$f(\stackrel{}{r},\stackrel{}{p})=f_0(z)e^{(ϵ_x+ϵ_y)/(k_BT(z))}e^{(p_z\overline{p}(z))^2/(2mk_BT(z))}Y(ϵ_{\mathrm{evap}}(z)ϵ_x)$$
(3)
where $`ϵ_x,ϵ_y`$ are the sum of kinetic energy and harmonic trapping potential energy along $`x`$ and $`y`$ respectively and $`Y`$ is the Heaviside function. The ansatz assumes a local thermal equilibrium with temperature $`T(z)`$. The temperature depends on $`z`$ only, not on $`x`$ and $`y`$ as, due to Eq.(2), the mean free path of the particles $`\mathrm{\Delta }v/\gamma _{\mathrm{coll}}`$ is much larger than the spatial transverse extension of the gas $`(k_BT/(m\omega _{}^2))^{1/2}\mathrm{\Delta }v/\omega _{}`$: transversally the gas is in the so-called collisionless regime. For the same reason the truncation of $`f`$ in position space is replaced by a truncation in energy space, with

$$ϵ_{\mathrm{evap}}(z)=\frac{1}{2}m\omega _{}^2\mathrm{\Lambda }^2(z).$$
(4)
Knowing the energy of the particle along $`x`$ allows to calculate the maximal excursion of the trajectory along $`x`$ (as atoms have in general the time to perform a full harmonic oscillation before experiencing a collision); if this maximal excursion exceeds $`\mathrm{\Lambda }(z)`$ the particle is evaporated.

There are a priori three unknown functions of $`z`$ in the ansatz Eq.(3): (i) the normalization factor $`f_0(z)`$ or equivalently the linear density $`\rho _{\mathrm{lin}}(z)`$, (ii) the mean momentum $`\overline{p}(z)`$ of the gas along $`z`$, and (iii) the temperature of the gas $`T(z)`$.

By multiplying Boltzmann’s equation by (i) unity, (ii) the momentum $`p_z`$, and (iii) the kinetic energy along $`z`$, $`p_z^2/(2m)`$, and by integrating over $`x,y,p_x,p_y,p_z`$ one gets three hydrodynamic type equations for $`\rho _{\mathrm{lin}}(z)`$, $`\overline{p}(z)`$ and $`k_BT(z)`$. These three equations contain the usual equations expressing conservation of probability, of momentum and of energy, plus extra terms describing the loss of particles, the change of momentum and energy due to the evaporation.

We have solved numerically the hydrodynamic type equations in steady state, assuming that the $`z`$-dependence of the parameter $`\mathrm{\Lambda }(z)`$ is adjusted to maintain a $`z`$-independent ratio $`\eta =ϵ_{\mathrm{evap}}(z)/(k_BT(z))`$. For the specific set of parameters of §1.1 we have to gain seven orders of magnitude on the phase space density to reach quantum degeneracy. The smallest spatial length of evaporation required is obtained for $`\eta 5`$ as shown in figure 1; it is equal to 7600 $`d_0`$ where $`d_0=\sqrt{\pi }/(2n_0\sigma )`$ is the mean free path at the entrance of the magnetic guide, that is $`11`$ meters for the considered parameters. After evaporative cooling along these 11 meters the flux of particles has been reduced by a factor 90 and the temperature has been decreased by a factor 4000.

We have also performed a numerical simulation of the full Boltzmann equation by a Monte Carlo algorithm using macro-atoms . The resulting numerical calculations take several days on a workstation. We have first simulated the one-dimensional evaporation scheme. In this way we have confirmed the accuracy of the predictions based on the ansatz Eq.(3). We have also performed simulations for the full two-dimensional evaporation. We have found typically that the length required to reach quantum degeneracy is reduced by a factor three as compared to one-dimensional evaporation, with the same loss of two orders of magnitude on the flux. For the specific set of parameters of §1.1 the required evaporation length reaches the experimentally reasonable value of 4 meters.

## 2 Coherence properties in the ideal Bose gas model

We now assume that evaporative cooling has allowed to reach the quantum degenerate regime at a certain distance from the entrance of the tube. We have not performed any kinetic study of the approach of quantum degeneracy, and we will here simply assume that the gas is at thermal equilibrium in the frame moving at the mean velocity of the gas. Such an assumption is reasonable if the temperature remains significantly larger than the quantum $`\mathrm{}\omega _{}`$ of transverse oscillation of the atoms of the guide. If the temperature was much smaller than $`\mathrm{}\omega _{}/k_B`$ the transverse degrees of freedom of the gas would be frozen in the ground state of the transverse harmonic oscillator; the gas would become a free one-dimensional Bose gas along $`z`$, that could not thermalize as colliding identical particles in one dimension simply exchange their momenta.

In this section we consider the model of the ideal Bose gas. The effects of atomic interactions are discussed in the next section.

### 2.1 Transverse Bose-Einstein condensation

Let us enclose the ideal Bose gas in a fictitious box of size $`L`$ along $`z`$, with periodic boundary conditions. Transversally the gas is confined by the harmonic potential Eq.(1). The one-particle eigenstates of the system are then labeled by three integers, the non-negative integers $`l_x,l_y`$ labeling the eigenstates of the harmonic oscillator along $`x`$ and $`y`$, and the integer $`l_z`$ labeling the momentum along $`z`$:

$$\mathrm{}k_z=\frac{2\pi }{L}l_z.$$
(5)
In the grand canonical ensemble the mean occupation number of the single particle state $`\stackrel{}{l}=(l_x,l_y,l_z)`$ is given by

$$n(\stackrel{}{l},\mu )=\left\{\mathrm{exp}\left[\beta \left(\frac{\mathrm{}^2k_z^2}{2m}+(l_x+l_y)\mathrm{}\omega _{}\mu \right)\right]1\right\}^1.$$
(6)
For convenience we have included the transverse zero-point motion energy in the chemical potential, so that $`\mu `$ varies from $`\mathrm{}`$ to 0.

It turns out that in our trapping geometry there is no Bose-Einstein condensation in the thermodynamical limit defined as $`L,N+\mathrm{}`$ with a fixed linear density $`N/L`$, $`N`$ being the mean number of particles.

Let us consider indeed a fixed value of $`L`$ and let us define the maximal mean number of particles $`N_{\mathrm{max}}^{}`$ that can be put in all states but the ground state of the trap at a fixed temperature. A Bose-Einstein condensate forms in the ground state of the trap when $`N`$ exceeds $`N_{\mathrm{max}}^{}`$. As each $`n(\stackrel{}{l},\mu )`$ (for $`\stackrel{}{l}\stackrel{}{0}`$) reaches its maximal accessible value for $`\mu =0`$, $`N_{\mathrm{max}}^{}`$ is given by

$$N_{\mathrm{max}}^{}=\underset{\stackrel{}{l}\stackrel{}{0}}{}n(\stackrel{}{l},\mu =0).$$
(7)
A lower bound on $`N_{\mathrm{max}}^{}`$ is obtained by restricting the sum to the states transversally in the ground state of the harmonic oscillator:

$$N_{\mathrm{max}}^{}\underset{l_z0}{}\left[\mathrm{exp}\left(\beta \frac{h^2l_z^2}{2mL^2}\right)1\right]^1.$$
(8)
When $`L`$ is large enough so that $`k_BTh^2/(2mL^2)`$ the exponential in the denominator can be expanded to first order, leading to

$$N_{\mathrm{max}}^{}\frac{\pi }{3}\left(\frac{L}{\lambda }\right)^2$$
(9)
where we have introduced the thermal de Broglie wavelength

$$\lambda =\frac{h}{(2\pi mk_BT)^{1/2}}.$$
(10)
One then realizes that $`N_{\mathrm{max}}^{}`$ grows faster than $`N`$ in the thermodynamical limit $`L+\mathrm{}`$, with $`N/L`$ fixed. This is connected to the known fact that there is no Bose-Einstein condensation in a one-dimensional homogeneous Bose gas in the thermodynamical limit.

On the other hand, transverse Bose-Einstein condensation  is taking place in our system. Let us calculate indeed the maximal number of atoms that can be put for a fixed temperature in the transversally excited single particle states:

$$N_{\mathrm{max}}^{}=\underset{(l_x,l_y)(0,0)}{}\underset{l_z}{}n(\stackrel{}{l},\mu =0).$$
(11)
Replacing in the large $`L`$ limit the sum over $`l_z`$ by an integral and expanding $`1/(\mathrm{exp}(x)1)`$ in powers of $`\mathrm{exp}(x)`$ we obtain:

$$N_{\mathrm{max}}^{}\frac{L}{\lambda }\underset{s1}{}s^{1/2}\left\{\left[1\mathrm{exp}\left(s\beta \mathrm{}\omega _{}\right)\right]^21\right\}\frac{L}{\lambda }\left(\frac{k_BT}{\mathrm{}\omega _{}}\right)^2\zeta (5/2).$$
(12)
The last equality is correct in the limit $`k_BT\mathrm{}\omega _{}`$. The function $`\zeta `$ is the Zeta function of Riemann, and $`\zeta (5/2)1.341`$. If the linear density of the gas $`\rho _{\mathrm{lin}}=N/L`$ exceeds the critical value

$$\rho _{\mathrm{lin}}^{(c)}=N_{\mathrm{max}}^{}/L\frac{1}{\lambda }\left(\frac{k_BT}{\mathrm{}\omega _{}}\right)^2\zeta (5/2)$$
(13)
the excess of density will accumulate in the transverse ground state of the trap. At linear densities much higher than $`\rho _{\mathrm{lin}}^{(c)}`$ the gas becomes almost monomode transversally: this is an interesting feature for atom optics applications as we have achieved in this way a Heisenberg limited transverse focalization of the beam.

We can also consider the maximal on-axis density of atoms in the excited transverse states. In this case it is more convenient to label the eigenstates of the two-dimensional harmonic oscillator by the angular momentum $`M`$ along $`z`$ and the radial quantum number $`l_\mathrm{r}`$. Then the wavefunctions on $`z`$ axis, that is in $`x=y=0`$, have a squared modulus equal to zero if $`M0`$ and equal to $`m\omega _{}/(\pi \mathrm{})`$ for $`M=0`$. We recall that the states with $`M=0`$ have an energy above the zero-point energy given by $`2l_\mathrm{r}\mathrm{}\omega _{}`$. We then obtain the maximal on-axis density of atoms in excited transverse states:

$$\rho _{\mathrm{axis}}^{(c)}=\frac{m\omega _{}}{\pi \mathrm{}L}\underset{l_\mathrm{r}1}{}\underset{l_z}{}\left\{\mathrm{exp}\left[\beta \left(\frac{h^2l_z^2}{2mL^2}+2l_\mathrm{r}\mathrm{}\omega _{}\right)1\right]\right\}^1.$$
(14)
With the same algebraic transformations as for Eq.(12) we obtain

$$\rho _{\mathrm{axis}}^{(c)}\frac{m\omega _{}}{\pi \mathrm{}\lambda }\underset{s1}{}s^{1/2}\left[\frac{1}{1\mathrm{exp}(2s\beta \mathrm{}\omega _{})}1\right]\frac{\zeta (3/2)}{\lambda ^3}$$
(15)
where $`\zeta (3/2)2.612`$ and we have used $`k_BT\mathrm{}\omega _{}`$. Transverse condensation therefore takes place when the usual Einstein’s condition $`\rho \lambda ^3=\zeta (3/2)`$ is satisfied on the axis of the trap!

### 2.2 In the quantum degenerate regime

We assume now that the gas is in the strongly degenerate regime, with a linear density $`\rho _{\mathrm{lin}}`$ larger than the critical value Eq.(13). The linear density of atoms in the excited transverse states has reached its saturated value $`\rho _{\mathrm{lin}}^{(c)}`$. This implies that $`|\mu |`$ is much smaller than $`\mathrm{}\omega _{}`$ ($`\mu `$ can be replaced by zero for the transversally excited states); as $`k_BT\mathrm{}\omega _{}`$ one has also

$$|\mu |k_BT$$
(16)
so that the occupation number of the absolute trap ground state $`\stackrel{}{l}=\stackrel{}{0}`$, though not of order $`N`$, is much larger than unity:

$$n(\stackrel{}{l}=0)=\frac{1}{\mathrm{exp}(\beta \mu )1}1.$$
(17)
Let us calculate the linear density of atoms in the transverse ground state of the trap:

$`\rho _{\mathrm{lin}}^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{l_z}{}}n(0,0,l_z;\mu )`$ (18)
$`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{l_z}{}}\left\{\mathrm{exp}[\beta (\mathrm{}^2k_z^2/(2m)\mu )]1\right\}^1.`$ (19)
As $`\beta |\mu |1`$ and the sum is one-dimensional the main contribution to the sum comes from states with kinetic energies on the order of $`|\mu |`$. Expanding the exponential in the Bose formula to first order we obtain a Lorentzian approximation for the occupation number as function of momentum:

$$n(0,0,l_z;\mu )\frac{k_BT}{\mathrm{}^2k_z^2/(2m)+|\mu |}.$$
(20)
Replacing finally the sum by an integral we obtain

$$\rho _{\mathrm{lin}}^{(0)}=\frac{1}{\lambda }\left(\frac{\pi k_BT}{|\mu |}\right)^{1/2}.$$
(21)
This allows to express the chemical potential as function of density, when combined with the relation

$$\rho _{\mathrm{lin}}=\rho _{\mathrm{lin}}^{(c)}+\rho _{\mathrm{lin}}^{(0)}.$$
(22)
Note that such a calculation would fail for a two-dimensional or three-dimensional free Bose gas, a Lorentzian momentum distribution being not normalizable in this case.

To characterize the coherence properties of the gas we use correlation functions for the atomic field operator $`\widehat{\mathrm{\Psi }}(x,y,z)`$ in direct analogy with the correlation functions considered in optics for the photonic field .

We define the first order correlation function as

$$g_1(z)=\widehat{\mathrm{\Psi }}^{}(0,0,z)\widehat{\mathrm{\Psi }}(0,0,0)$$
(23)
where the expectation value is taken in thermal equilibrium; this function is sensitive to the coherence of the atomic field between two points on the axis of the guide separated by a distance $`|z|`$. It can be written as the sum of the contributions of the transversally excited states and of the states in the transverse ground state. The two contributions behave in a very different way in the degenerate limit, that is in the limit $`\mu 0`$. This is illustrated in figure 2 for a moderately degenerate regime.

By calculations similar to the ones leading to Eqs.(14,15) we find for the contribution of the transverse excited states to $`g_1`$ in the thermodynamical limit:

$$g_1^{}(z)=\frac{m\omega _{}}{\pi \mathrm{}\lambda }\underset{s1}{}\frac{e^{\beta \mu s}}{s^{1/2}}\left[\frac{1}{1\mathrm{exp}(2s\beta \mathrm{}\omega _{})}1\right]\mathrm{exp}\left(\pi \frac{z^2}{s\lambda ^2}\right).$$
(24)
Using furthermore the fact that $`|\mu |<\mathrm{}\omega _{}k_BT`$ we set $`\mu =0`$ and we expand the exponential function between square brackets to first order in $`\beta \mathrm{}\omega _{}`$:

$$g_1^{}(z)\frac{1}{\lambda ^3}\underset{s1}{}\frac{1}{s^{3/2}}\mathrm{exp}\left(\pi \frac{z^2}{s\lambda ^2}\right).$$
(25)
The maximal value of $`g_1^{}`$ is obtained in $`z=0`$ and is equal to $`\rho _{\mathrm{axis}}^{(c)}`$. The half-width of $`g_1^{}`$ is on the order of $`0.75\lambda `$.

The contribution to $`g_1`$ of the states in the transverse ground state is given by

$$g_1^{(0)}(z)=\frac{m\omega _{}}{\pi \mathrm{}L}\underset{l_z}{}n(0,0,l_z;\mu )e^{ik_zz}.$$
(26)
In the thermodynamical limit, we replace the sum by an integral. We use the Lorentzian approximation for the occupation numbers Eq.(20) and calculate its Fourier transform. This leads to a correlation function being an exponential function of $`|z|`$:

$$g_1^{(0)}(z)\frac{m\omega _{}}{\pi \mathrm{}}\rho _{\mathrm{lin}}^{(0)}e^{z/l_c}$$
(27)
with a coherence length

$$l_c=\frac{\mathrm{}}{(2m|\mu |)^{1/2}}=\frac{\lambda ^2}{2\pi }\rho _{\mathrm{lin}}^{(0)}.$$
(28)
The exponential decay of $`g_1`$ found here is to be contrasted with $`g_1`$ going to a finite value in the case of a three dimensional homogeneous Bose-Einstein condensate . The coherence length $`l_c`$ can be larger than the thermal de Broglie wavelength $`\lambda `$. It can be expressed in terms of the on-axis density using the relation $`g_1^{(0)}(z=0)=\rho _{\mathrm{axis}}\rho _{\mathrm{axis}}^{(c)}\rho _{\mathrm{axis}}^{(0)}`$:

$$l_c=\frac{\lambda ^2\mathrm{}}{2m\omega _{}}\rho _{\mathrm{axis}}^{(0)}.$$
(29)
A related issue is the expression connecting the on-axis and the linear density of the atoms in the transverse Bose-Einstein condensate:

$$\frac{\rho _{\mathrm{axis}}^{(0)}}{\rho _{\mathrm{axis}}^{(c)}}\frac{2\zeta (5/2)}{\zeta (3/2)}\frac{k_BT}{\mathrm{}\omega _{}}\frac{\rho _{\mathrm{lin}}^{(0)}}{\rho _{\mathrm{lin}}^{(c)}}.$$
(30)
This expression holds in the regime $`k_BT\mathrm{}\omega _{}`$ so that a modest value of $`\rho _{\mathrm{lin}}`$ above the critical value $`\rho _{\mathrm{lin}}^{(c)}`$ may actually correspond to a strongly degenerate regime $`\rho _{\mathrm{axis}}\lambda ^32.612`$.

We define the second order correlation function of the atomic field as

$$g_2(z)=\widehat{\mathrm{\Psi }}^{}(0,0,z)\widehat{\mathrm{\Psi }}^{}(0,0,0)\widehat{\mathrm{\Psi }}(0,0,0)\widehat{\mathrm{\Psi }}(0,0,z).$$
(31)
From a field point of view $`g_2`$ is the correlation function of the intensity of the field; it is a measure of the intensity fluctuations of the field. From a corpuscular point of view $`g_2(z)`$ is proportional to the pair distribution function of the atoms in the gas, that is the probability density to find a pair of atoms separated by a distance $`|z|`$ in the gas.

For the ideal Bose gas in the grand canonical ensemble the use of Wick’s theorem readily allows one to express $`g_2`$ in terms of $`g_1`$ by performing all the possible binary contractions of the field operators:

$$g_2(z)=g_1^2(0)+g_1^2(z).$$
(32)
We find the unpleasant result that the gas is subject to large intensity fluctuations over a length scale on the order of the coherence length $`l_c`$. In particular

$$\frac{g_2(0)}{g_1^2(0)}=2.$$
(33)
This value of two is typical of a bosonic bunching effect that manifests itself in an Hanbury-Brown and Twiss experiment in optics with thermal sources. It has to be contrasted with $`g_2g_1^2`$ obtained with laser light for photons or with almost pure Bose-Einstein atomic condensates . The output of our magnetic guide cannot be termed an ‘atom-laser’ if $`g_2(z)`$ significantly differs from $`g_1^2(0)`$ over the coherence length of the field. Fortunately we shall see in the next sections that atomic interactions can improve the situation.

Remark:

A careful reader may argue that it is dangerous to use the grand canonical ensemble for the ideal Bose gas to calculate fluctuations of the number of particles, and therefore of the field intensity. This fear is justified when a condensate is formed, a well known problem in three-dimensions: the number of particles in the condensate has then unphysically large fluctuations. This problem does not take place here in the thermodynamical limit where no Bose-Einstein condensate is formed. More precisely one can deduce from $`g_2`$ the standard deviation of the total number of particles:

$$\frac{\mathrm{\Delta }N}{N}=\left(\frac{l_c}{L}\right)^{1/2}.$$
(34)
The relative fluctuations in the number of particles become small as soon as the length $`L`$ of the gas becomes much larger than the coherence length $`l_c`$.

To illustrate the strong intensity fluctuations of the ideal Bose gas in a dramatic way we introduce the Sudarshan-Glauber P representation of the many-body density operator :

$$\widehat{\sigma }=𝒟\mathrm{\Psi }|\text{coh}:\mathrm{\Psi }\text{coh}:\mathrm{\Psi }|P(\{\mathrm{\Psi }\},\{\mathrm{\Psi }^{}\}).$$
(35)
This expression is a functional integral $`𝒟\mathrm{\Psi }`$ over the real part and the imaginary part of the c-number field $`\mathrm{\Psi }(\stackrel{}{r})`$. It involves the coherent or Glauber state of the atomic field associated to $`\mathrm{\Psi }`$:

$$|\text{coh}:\mathrm{\Psi }\mathrm{exp}\left[\frac{1}{2}d^3\stackrel{}{r}|\mathrm{\Psi }(\stackrel{}{r})|^2\right]\mathrm{exp}\left[d^3\stackrel{}{r}\mathrm{\Psi }(\stackrel{}{r})\widehat{\mathrm{\Psi }}^{}(\stackrel{}{r})\right]|\text{vacuum}.$$
(36)
Here the many-body density operator $`\widehat{\sigma }`$ is the grand canonical thermal density operator

$$\widehat{\sigma }\mathrm{exp}[\beta (\widehat{H}\mu \widehat{N})]$$
(37)
where $`\widehat{H}`$ is the Hamiltonian of the gas containing the kinetic energy and the trapping potential energy, and $`\widehat{N}`$ is the operator total number of particles. The Glauber distribution function $`P`$ can then be calculated exactly . One expands the field $`\mathrm{\Psi }`$ on the eigenmodes of the trap:

$$\mathrm{\Psi }(\stackrel{}{r})\underset{\stackrel{}{l}}{}\alpha _\stackrel{}{l}\varphi _{l_x}(x)\varphi _{l_y}(y)\frac{1}{L^{1/2}}e^{ik_zz}$$
(38)
where $`\varphi _n`$, $`n=0,1,\mathrm{}`$ are the normalized eigenfunctions of the 1D harmonic oscillator of frequency $`\omega _{}`$ and where the plane waves along $`z`$ have a wavevector given by Eq.(5). Then the Glauber distribution function is simply a product over all modes of Gaussian distributions with squared widths given by the occupation number of the modes:

$$P(\{\mathrm{\Psi }\},\{\mathrm{\Psi }^{}\})\mathrm{exp}\left[\underset{\stackrel{}{l}}{}\frac{|\alpha _\stackrel{}{l}|^2}{n(\stackrel{}{l})}\right].$$
(39)
As this distribution $`P`$ is positive the thermal equilibrium $`\widehat{\sigma }`$ can be viewed exactly as a statistical mixture of coherent states. One can then imagine that a given experimental realization of the Bose gas is a coherent state characterized by a field $`\mathrm{\Psi }`$. This field $`\mathrm{\Psi }`$ is stochastic as it varies in an unpredictable way from one experimental realization to the other. As the coherent state $`|\text{coh}:\mathrm{\Psi }`$ is an eigenstate of $`\widehat{\mathrm{\Psi }}(\stackrel{}{r})`$ with the eigenvalue $`\mathrm{\Psi }(\stackrel{}{r})`$ one can check that the correlation functions $`g_1`$ and $`g_2`$ are equal to the following averages over all possible realizations of the c-number field:

$`g_1(z)`$ $`=`$ $`\mathrm{\Psi }^{}(z)\mathrm{\Psi }(0)_{\mathrm{stoch}}`$ (40)
$`g_2(z)`$ $`=`$ $`|\mathrm{\Psi }(z)|^2|\mathrm{\Psi }(0)|^2_{\mathrm{stoch}}.`$ (41)
We have plotted in figure 3 the intensity of the field as function of position for two numerically generated realizations of $`\mathrm{\Psi }`$. Figure 3a corresponds to a non-degenerate situation; the only spatial scale for the intensity fluctuations is the thermal de Broglie wavelength $`\lambda `$. Figure 3b corresponds to a strongly degenerate regime; there are clearly two spatial scales for the intensity fluctuations, one on the order of $`\lambda `$ coming from the transversally non-condensed fraction and the other one on the order of $`l_c`$ due to the gas in the transverse ground state of the trap. The large intensity fluctuations at the scale of $`l_c`$ manifest themselves as “droplets” in the atomic density.

## 3 The interacting case: a one-dimensional classical field model

### 3.1 Model Hamiltonian for the transversally Bose condensed gas

We have seen on the ideal Bose gas model that transverse Bose-Einstein condensation takes place in the magnetic guide at sufficient high density. Although we have not performed any detailed analysis we expect the same phenomenon to occur for the interacting Bose gas in the weakly interacting regime.

At the presently considered low temperatures, interaction between the atoms takes place in $`s`$-wave mainly, and the relative wavevector of two colliding atoms is much smaller than the inverse of the scattering length $`a_{3d}`$ of the interaction potential. The interaction potential is then commonly replaced by an effective low energy interaction potential $`g_{3d}\delta (\stackrel{}{r}_1\stackrel{}{r}_2)`$  with a coupling constant

$$g_{3d}=\frac{4\pi \mathrm{}^2}{m}a_{3d}.$$
(42)
The resulting Hamiltonian is expressed in terms of the field operator $`\widehat{\mathrm{\Psi }}(x,y,z)`$ as

$`\widehat{H}_{3d}={\displaystyle 𝑑x𝑑y_0^L𝑑z\left[\frac{\mathrm{}^2}{2m}|\stackrel{}{\text{grad}}\widehat{\mathrm{\Psi }}|^2+\frac{1}{2}g_{3d}\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}\widehat{\mathrm{\Psi }}+(U(x,y)\mu )\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}\right]}`$ (43)
with the trapping potential given in Eq.(1) and with periodic boundary conditions along $`z`$. We have included for convenience the term $`\mu \widehat{N}`$ where $`\mu `$ is the chemical potential so that the thermal equilibrium density operator is simply $`\mathrm{exp}[\beta \widehat{H}_{3d}]`$ in the grand canonical ensemble.

In our magnetic guide geometry the situation is particularly simple if the typical interaction energy per particle $`\rho _{\mathrm{axis}}g_{3d}`$ is smaller than the quantum of transverse oscillation $`\mathrm{}\omega _{}`$. Transverse Bose-Einstein condensation will then take place in a transverse wavefunction $`\varphi _{}(x,y)`$ close to the ground state of the harmonic oscillator. If we wish to describe only the one-dimensional Bose gas of atoms in this transverse condensate, assuming that the remaining atoms have a much smaller density, we can neglect the contribution of the transverse modes to the field operator by setting

$$\widehat{\mathrm{\Psi }}(x,y,z)\varphi _{}(x,y)\widehat{\psi }(z).$$
(44)
This will eliminate the contribution of the atoms in the transverse excited states of the trap; we keep in mind however that these atoms are essential to ensure thermalization, with a temperature $`k_BT>\mathrm{}\omega _{}`$. The reduced field operator $`\widehat{\psi }`$ obeys the usual bosonic commutation relations of a one-dimensional bosonic field. Inserting the approximate expression for $`\widehat{\mathrm{\Psi }}`$ in the Hamiltonian results in the following model Hamiltonian for the transversally condensed Bose gas:

$$\widehat{H}=_0^L𝑑z\left[\frac{\mathrm{}^2}{2m}|_z\widehat{\psi }|^2+\frac{1}{2}g\widehat{\psi }^{}\widehat{\psi }^{}\widehat{\psi }\widehat{\psi }\mu \widehat{\psi }^{}\widehat{\psi }\right].$$
(45)
It corresponds to a one-dimensional Bose gas with a contact interaction potential $`g\delta (z_1z_2)`$ between particles with an effective coupling constant

$$g=g_{3d}𝑑x𝑑y|\varphi _{}(x,y)|^4.$$
(46)
### 3.2 Classical field approximation

The model Hamiltonian Eq.(45) leads to an exactly solvable $`N`$-body problem: exact eigenenergies and eigenfunctions are known . It is however not easy to extract information from the exact solution, even at zero temperature.

We use here a simpler approach, valid in a sufficiently high temperature regime. The idea is to write the thermal equilibrium density operator (up to a normalization factor) as the result of a fictitious time evolution:

$$\frac{d}{d\tau }\widehat{\sigma }=\frac{1}{2}\left(\widehat{H}\widehat{\sigma }+\widehat{\sigma }\widehat{H}\right)$$
(47)
with the initial condition $`\widehat{\sigma }(\tau =0)`$ equal to the identity operator. This evolution corresponds to the so-called imaginary time evolution. It leads to a density operator at ‘time’ $`\tau `$ given by:

$$\widehat{\sigma }=e^{\tau \widehat{H}}$$
(48)
which is (up to a normalization factor) the thermal equilibrium density operator at temperature $`k_BT=1/\tau `$. High temperatures correspond to low values of $`\tau `$ that is to short ‘time’ evolution.

To take advantage of the fact that the ‘time’ evolution is short it is better to rewrite Eq.(47) using some of the representations of the density operator introduced in quantum optics. We use here the Glauber P distribution already introduced in section 2. It is defined as in Eq.(35) with the difference that $`\psi `$ is now a function of the coordinate $`z`$ only. So we apply the Glauber transform to both sides of Eq.(47). The transform of the right-hand is performed with the following rules:

$`\text{Glaub}[\widehat{\psi }(z)\widehat{\nu }]`$ $`=`$ $`\psi (z)\text{Glaub}[\widehat{\nu }]`$ (49)
$`\text{Glaub}[\widehat{\psi }^{}(z)\widehat{\nu }]`$ $`=`$ $`\left(\psi ^{}(z)_{\psi (z)}\right)\text{Glaub}[\widehat{\nu }]`$ (50)
where $`\widehat{\nu }`$ is any operator and $`\text{Glaub}[\widehat{\nu }]`$ stands for the Glauber P distribution of $`\widehat{\nu }`$. The first rule Eq.(49) comes from the fact that the coherent state $`|\text{coh}:\psi `$ is an eigenstate of $`\widehat{\psi }(z)`$ with the eigenvalue $`\psi (z)`$. The second rule Eq.(50) involves a functional derivative with respect to the field $`\psi `$, the fields $`\psi `$ and $`\psi ^{}`$ being formally considered as independent variables. Its derivation closely follows the one for single mode fields in . One first uses the following identity:

$$_{\psi (z)}\left[|\text{coh}:\psi \text{coh}:\psi |\right]=(\widehat{\psi }^{}(z)\psi ^{}(z))|\text{coh}:\psi \text{coh}:\psi |$$
(51)
as can be checked from the definition Eq.(36) transposed to the one-dimensional case. Then one integrates by part in the functional integral over the field to convert the derivative over the dyadic $`|\text{coh}:\psi \text{coh}:\psi |`$ into a derivative of the Glauber P distribution. We finally obtain the Fokker-Planck type equation for the fictitious time evolution of the Glauber distribution:

$$_\tau P=EP𝑑z\left[_\psi (F(z)P)+\frac{g}{4}_\psi ^2(\psi ^2P)+\text{c.c.}\right].$$
(52)
The first term in Eq.(52) involves a multiplication of $`P`$ by a functional $`E`$ which is simply the Gross-Pitaevskii energy functional , obtained by replacing in the Hamiltonian the field operator by the c-number field $`\psi `$:

$$E[\{\psi \},\{\psi ^{}\}]=_0^L𝑑z\left[\frac{\mathrm{}^2}{2m}|_z\psi |^2+\frac{1}{2}g|\psi |^4\mu |\psi |^2\right].$$
(53)
If this term was alone in the evolution equation for $`P`$ we could readily integrate it to obtain:

$$P_{\mathrm{class}}[\{\psi \},\{\psi ^{}\}]=e^{\tau E[\{\psi \},\{\psi ^{}\}]}.$$
(54)
We have termed this solution $`P_{\mathrm{class}}`$ as it corresponds to the thermal Boltzmann distribution for a classical field $`\psi `$ with an energy given by $`E[\{\psi \},\{\psi ^{}\}]`$! This is nineteenth century equilibrium physics for fields.

The next term in Eq.(52) can be termed a force term by analogy with the Fokker-Planck equation, as it involves a first order derivative in $`\psi `$. The ‘force’ functional is given by:

$$F[\{\psi \},\{\psi ^{}\}](z)=\frac{1}{2}_{\psi ^{}(z)}E[\{\psi \},\{\psi ^{}\}]=\frac{1}{2}\left[\frac{\mathrm{}^2}{2m}_z^2+g|\psi |^2\mu \right]\psi (z).$$
(55)
At sufficiently short ‘time’ $`\tau `$, that is at high temperature, the effect of the force term during the evolution ‘time’ $`\tau `$ is to shift the field $`\psi (z)`$ from its initial (random) value by the amount $`F(z)\tau `$. One may hope that this shift is negligible for high enough temperature.

Let us calculate this shift for the ideal Bose gas. The field $`\psi `$ is then conveniently expanded on plane waves with momenta given by Eq.(5):

$$\psi (z)=\underset{k_z}{}\alpha _{k_z}\frac{e^{ik_zz}}{L^{1/2}}.$$
(56)
Taking now $`k_z`$ and $`\alpha _{k_z}`$ as coordinate and field variables (rather than $`z`$ and $`\psi (z)`$) we can write:

$$E=\underset{k_z}{}\epsilon _{k_z}\alpha _{k_z}^{}\alpha _{k_z}\text{and}F(k_z)=\frac{1}{2}_{\alpha _{k_z}^{}}E=\frac{1}{2}\epsilon _{k_z}\alpha _{k_z}$$
(57)
where we have introduced the mode eigenenergy $`\epsilon _{k_z}=\mathrm{}^2k_z^2/(2m)\mu `$. The shift $`F(k_z)\tau `$ has a negligible effect on field mode $`k_z`$ if it is small as compared to $`\alpha _{k_z}`$; this leads to the condition

$$k_BT\epsilon _{k_z}.$$
(58)
The classical field approximation $`P_{\mathrm{class}}`$ is therefore an acceptable approximation for the modes with an energy much smaller than $`k_BT`$. The occupation number of such modes can then be obtained from the classical energy equipartition formula: a mode of a complex field corresponds to a one-dimensional harmonic oscillator so it has a mean energy equal to $`k_BT`$, and

$$|\alpha _{k_z}|^2=\frac{k_BT}{\epsilon _{k_z}}.$$
(59)
This formula coincides indeed with the quantum Bose formula Eq.(6) in the limit Eq.(58), that is in the limit of a large occupation number of the mode.

So when can we use the classical field approximation Eq.(54)? The answer depends on the observable quantity we wish to calculate.

For the calculation of the mean energy the classical field approximation is never acceptable: in the absence of energy cut it predicts an infinite mean energy, the well known blackbody radiation catastrophe. To save the situation one has to introduce an energy cut $`\epsilon _{\mathrm{cut}}`$ on the order of $`k_BT`$ to reproduce ‘by hand’ the fact that the Bose formula gives an exponentially small occupation number to modes with eigenenergy much larger than $`k_BT`$. The mean energy is then finite but depends on the precise value of $`\epsilon _{\mathrm{cut}}`$.

The conclusion is different for the calculation of the correlation functions $`g_1`$ and $`g_2`$. We have actually already used the classical field approximation in the derivation of $`g_1`$, see the approximation Eq.(20)! This did not lead to any divergence, a fortunate feature peculiar to the free one-dimensional Bose gas. The classical field predictions for $`g_1`$ and $`g_2`$ therefore do not depend on the energy cut $`\epsilon _{\mathrm{cut}}`$ provided that the energy cut is large enough. For the ideal Bose gas the condition is $`\epsilon _{\mathrm{cut}}|\mu |`$, as the kinetic energy width of the Lorentzian Eq.(20) is $`|\mu |`$; as $`\epsilon _{\mathrm{cut}}k_BT`$ the classical field calculation of $`g_{1,2}`$ requires $`k_BT|\mu |`$: we recover Eq.(16).

The validity conditions of the classical field approximation for the interacting case are more subtle to derive and will be discussed in §3.5.

Finally, the last term in Eq.(52) can be termed a diffusion term by analogy with the Fokker-Planck equation, as it involves second order derivatives in $`\psi `$. The corresponding ‘diffusion’ matrix in point $`z`$, given by:

$$D=\frac{g}{4}\left(\begin{array}{cc}0& \psi ^2\\ \psi ^2& 0\end{array}\right)$$
(60)
is however a non positive matrix: one can check that the field quadrature along $`\psi `$ is squeezed by the ‘diffusion’ while the field quadrature orthogonal to $`\psi `$ gets anti-squeezed. This non-positivity of the ‘diffusion’ matrix makes it impossible to perform a stochastic, Brownian type simulation of Eq.(52), which would have provided an exact numerical solution to the problem.

### 3.3 How to calculate the correlation functions in the classical field approximation

We take here as Glauber P distribution the classical approximation Eq.(54), without introducing any energy cut, and we wish to calculate the stochastic averages Eqs.(40,41). This amounts to calculating ratios of functional integrals over paths parametrized by $`z[0,L]`$.. For example the first order correlation function of the field $`\widehat{\mathrm{\Psi }}(z)`$ in the approximation Eq.(44) is given by

$$g_1(z)=|\varphi _{}(0,0)|^2\frac{𝑑\psi (0)_{\mathrm{\Gamma }_{\psi (0)}}𝒟\psi \psi ^{}(z)\psi (0)e^{\beta E[\{\psi \},\{\psi ^{}\}]}}{𝑑\psi (0)_{\mathrm{\Gamma }_{\psi (0)}}𝒟\psi e^{\beta E[\{\psi \},\{\psi ^{}\}]}}.$$
(61)
In the above formula the functional integrals are performed over all possible closed paths, as the gas is subject to spatially periodic boundary conditions; we have split the functional integral as a regular integral over the value of the path $`\psi (0)`$ in $`z=0`$ and a functional integral over the set $`\mathrm{\Gamma }_{\psi (0)}`$ of all paths starting with the value $`\psi (0)`$ in $`z=0`$ and ending with the same value in $`z=L`$. We explain here how to calculate these functional integrals. The reader not interested in technicalities may jump directly to §3.4.

To calculate functional integrals like Eq.(61) it is of course possible to use a Monte Carlo method to sample the distribution of $`\psi `$, e.g. by representing $`\psi `$ on a finite spatial grid with step $`dz`$ and by evolving $`\psi `$ one time step $`dt`$ after the other according to the stochastic evolution

$$d\psi (z)=dt\left[\frac{\mathrm{}^2}{2m}_z^2+g|\psi |^2\mu \right]\psi (z)+\left(\frac{k_BT}{dz}\right)^{1/2}d\xi (z)$$
(62)
where the $`d\xi (z)`$’s are statistically independent complex noises of variance $`2dt`$.

There exists however a more elegant and much faster solution. One can use the link between path integrals and quantum mechanics propagator put forward by Feynman. The functional integral over the classical complex field then corresponds to a propagator in imaginary time of the quantum mechanical problem of a particle in two dimensions. The point by point analogy between the two problems is specified in the translation table:

| classical field problem |  | quantum mechanical analogy |
| --- | --- | --- |
| path integral | $``$ | quantum propagator |
| abscissa $`z`$ | $``$ | time $`t`$ |
| Re ($`\psi (z)`$) | $``$ | position $`x(t)`$ |
| Im ($`\psi (z)`$) | $``$ | position $`y(t)`$ |
| $`_{\psi (0)=\psi _i}^{\psi (z)=\psi _f}𝒟\psi e^{\beta E[\{\psi \},\{\psi ^{}\}]}`$ | $``$ | $`x_f,y_f|e^{t/\mathrm{}}|x_i,y_i`$ . |

We have to identify the Hamiltonian $``$ of the equivalent quantum mechanics problem. We postulate the following form:

$$=\frac{p_x^2+p_y^2}{2M}+V(x,y).$$
(63)
The imaginary time propagator is then expressed in terms of the path integral :

$$x_f,y_f|e^{t/\mathrm{}}|x_i,y_i=_{x(0)=x_i}^{x(t)=x_f}𝒟x(\tau )_{y(0)=y_i}^{y(t)=y_f}𝒟y(\tau )e^{S[\{x\},\{y\}]/\mathrm{}}$$
(64)
where the action $`S`$ is a functional of the path $`x(\tau ),y(\tau )`$:

$$S[\{x\},\{y\}]=_0^t𝑑\tau \left[\frac{1}{2}M\left(\frac{dx}{d\tau }\right)^2+\frac{1}{2}M\left(\frac{dy}{d\tau }\right)^2+V(x(\tau ),y(\tau ))\right].$$
(65)
One identifies this action with $`\mathrm{}\beta E`$ and one uses the translation table to obtain the values of the parameters of the equivalent quantum mechanics problem: a mass

$$M=\frac{\mathrm{}^3}{mk_BT}$$
(66)
and a potential

$$V(x,y)=\mathrm{}\beta \left[\frac{g}{2}(x^2+y^2)^2\mu (x^2+y^2)\right].$$
(67)
This potential is rotationally invariant, as a consequence of the $`U(1)`$ symmetry of the Hamiltonian Eq.(45). We can therefore classify the eigenstates of $`V`$ with two quantum numbers, an angular momentum $`m`$ and a radial quantum number $`n=0,1,\mathrm{}`$. We call $`\varphi _n^m`$ the corresponding normalized eigenvector with eigenvalue $`\epsilon _n^m`$. As usual the absolute ground state of $``$ is of angular momentum $`m=0`$ and radial quantum number $`n=0`$.

We translate the functional integrals of Eq.(61) into quantum propagators. In particular we note that the integral over $`\psi (0)`$ in Eq.(61) corresponds to an integral over all possible initial coordinates of the particle, that is to a trace over all possible initial quantum states of the particle. We finally obtain:

$$g_1(z)=|\varphi _{}(0,0)|^2\frac{\text{Tr}\left[e^{(Lz)/\mathrm{}}(xiy)e^{z/\mathrm{}}(x+iy)\right]}{\text{Tr}\left[e^{L/\mathrm{}}\right]}$$
(68)
and a similar expression for $`g_2`$.

Physically, as there is no Bose-Einstein condensation along $`z`$, one expects that the length of magnetic guide $`L`$ in the experiment is much larger than any correlation length of the gas. One can then take the thermodynamical limit along $`z`$, putting $`L`$ to infinity while keeping a constant chemical potential $`\mu `$ (this ensures that the mean linear density $`N/L`$ tends to constant) . In this case $`\mathrm{exp}[L/\mathrm{}]`$ becomes proportional to the projector on the absolute ground state of $``$:

$$e^{L/\mathrm{}}e^{L\epsilon _0^{m=0}/\mathrm{}}|\varphi _0^{m=0}\varphi _0^{m=0}|.$$
(69)
The thermodynamical limit approximation greatly simplifies the expressions for the correlation functions, as the trace Tr can be restricted to the ground state of $``$:

$`g_1(z)`$ $`=`$ $`|\varphi _{}(0,0)|^2\varphi _0^{m=0}|(xiy)e^{z(\epsilon _0^{m=0})/\mathrm{}}(x+iy)|\varphi _0^{m=0}`$ (70)
$`g_2(z)`$ $`=`$ $`|\varphi _{}(0,0)|^4\varphi _0^{m=0}|(x^2+y^2)e^{z(\epsilon _0^{m=0})/\mathrm{}}(x^2+y^2)|\varphi _0^{m=0}.`$ (71)
The operator $`x+iy`$ maps the absolute ground state to a state with angular momentum equal to unity. If we restrict for simplicity to the large $`z`$ limit, the operator $`e^z`$ in Eq.(70) becomes proportional to the projector on the ground state $`\varphi _0^{m=1}`$ of $``$ with angular momentum $`m=1`$ so that

$$g_1(z)|\varphi _{}(0,0)|^2a_1\mathrm{exp}(\kappa _1z)$$
(72)
with

$$\kappa _1=(\epsilon _0^{m=1}\epsilon _0^{m=0})/\mathrm{}\text{and}a_1=|\varphi _0^{m=1}|x+iy|\varphi _0^{m=0}|^2.$$
(73)
A similar analysis can be applied to the correlation function $`g_2`$, with the difference that the operator $`x^2+y^2`$ in Eq.(71) maps the absolute ground state to a rotationally invariant state. In the large $`z`$ limit we keep the contributions of the first two eigenstates with $`m=0`$ to obtain

$$g_2(z)g_1^2(0)+|\varphi _{}(0,0)|^4a_2\mathrm{exp}(\kappa _2z)$$
(74)
with

$$\kappa _2=(\epsilon _1^{m=0}\epsilon _0^{m=0})/\mathrm{}\text{and}a_2=|\varphi _1^{m=0}|x^2+y^2|\varphi _0^{m=0}|^2.$$
(75)
### 3.4 Results of the classical field approximation

We wish to calculate the correlation functions $`g_1`$ and $`g_2`$ using the formalism of §3.3. One has then to solve the quantum mechanics equivalent problem of a particle in two-dimension with the Hamiltonian Eq.(63). This Hamiltonian can be diagonalized numerically. We wish to express the results in terms of the linear density of the transversally condensed Bose gas, rather than in terms of the chemical potential $`\mu `$. For a given linear density $`\rho _{\mathrm{lin}}^{(0)}`$ we therefore have to adjust $`\mu `$ in order to satisfy

$$\rho _{\mathrm{lin}}^{(0)}=\widehat{\psi }^{}(0)\widehat{\psi }(0).$$
(76)
The problem can be simplified by an efficient parameterization. We express the coordinates $`x`$ and $`y`$ in the quantum mechanics analogy (which correspond to the real and imaginary part of $`\psi `$) in units of the square root of the linear density. We express the physical length $`z`$ in units of $`\rho _{\mathrm{lin}}^{(0)}\lambda ^2/(2\pi )`$ where $`\lambda `$ is the thermal de Broglie wavelength Eq.(10). We then find that once $`\mu `$ has been eliminated there is a single parameter left in the classical field theory:

$$\chi =\frac{\rho _{\mathrm{lin}}^{(0)}g}{k_BT}\frac{\left(\rho _{\mathrm{lin}}^{(0)}\lambda \right)^2}{2\pi }.$$
(77)
We plot in figure 4 as function of $`\chi `$ the coefficients $`\kappa _{1,2}`$ giving the long distance behavior of the correlation functions $`g_{1,2}`$ as defined in Eqs.(73,75). In the limit of a vanishing $`\chi `$ we recover the results of the ideal Bose gas, Eq.(28), with $`\kappa _1=1/l_c`$, and Eq.(32) leading to $`\kappa _2=2\kappa _1`$. For an increasing interaction strength, $`\chi `$ increases: the coherence length $`1/\kappa _1`$ has a modest increase by up to a factor two; the correlation length of the intensity fluctuations $`1/\kappa _2`$ is dramatically reduced by the atomic interactions and becomes much smaller than the coherence length, a positive point already !

What happens to the amplitude of the intensity fluctuations ? In an ideal ‘atom-laser’ there is no fluctuation of the field intensity; the deviation from this ideal situation can be measured by the ratio of the maximal to the minimal value of $`g_2`$, that is

$$\frac{g_2(0)}{g_1^2(0)}=\frac{\varphi _0^{(m=0)}|(x^2+y^2)^2|\varphi _0^{(m=0)}}{\left(\varphi _0^{(m=0)}|x^2+y^2|\varphi _0^{(m=0)}\right)^2}.$$
(78)
We have plotted this quantity in figure 5. It is equal to two for the ideal Bose gas, as predicted by Eq.(32). It is sharply reduced by the atomic interactions for low values of $`\chi `$ then it slowly goes to unity for large values of $`\chi `$.

Can we understand the origin of the reduction of intensity fluctuations using the quantum mechanics analogy ? According to Eq.(78) this amounts to understanding the fluctuations of the distance of the quantum mechanical particle from the origin in the ground state of $``$. By inspection of Eq.(67) giving the trapping potential seen by the quantum mechanical particle we realize that there are two situations depending on the sign of the chemical potential. For a negative chemical potential, as in the case of the ideal Bose gas, the potential $`V`$ has an absolute minimum in $`x=y=0`$ so that $`|\varphi _0^{m=0}|^2`$ is localized around the origin (see figure 6a): large fluctuations of the intensity of the field are expected. For a positive chemical potential, which is the case for strong enough repulsive interactions, the potential $`V(x,y)`$ has a Mexican hat shape: it has a local maximum at the origin and a global minimum on a finite circle (see figure 6b). The ground state wavefunction then tends to be localized around the circle. The critical regime for the apparition of the Mexican hat potential is such that $`\mu =0`$; we find numerically that this corresponds to

$$\chi _c0.28.$$
(79)
This low value of $`\chi `$ explains why a sharp variation is obtained at the scale of figure 5.

In the large $`\chi `$ limit the quantum mechanical particle will get more deeply bound at the bottom of the Mexican hat potential. Writing Schrödinger’s equation for $`\varphi _n^{m=0,1}`$ in polar coordinates and treating perturbatively the deviation of the Mexican hat plus centrifugal potential from a harmonic approximation we obtain after some algebra the large $`\chi `$’s expansions:

$$\begin{array}{ccccccc}\kappa _1& =& \frac{2\pi }{\rho _{\mathrm{lin}}^{(0)}\lambda ^2}\left[\frac{1}{2}+\frac{1}{2\chi ^{1/2}}+\mathrm{}\right]& & \kappa _2& =& \frac{2\pi }{\rho _{\mathrm{lin}}^{(0)}\lambda ^2}\left[2\chi ^{1/2}+\mathrm{}\right]\\ a_1& =& \rho _{\mathrm{lin}}^{(0)}\left[1+\mathrm{}\right]& & a_2& =& \left(\rho _{\mathrm{lin}}^{(0)}\right)^2\left[\frac{1}{\chi ^{1/2}}+\mathrm{}\right]\\ \frac{g_2(0)}{g_1^2(0)}& =& 1+\frac{1}{\chi ^{1/2}}+\mathrm{}& & \mu & =& \rho _{\mathrm{lin}}^{(0)}g\left[1+\frac{1}{2\chi ^{1/2}}+\mathrm{}\right]\end{array}$$
(80)
We note that the correlation length of the intensity fluctuations $`1/\kappa _2`$ becomes proportional in the large $`\chi `$ limit to the so-called healing length $`\xi `$ of the gas, a crucial parameter in the theory of Bose-Einstein condensates :

$$\kappa _2^1\frac{1}{2}\left(\frac{\mathrm{}^2}{m\rho _{\mathrm{lin}}^{(0)}g}\right)^{1/2}=\frac{\xi }{\sqrt{2}}.$$
(81)
### 3.5 Validity conditions of the classical field approximation

The validity condition of the classical field approximation in general depends on the observable to be calculated. Here the relevant quantities are the first and second order correlation functions $`g_{1,2}`$.

Let us recall briefly what happens in the ideal Bose gas case. The correlations functions are then characterized by a single length, the coherence length $`l_c=1/\kappa _1=2/\kappa _2`$. The plane wave modes of the field contributing to $`g_{1,2}`$ have therefore wave vectors on the order of $`\kappa _1`$ or less; for these modes the validity condition Eq.(58) reads $`k_BT\epsilon _{\kappa _1}`$. The mode eigenenergy $`\epsilon _{\kappa _1}`$ is equal to $`\mathrm{}^2\kappa _1^2/(2m)\mu `$; using the value Eq.(28) of the coherence length $`l_c=1/\kappa _1`$ we arrive at

$$k_BT|\mu |.$$
(82)
This result was already obtained in Eq.(16). It can be rewritten with Eq.(21) as a hierarchy among the three relevant lengths of the problem:

$$\left(\rho _{\mathrm{lin}}^{(0)}\right)^1<\lambda <l_c=\kappa _1^1.$$
(83)
What happens in the interacting regime $`\chi >1`$ ? A difficulty is that the Hamiltonian for the classical field $`\psi `$ in Eq.(53) is no longer quadratic in $`\psi `$ so that it is not straightforward to calculate energy eigenmodes. Can we find some good quadratic approximation to it?

One could think to use the Bogolubov approach . In this approach one identifies the field $`\psi _0`$ minimizing the energy Eq.(53); in the present homogeneous case $`\psi _0`$ is $`z`$-independent, $`\psi _0=(\mu /g)^{1/2}\mathrm{exp}(i\theta _0)`$ where $`\theta _0`$ is an arbitrary constant phase. Then one splits the field as $`\psi (z)=\psi _0+\delta \psi (z)`$. Under the assumption of $`|\delta \psi ||\psi _0|`$ one neglects in Eq.(53) the terms cubic and quartic in $`\delta \psi `$ which leads to a quadratic Hamiltonian that can be diagonalized. This approach is not well suited to the present situation in the large $`L`$ limit; it predicts

$$\frac{\delta \psi ^{}(z)\delta \psi (z)}{|\psi _0|^2}\frac{1}{6}\kappa _1L.$$
(84)
Although the emergence of a coefficient proportional to $`\kappa _1`$ is promising, the Bogolubov quadratization procedure is not justified when $`L`$ exceeds the coherence length of the gas. This is physically not surprising: the finite range of first order field coherence in our one-dimensional geometry is precisely due to large fluctuations of $`\psi (z)`$ away from $`\psi _0`$.

Fortunately it is possible to adapt Bogolubov’s idea taking advantage of the weak intensity fluctuations of the field in the large $`\chi `$ regime. We split the field in a modulus and a phase factor

$$\psi (z)=\rho ^{1/2}(z)e^{i\theta (z)}.$$
(85)
We recall that the intensity $`\rho (z)`$ and the phase $`\theta (z)`$ of a field are Hamiltonian conjugate variables in the same way as $`\psi (z)`$ and $`\psi ^{}(z)`$ are. The intensity $`\rho (z)`$ has only small fluctuations away from the most probable value $`\rho _0=|\psi _0|^2=\mu /g`$:

$$\rho (z)=\rho _0+\delta \rho (z)\text{with}|\delta \rho (z)|\rho _0.$$
(86)
The mean value of $`\rho (z)`$ is simply the linear density of atoms in the transverse ground state $`\rho _{\mathrm{lin}}^{(0)}`$; to lowest order, $`\rho _0`$ can be identified with this mean density. The phase $`\theta `$ has on the contrary large fluctuations away from any fixed constant phase $`\theta _0`$ over distances larger than a few coherence lengths. So contrarily to Bogolubov’s method we do not assume $`\theta (z)\theta _0`$ to be small.

The potential energy density of the field $`g|\psi |^4/2\mu |\psi |^2`$ leads to a constant term and a term quadratic in $`\delta \rho `$. In the kinetic energy density of the field we replace $`\rho (z)`$ by its lowest order approximation giving a non-zero contribution. Up to a constant, this leads to a quadratic approximation for the Hamiltonian Eq.(53):

$$E_{\mathrm{quad}}=_0^L𝑑z\left\{\frac{\mathrm{}^2}{2m}\left[\rho _0\left(_z\theta \right)^2+\frac{1}{4\rho _0}\left(_z\delta \rho \right)^2\right]+\frac{g}{2}\delta \rho ^2(z)\right\}.$$
(87)
One has finally to diagonalize this quadratic Hamiltonian using canonically conjugate variables. An easy way is to expand the field variables on the eigenmodes of the linear time evolution equations:

$`\mathrm{}_t\delta \rho (z)`$ $`=`$ $`_{\theta (z)}E_{\mathrm{quad}}={\displaystyle \frac{\mathrm{}^2}{m}}\rho _0_z^2\theta (z)`$ (88)
$`\mathrm{}_t\theta (z)`$ $`=`$ $`_{\delta \rho (z)}E_{\mathrm{quad}}={\displaystyle \frac{\mathrm{}^2}{4m\rho _0}}_z^2\delta \rho (z)+g\delta \rho (z).`$ (89)
These well known equations have the form of linearized hydrodynamics equations for a superfluid . Eigenmodes are plane waves of wavevector $`k_z`$ obeying Eq.(5) and with eigenenergy

$$\epsilon _{k_z}=\left[\frac{\mathrm{}^2k_z^2}{2m}\left(\frac{\mathrm{}^2k_z^2}{2m}+2\rho _0g\right)\right]^{1/2}.$$
(90)
This is the famous Bogolubov spectrum; the field energy will then formally appear as a sum of decoupled harmonic oscillators, corresponding to an ideal Bose gas of quasi-particles .

We can finally reproduce the reasoning performed above in the ideal Bose gas case to identify the validity condition of the classical field approximation. The temperature must be much larger than the mode eigenenergies $`\epsilon _{k_z}`$ at wavevectors $`k_z=\kappa _1,\kappa _2`$ relevant for the correlation functions $`g_1,g_2`$:

$$k_BT\epsilon _{\kappa _1},\epsilon _{\kappa _2}.$$
(91)
As $`\kappa _2\kappa _1`$ in the large $`\chi `$ regime, the condition involving $`\kappa _2`$ is the most stringent one. From Eq.(81) we find formally the same condition as the ideal Bose gas case

$$k_BT\mu $$
(92)
with now a different expression for the chemical potential, $`\mu \rho _{\mathrm{lin}}^{(0)}g`$. We note that in the large $`\chi `$ regime and the high temperature regime Eq.(92) we have the following hierarchy among the various relevant scales of the problem:

$$\left(\rho _{\mathrm{lin}}^{(0)}\right)^1<\lambda <\kappa _2^1\xi <l_c=\kappa _1^1<a$$
(93)
where the ‘one-dimensional scattering length’ $`a`$ is defined as

$$g=\frac{\mathrm{}^2}{ma}.$$
(94)
The property that the microscopic scale $`a`$ is larger than the other physical lengths of the problem (with the exception of $`L`$ of course!) ensures that the one-dimensional Bose gas is a weakly interacting Bose gas. The fact that this is a necessary condition to a classical field approximation (like the Gross-Pitaevskii equation) is known for three-dimensional Bose gases, with the difference that the three-dimensional scattering length $`a_{3d}`$ has then to be smaller than the macroscopic scales of the gas . In particular the small gaseous parameter is $`a_{3d}/\xi `$ for three-dimensional Bose gases whereas it is $`\xi /a`$ for one-dimensional Bose gases . This is not surprising if one realizes that the weak interaction limit $`g0`$ leads to $`a\mathrm{}`$ in one-dimension.

### 3.6 Comparison with results in the literature

The model of the one-dimensional interacting Bose gas has already been studied by several authors without performing the classical field approximation.

A first line of thought deals with the exact solvable model of a contact interaction potential . To our knowledge a full calculation of $`g_1(z)`$ and $`g_2(z)`$ at finite temperature has not been performed. The value of $`g_2(0)`$ can be found in  at zero temperature. Exact results have been obtained in the strongly interacting regime $`g+\mathrm{}`$: the correlation function $`g_1(z)`$ has been calculated in  at zero temperature; the correlation functions $`g_1(z)`$ and $`g_2(z)`$ at finite temperature have been obtained in . These regimes are different from the finite temperature, weakly interacting Bose gas considered in the present paper. The extension of the calculation of $`g_2(0)`$ to finite temperature seams feasible but we do not know any reference.

A second line of thought is to consider the regime of low intensity fluctuations of the field (our large $`\chi `$ regime): one may take advantage of the weakness of the intensity fluctuations by a linearization of the equations of motion in the hydrodynamic point of view , a Bogolubov type approach , a quadratization of action in a path integral formulation .

We have checked that the classical field predictions in the asymptotic limit $`\chi 1`$ reproduce the results of e.g.  when $`(k_BT)^{1/2}(\rho _{\mathrm{lin}}^{(0)}g)^{1/2}`$. The advantage of the classical field approximation is that it is not restricted to the large $`\chi `$ regime so that the transition from the ideal Bose gas to the strongly interacting case can be studied.

In preparing this article we have discovered that a one-dimensional classical field model very similar to our model has been studied in , with a different physical motivation.

## 4 Conclusion and perspectives

In this paper we have discussed several aspects of our proposal for the production of a continuous ‘atom-laser’ source, consisting in evaporatively cooling an atomic beam in a long magnetic guide. From the classical Boltzmann equation and for expected typical parameters we have estimated the length required to reach the quantum degenerate regime, a few meters for two-dimensional evaporation; the corresponding loss on the atomic flux is only two orders of magnitude.

We have also characterized the coherence properties of the output beam once quantum degeneracy is reached. The gas is expected to experience transverse Bose-Einstein condensation in the guide, leading to a transversally monomode output beam. We have therefore introduced a one-dimensional model for the gas, that we have solved in a classical field approximation. The coherence length of the field along the axis of the magnetic guide can be much larger than the thermal de Broglie wavelength, typically by one order of magnitude. The intensity fluctuations of the beam are very large for an ideal Bose gas but are strongly reduced by repulsive atomic interactions when the healing length $`\xi `$ of the gas becomes smaller than the coherence length of the gas.

Possible extensions of this work are the discussion of superfluidity properties of the ‘atom-laser’ and a complete three-dimensional modeling of the interacting gas in the magnetic guide.

We acknowledge useful discussions with Gora Shlyapnikov, Philippe Grangier, Gordon Baym and Tony Leggett. We thank Alice Sinatra for helpful comments on the manuscript. Laboratoire Kastler Brossel is a unité de recherche de l’École normale supérieure et de l’Université Pierre et Marie Curie, associée au CNRS.


format: text

text:
# Many-spin effects and tunneling splittings in Mn12 magnetic molecules

## I Introduction

Molecular magnets genintro  have proven to be very suitable systems for the study of mesoscopic tunneling effects in magnetic materials. A number of impressive experimental results have been obtained recently, such as thermally-assisted tunneling jumps ; jumps1 , ground state - to - ground state tunneling fe8tun ; sangregorio  and topological phase effects in spin tunneling fe8tun . Among others, the molecular magnet $`\mathrm{Mn}_{12}\mathrm{O}_{12}(\mathrm{CH}_3\mathrm{COO})_{16}(\mathrm{H}_2\mathrm{O})_4`$ (below referred to as Mn12) has received special attention. The effect of resonant magnetization tunneling has been first observed and studied in detailed experiments jumps ; jumps1  on Mn12, and, at present, a substantial amount of reliable experimental data has been collected. Quantitative analysis of these experiments is a challenging theoretical problem involving fundamental issues about tunneling phenomena in mesoscopic magnetic systems. The basic prerequisite for solving this problem is our ability to evaluate accurately and reliably the energy splittings occuring as a result of tunneling between two (quasi)degenerate levels vanhemmen . At present, carefully designed magnetic relaxation experiments at low and ultralow temperatures (tens or hundreds of milliKelvins) can detect jumps ; jumps1  the changes in relaxation time caused by the splittings of order 10-2–10-4 K, and even smaller fe8tun , of order 10-6–10-7 K. The relaxation time data obtained in these experiments give information (although indirect) about the splitting values, so that predictions of the theoretical models can be compared with experimental results.

Conventionally, the molecular magnet Mn12 is considered as a large single spin $`𝒮=10`$ with quasidegenerate levels $`𝒮_z=+M`$ and $`𝒮_z=M`$ split because of tunneling. However, the single-spin Hamiltonian is a phenomenological construct; in reality, this is a many-spin system, consisting of 12 manganese ions coupled by exchange interactions. Here, using Mn12 as a well-studied example, we address the problem of reliable many-spin calculation of the tunneling splittings in molecular magnets. Such a calculation is a very complicated task. For example, the Hilbert space of the spin Hamiltonian describing a molecule of Mn12 consists of 108 levels, while the smallest tunneling splittings in Mn12 are of order of 10-10 Kelvin (as measured in Ref. barbara  for $`m=\pm 10`$). The brute-force direct calculation of tiny tunneling splittings in this system, even for several low-lying states, is beyond the capabilities of modern computers. The general strategy to solve this problem is to truncate the full Hilbert space thus reducing consideration to a much smaller number of relevant energy levels. This idea, implemented in a rather sophisticated way, forms a basis of several approaches for the evaluation of tunneling phenomena, such as quantum Monte-Carlo methods qmcmod , stochastic diagonalization deraedt , and instanton calculations spindyn .

To our knowledge, all calculations of the tunneling splittings in molecular magnets starting from realistic Hamiltonians have employed truncation of the Hilbert space in a much more straightforward, and much less justified manner. High-energy basis states, assumed to be irrelevant, are being explicitely excluded from consideration, and only the low-energy part of the spectrum is being taken into account mn12spl . In the present paper, we calculate tunneling splittings using the many-spin model of Mn12, examining the accuracy and reliability of this straightforward scheme. We demonstrate that, because of strong Dzyaloshinsky-Morya interactions present in Mn12, the splitting values obtained in this way are unreliable. We also consider the sensitivity of the calculated splitting values to variation in the Hamiltonian parameters, and determine the accuracy needed for reliable splittings calculation.

The paper is organized as follows. In Section II, we discuss the 8-spin model of Mn12 and the methods used to calculate tunneling splittings based on this model. We also consider the stability of the results with respect to possible limitations of the model Hamiltonians. In Section III we consider the reasons for the failure of the energy-based truncation scheme in the splittings calculations. Our conclusions can be found in the Summary.

## II  8-spin model of Mn12 and calculations of the tunneling splittings

The cluster Mn12 consists of eight Mn3+ ions having spin 2 and four Mn4+ ions having spin 3/2, coupled by exchange interactions. The total number of spin states in Mn12 is 108, and a corresponding Hamiltonian matrix is rather large to be treated by modern computers. To overcome this difficulty, we can employ the natural hierarchy of interactions present in Mn12. The antiferromagnetic exchange interactions $`J_1220`$ K between Mn3+ and Mn4+ ions are significantly stronger than all the others gat , so corresponding pairs of Mn3+ and Mn4+ ions can be considered as stiff dimers with the total spin $`s=1/2`$, thus giving rise to the 8-spin model of Mn12. The range of validity of the 8-spin model, and the corresponding 8-spin Hamiltonian of Mn12 have been considered in Ref. mn12kat, . After examination of different possible interactions, the following Hamiltonian has been proposed:

$``$ $`=`$ $`J\left({\displaystyle \underset{i}{}}𝐬_i\right)^2J^{}{\displaystyle \underset{k,l}{}}𝐬_k𝐒_lK_z{\displaystyle \underset{i=1}{\overset{4}{}}}\left(S_i^z\right)^2`$
$`+{\displaystyle \underset{i,j}{}}𝐃^{i,j}[𝐬_i\times 𝐒_j].`$
Here, $`𝐒_i`$ and $`𝐬_i`$ are the spin operators for the large spins $`S=2`$ and small dimer spins $`s=1/2`$, correspondingly (the subscript $`i`$ indexes the spins). The first two terms describe isotropic Heisenberg exchange between the spins. The third term describes the single-ion easy-axis anisotropy of large spins. The fourth term represents the antisymmetric Dzyaloshinsky-Morya (DM) interactions in Mn12, where $`𝐃^{i,j}`$ is the Dzyaloshinsky-Morya vector describing the DM-interaction between $`i`$-th small spin and $`j`$-th large spin. Existence of DM-interactions in Mn12 has been suggested in Ref. barbara1 , and their magnitude has been estimated in Ref. mn12kat  based on the neutron scattering data hennion . The molecules of Mn12 possess a fourfold rotational-reflection axis (symmetry $`S_4`$) imposing restrictions on the DM-vectors $`𝐃^{i,j}`$, so that Dzyaloshinsky-Morya interactions can be described by only three parameters $`D_xD_x^{1,8}`$, $`D_yD_y^{1,8}`$, and $`D_zD_z^{1,8}`$.

It has been demonstrated mn12kat  that the above model satisfactorily describes a rather wide range of experimental data, such as the splitting of the neutron scattering peaks, results of EPR measurements and the temperature dependence of magnetic susceptibility. Here, for calculations we use the parameter set A from Ref. mn12kat, :

set A: (2)
$`J=0,J^{}=105\text{ K},K_z=5.69\text{ K}`$
$`D_x=25\text{ K},D_y=0,D_z=1.2\text{ K}.`$
which also gives a good description of the response of Mn12 molecules to a transverse magnetic field (external field applied perpendicular to the easy axis of the molecule). However, this set of parameters should not be considered as being accurately determined, since the amount of the experimental information available is not yet sufficient to achieve particularly reliable parameters. In the Hamiltonian (II), only the fourth term, representing the Dzyaloshinsky-Morya (DM) interactions, can lead to tunneling dm : the first three terms conserve the $`z`$-projection of the total spin $`𝒮_z`$ and can not induce tunneling between levels with different $`𝒮_z`$, while the DM-term mixes levels with different $`𝒮_z`$. In what follows, we will label the energy levels by the value of $`𝒮_z`$. Although it is not an exact quantum number, we can formally consider the DM-interaction as a perturbation, and use perturbation theory terminology.

The following values of the tunneling splittings corresponding to the parameter set (2) have been obtained by the diagonalization of the full Hamiltonian matrix (of the size $`10^4\times 10^4`$) using quadruple precision arithmetics:

$`\mathrm{\Delta }E(\pm 10)`$ $`=`$ $`1.1810^{15}\mathrm{K},\mathrm{\Delta }E(\pm 8)=1.0610^{11}\mathrm{K},`$ (3)
$`\mathrm{\Delta }E(\pm 6)`$ $`=`$ $`3.8710^8\mathrm{K},\mathrm{\Delta }E(\pm 4)=2.0810^6\mathrm{K},`$
$`\mathrm{\Delta }E(\pm 2)`$ $`=`$ $`4.1710^2\mathrm{K}.`$
The splittings for odd values of $`𝒮_z`$ are not shown: they constantly remain at the level of the numerical precision of the calculations (of order of 10-19 K) cray . In Mn12, these splittings should be zero since the fourfold symmetry of the molecule imposes certain restrictions on the symmetry of the spin Hamiltonian and makes some matrix elements vanish. In the single-spin model of Mn12 this property of the spin Hamiltonian is introduced explicitly, by retaining only those operators which possess the required fourfold symmetry. In the many-spin simulations, we obtain the same result independently.

The first question to pose concerns the accuracy of the level splitting evaluation. Parameters of the Hamiltonian are determined with some finite precision, and a small error (say, of the order of several Kelvin) affects the level energy by an amount of order of Kelvin, which is much larger than the very small value of tunneling splitting (of order of 10-12 K). Does it deprive the calculational results of all meaning? To answer this question, we note that the levels $`|𝒮_z=+M`$ and $`|𝒮_z=M`$ are degenerate due to exact symmetry properties of the spin Hamiltonian, and, in the absence of the DM-term, would be degenerate at any value of parameters. Therefore, the tunneling splittings $`\mathrm{\Delta }E_{+M,M}`$ are governed only by the strength of the interaction which breaks the symmetry, i.e. the DM-interaction. If the parameters of the Hamiltonian are determined with reasonably small relative error, and if the numerical calculation is done with sufficient precision, then the relative error of the level splittings will also be small. This conclusion is supported by our calculations: a 10% variation in the Hamiltonian parameters leads to the variation in the splitting values at most by a factor of ten, so that accurate determination of the Hamiltonian parameters is necessary for reliable calculation of the tunneling splittings. If only a logarithmic accuracy in the splitting values is needed, then the 10% uncertainty in the Hamiltonian parameters is sufficient.

However, there is another, much more important source of possible error. The description of the Mn12 molecule by the 8-spin model requires a full, high-precision diagonalization of the Hamiltonian matrix with dimensions $`10^4\times 10^4`$. Solving this problem is rather time-consuming. Matrices of that size can be processed very effectively using Lanczos-type methods, but the application of these methods to the tunneling splitting calculations constitutes quite a difficult problem by itself. A very large number of iterations is needed to achieve the necessary precision and in addition the precision is hard to control when the level separation is very small, so that special techniques are necessary.

Therefore it is natural first to explore another approach, namely, to omit high-energy basis states, retaining only the low-lying part of the spectrum where basis levels have energies less than some threshold value Ecut. This approach has been adopted extensively and in fact, we are not aware of any calculations of tunnel-splitting of magnetic molecules done in a different way: calculations based on both the single-spin and the many-spin model mn12spl  have employed this method. In this paper, we assess the validity of this energy-based truncation approach by considering the dependence of the tunneling splittings $`\mathrm{\Delta }E_{+M,M}`$ for different pairs of degenerate levels $`|𝒮_z=+M`$ and $`|𝒮_z=M`$ on the number of lowest levels $`N_{\mathrm{low}}`$ actually used in calculations (or, in other words, their dependence on the energy threshold $`E_{\mathrm{cut}}`$).

A brief description of the basis states is in order. We first consider the first two exchange terms in the Hamiltonian of Eq. (II) and diagonalize within the manifold of all the 8-spin configurations yielding states with $`𝒮_z=0`$; there are 1286 energy eigenvalues corresponding to eigenvectors with $`𝒮`$ ranging from 0 to 10. The distribution of states is: (10,1), (9,7), (8,24), (7,56), (6,104), (5,164), (4,220), (3,248), (2,232), (1,168), (0,62), where the first number in parenthesis is the value of $`𝒮`$ and the second is the number of levels with this value of $`𝒮`$. With the $`2𝒮+1`$ degeneracies included, there are exactly 10000 states. These are the basis states which are then used to diagonalize the full Hamiltonian, including anisotropy and DM terms.

The initial increase in the number of basis states considered, $`N_{\mathrm{low}}`$, leads to an overall increase in $`\mathrm{\Delta }E_{+M,M}`$ accompanied by oscillations (see Fig. 1). After $`N_{\mathrm{low}}`$ achieves the value of about 700, the oscillations have become small and $`\mathrm{\Delta }E_{+M,M}`$ versus $`N_{\mathrm{low}}`$ exhibits a plateau. This saturation lead in Ref. mn12spl,  to the conclusion that the resulting values give the actual splittings with sufficient accuracy. But this conclusion is wrong. A further increase of the number of levels leads to a resurrection of the oscillations at $`N_{\mathrm{low}}1200`$, with a quite pronounced jump in $`\mathrm{\Delta }E_{+M,M}`$ for $`N_{\mathrm{low}}1700`$. For a larger number of levels, the situation repeats itself: the values of the splittings reach another plateau, then oscillations appear again with a subsequent jump, etc. We have traced this behavior up to $`N_{\mathrm{low}}3000`$, which is already 1/3 of the total number of levels. The observed behavior of $`\mathrm{\Delta }E_{+M,M}`$ is, in our opinion, a very clear signal that energy-based truncation of the Hilbert space is not a good strategy for the computation of tunneling splittings: it gives unreliable results.

The rather sharp jumps in the tunneling splittings as discussed above and illustrated in Fig. 1 are associated with the inclusion of basis states with large $`𝒮`$ values. Because of the selection rule for the DM term ($`𝒮𝒮\pm 1`$), the $`𝒮=10`$ ground state only couples with $`𝒮=9`$ states. States with smaller $`𝒮`$ values affect the splittings more indirectly by coupling with other states which eventually couple to the ground state. While the states with large $`𝒮`$ cause jumps in the splitting values, there are few of them, and the smaller coupling of smaller $`𝒮`$ states still is significant because of the cumulative effect of so many states (see the distribution given above). Therefore, the evaluation of tunneling splittings for a general system possessing strong DM interactions requires consideration of sufficiently large portion of Hilbert space.

It is noteworthy that the same truncation method works rather well for calculations of the energies of well-separated levels. To compare the model against most of the experiments, it suffices to know the positions of the levels with much less precision, usually an error less than 0.1 K is already adequate. This level of precision can be obtained by taking into account $`N_{\mathrm{low}}1000`$ levels (i.e., 1/10 of the total Hilbert space). Even using $`N_{\mathrm{low}}500`$, the error in the level position is less than 1 K even for the states of energy about 60 K. Therefore, the matrix-truncation approach is adequate for fitting the model parameters to experimental data. But the calculations of the tunneling splittings should be done using the full Hamiltonian matrix.

## III Discussion

We have shown that truncating the Hilbert space leads to large errors in the calculated values of tunneling splittings. But actually, any sensible Hamiltonian is inevitably obtained due to some truncation of the Hilbert space. For example, the Hamiltonian (II) can be considered as a result of the two-step procedure yosida : (i) projection of the real many-electron Hamiltonian onto the subspace of suitably chosen single-electron orbital states, yielding a general spin Hamiltonian of the molecule; and (ii) projection of the resulting spin Hamiltonian onto the subspace of the 8-spin model. This procedure is usually justified (at least, at the heuristic level) by invoking some kind of perturbation or WKB-theory arguments, and corresponds to an intuitive idea of different, practically independent energy scales.

However, in the case of the tunneling splittings, we see that very different energy scales significantly affect each other. Why do the same arguments not work if we truncate the 8-spin Hamiltonian? In our opinion, this takes place because the conditions of the applicability of WKB-reasoning (or similar arguments based on perturbation theory) are not satisfied. The spin of the system $`𝒮=10`$ is too small, so that the instanton action vanhemmen  on the trajectories corresponding to the 8-spin model is not large enough. Indeed, for systems with well-separated levels, the quasiclassical approximation usually already works reasonably for a total spin $`𝒮2`$–3. However, as has been demonstrated levine , to apply the same type of arguments to the splitting calculations, the (normalized) instanton action $`S_I`$ should exceed the value of 12. For the model employed in Ref. levine, , this corresponds to the system with a total spin (more exactly, with the total antiferromagnetic vector) of order of several thousand. Thus, the tunneling splittings, in general, appear to be much more sensitive to the method of calculation than the level energies themselves, and conditions for applicability of the conventional WKB-reasoning are considerably more stringent (though for Mn12 they can of course be different from the condition $`S_I>12`$). Qualitatively this agrees with our observations (see Section II). Even a rather severe truncation of the Hilbert space has a minor effect on the level energies, while correct values of the tunneling splittings require a diagonalization of the full Hamiltonian.

Briefly, these arguments can be expressed in a rather obvious form: the 8-spin model is not “macroscopic enough” to justify the truncation of the Hilbert space by some WKB or similar perturbation approach. In this case the intuitive picture of different independent energy scales is misleading.

This conclusion raises important questions, namely, is the 8-spin model, being the result of the truncation of, e.g., 12-spin Hamiltonian, sufficient to predict reliably the tunneling splittings (or, in other words, is the 12-spin model “macroscopic enough” to be truncated)? What is the minimal model allowing the splittings to be calculated correctly? We believe that these are key questions, not only for Mn12 but for the whole class of magnetic molecules. For this purpose, ab-initio calculations of the exchange and anisotropic intramolecular interactions in Mn12 could be very useful. Also, reliable experimental data for the tunneling splittings would obviously be of great value for further development.

## IV Summary

We have calculated the tunneling splittings in Mn12 on the basis of the 8-spin model proposed earlier mn12kat . We have shown that rather accurate knowledge of the Hamiltonian parameters is needed for the accurate splitting calculations; although, for logarithmic accuracy, 10% error in the parameters can be tolerated. Furthermore, we have demonstrated that a reliable calculation of the tunneling splittings for a system with strong DM interactions requires the use of the full Hamiltonian matrix. We have explicitely shown that an energy-based Hilbert space truncation scheme can be successfully used for the determination of the level energies, but leads to erroneous results when applied to the splitting calculations.

## Acknowledgments

This work was partially carried out at the Ames Laboratory, which is operated for the U. S. Department of Energy by Iowa State University under Contract No. W-7405-82 and was supported by the Director of the Office of Science, Office of Basic Energy Research of the U. S. Department of Energy. Support from the Dutch “Stichting Nationale Computer Faciliteiten (NCF)” and the Dutch “Stichting voor Fundamenteel Onderzoek der Materie (FOM)” is gratefully acknowledged.