# Liquid-Gas phase transition in Bose-Einstein Condensates

\[

## Abstract

We study the effects of a repulsive three-body interaction on a system of trapped ultra-cold atoms in a Bose-Einstein condensed state. The corresponding $`s`$wave non-linear Schrödinger equation is solved numerically and also by a variational approach. A first-order liquid-gas phase transition is observed for the condensed state up to a critical strength of the effective three-body force.   
PACS 03.75.Fi, 36.40.Ei, 05.30.Jp, 34.10.+x

\]

The experimental evidences of Bose-Einstein condensation (BEC) in magnetically trapped weakly interacting atoms  brought a considerable support to the theoretical research on bosonic condensation. The nature of the effective atom-atom interaction determines the stability of the condensed state: the two-body pseudopotential is repulsive for a positive $`s`$wave atom-atom scattering length and it is attractive for a negative scattering length . The ultra-cold trapped atoms with repulsive two-body interaction undergoes a Bose-Einstein phase-transition to a stable condensed state, in a number of cases found experimentally, as for <sup>87</sup>Rb , for <sup>23</sup>Na  and <sup>7</sup>Li . However, a condensed state of atoms with negative $`s`$wave atom-atom scattering length would be unstable for a large number of atoms .

It was indeed observed in the <sup>7</sup>Li gas , for which the $`s`$wave scattering length is $`a=(14.5\pm 0.4)`$ Å, that the number of allowed atoms in the condensed state was limited to a maximum value between 650 and 1300, which is consistent with the mean-field prediction . An earlier experiment  suggested that the number of atoms in the condensate state was significantly larger than the theoretical predictions with two-body pseudopotential. This is consistent with an addition of a repulsive three-body interaction, which can extend considerably the region of stability for the condensate even for a very weak three-body force.

It was reported in Ref.  that a sufficiently dilute and cold bosonic gas exhibits similar three-body dynamics for both signs of the $`s`$wave atom-atom scattering length and the long-range three-body interaction between neutral atoms is effectively repulsive for either sign of the scattering length. It was suggested that, for a large number of bosons the three-body repulsion can overcome the two-body attraction, and a stable condensate will appear in the trap . Singh and Rokhsar have also observed that above the critical value $`n`$ (which is proportional to their $`\gamma _c`$) the only local minimum is a dense gas state, where the neglect of three-body collisions fails.

In this work, using the mean-field approximation, we investigate the competition between the leading term of an attractive two-body interaction, which is originated from a negative two-atom $`s`$wave scattering length, and a repulsive three-body interaction, which can happen in the Efimov limit  ($`|a|\mathrm{}`$) as discussed in Ref.  <sup>*</sup><sup>*</sup>*The physics of three-atoms in the Efimov limit is discussed in Ref. . This reference extends a previous study of universal aspects of the Efimov effect . The relevance of three-body effects in BEC was also previously reported in Refs. , where it is discussed the stability of the numerical solutions. We show that a kind of liquid-gas phase transition appears inside the Bose condensate.

The Ginzburg - Pitaevskii - Gross (GPG) nonlinear Schrödinger equation (NLSE)  is extended to include the effective potential coming from the three body interaction and then solved numerically in the $`s`$wave channel. The dimensionless parameters are related to the two-body scattering length, the strength of the three-body interaction and the number of atoms in the condensed state. As particularly observed in Ref. , to incorporate all two-body scattering processes in such many particle system, the two-body potential should be replaced by the many-body $`T`$matrix. Usually, at very low energies, this is approximated by the two-body scattering matrix, which is directly proportional to the scattering length $`a`$ . To obtain the desired equation, we first consider the effective Lagrangian density, which describes the condensed wave-function in the Hartree approximation, implying in the GPG energy functional  for the trial wave function $`\mathrm{\Psi }`$:

$``$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}\left(\mathrm{\Psi }^{}{\displaystyle \frac{\mathrm{\Psi }}{t}}{\displaystyle \frac{\mathrm{\Psi }^{}}{t}}\mathrm{\Psi }\right)+{\displaystyle \frac{\mathrm{}^2}{2m}}\mathrm{\Psi }^{}^2\mathrm{\Psi }`$ (2)   
$`{\displaystyle \frac{m}{2}}\omega ^2r^2|\mathrm{\Psi }|^2+_\mathrm{I}.`$   
In our description, the atomic trap is given by a rotationally symmetric harmonic potential, with angular frequency $`\omega `$, and $`_\mathrm{I}`$ gives the effective atom interactions up to three particles.

The effective interaction Lagrangian density for ultra-low temperature bosonic atoms, including two  and three-body effective interaction at zero energy, is written as:

$$_\mathrm{I}=\frac{2\pi \mathrm{}^2a}{m}\left|\mathrm{\Psi }\right|^4\frac{2\lambda _3}{3!}\left|\mathrm{\Psi }\right|^6,$$
(3)   
where $`\lambda _3`$ is the strength of the three-body effective interaction and $`a`$ the scattering length.

The NLSE, which describes the condensed wave-function in the mean-field approximation, is variationally obtained from the effective Lagrangian given in Eq. (2). By considering a stationary solution, $`\mathrm{\Psi }(\stackrel{}{r},t)=e^{i\mu t/\mathrm{}}`$ $`\psi (\stackrel{}{r})`$ where $`\mu `$ is the chemical potential and $`\psi (\stackrel{}{r})`$ is normalized to 1 and by rescaling the NLSE for the $`s`$wave solution, we obtain

$$\left[\frac{d^2}{dx^2}+\frac{1}{4}x^2\frac{|\mathrm{\Phi }(x)|^2}{x^2}+g_3\frac{|\mathrm{\Phi }(x)|^4}{x^4}\right]\mathrm{\Phi }(x)=\beta \mathrm{\Phi }(x)$$
(4)   
for $`a<0`$, where $`x\sqrt{2m\omega /\mathrm{}}r`$ and $`\mathrm{\Phi }(x)N^{1/2}\sqrt{8\pi |a|}r\psi (\stackrel{}{r})`$. The dimensionless parameters, related to the chemical potential and the three-body strength are, respectively, given by $`\beta \mu /\mathrm{}\omega `$ and $`g_3\lambda _3\mathrm{}\omega m^2/(4\pi \mathrm{}^2a)^2`$. The normalization for $`\mathrm{\Phi }(x)`$ reads $`_0^{\mathrm{}}𝑑x|\mathrm{\Phi }(x)|^2=n`$ where the reduced number $`n`$ is related to the number of atoms $`N`$ by $`n2N|a|\sqrt{2m\omega /\mathrm{}}.`$ The boundary conditions  in Eq.(4) are given by $`\mathrm{\Phi }(0)=0`$ and $`\mathrm{\Phi }(x)`$ $`C\mathrm{exp}(x^2/4+[\beta \frac{1}{2}]\mathrm{ln}(x))`$ when $`x\mathrm{}.`$

The above equation, (4) will be treated by numerical procedures for non-linear differential equations, employing the Runge-Kutta (RK) and shooting methods. However, it will be helpful first to consider a variational procedure , using a trial Gaussian wave-function for $`\psi (\stackrel{}{r})`$

$$\psi _{var}(\stackrel{}{r})=\left(\frac{1}{\pi \alpha ^2}\frac{m\omega }{\mathrm{}}\right)^{\frac{3}{4}}\mathrm{exp}\left[\frac{r^2}{2\alpha ^2}\left(\frac{m\omega }{\mathrm{}}\right)\right],$$
(5)   
where $`\alpha `$ is a dimensionless variational parameter. The corresponding root-mean-square radius is proportional to the variational parameter $`\alpha `$, as $`r^2_{var}`$ $`=\alpha ^23\mathrm{}/(2m\omega )`$, while the central density is given by $`\rho _{c,var}(\alpha )`$ $`=\alpha ^3\left(m\omega /\pi \mathrm{}\right)^{3/2}`$. The expression for the total variational energy is given by

$$E_{var}(\alpha )=\mathrm{}\omega N\left[\frac{3}{4}\left(\alpha ^2+\frac{1}{\alpha ^2}\right)\frac{n}{4\sqrt{\pi }\alpha ^3}+\frac{2n^2g_3}{9\sqrt{3}\pi \alpha ^6}\right].$$
(6)   
In the same way, we can obtain the corresponding chemical potential , Eq. ( 4):

$$\mu _{var}(\alpha )=\mathrm{}\omega \left[\frac{3}{4}\left(\alpha ^2+\frac{1}{\alpha ^2}\right)\frac{n}{2\sqrt{\pi }\alpha ^3}+\frac{2n^2g_3}{3\sqrt{3}\pi \alpha ^6}\right].$$
(7)   
The variational solutions of $`E_{var}(\alpha )`$ are given, as a function of $`n`$ and $`g_3`$ (where $`a<0`$ and $`g_3>0`$), by finding the extrema of Eq. (6) with respect to $`\alpha `$.

In Fig. 1, we first illustrate the variational procedure considering an arbitrarily small three-body interaction, chosen as $`g_3=0.005`$. In the upper part of the figure, we show five small plots for the total variational energy $`E`$, in terms of the variational width $`\alpha `$. Each one of the small plots corresponds to particular values of $`n`$. For each number $`n`$ we report the energy of the variational extrema in the lower part of figure 1. In region (I) where the number of atoms is still small, the attractive two body force dominates over the repulsive three-body force and just one minima of the energy as a function of the variational parameter $`\alpha `$ is found. That is also the case for $`g_3=0`$. When the number of atoms is further increased (region (II)) two minima appear in the energy $`E\left(\alpha \right).`$ An unstable maximum is also found between the two minima. The lower energy minimum is stable while the solution corresponding to the smaller $`\alpha `$ is metastable. This solution has a higher density and, consequently, its metastability is justified by the repulsive three-body force acting at higher densities. The minimum number $`n`$ for the appearance of the metastable state is characterized by an inflection point in the energy as a function of $`\alpha `$. The value of $`n`$ at the inflection point corresponds to the beak in the plot of extremum energy versus $`n`$ because for larger $`n`$ three variational solutions are found as depicted in the lower part of figure 1. The attractive two-body and trap potentials dominate the condensed state in the low-density stable phase up to the crossing point (C). At this point, the denser metastable solution becomes degenerate in energy with the lower-density stable solution and a first order phase transition takes place. Since the two solutions differ by their density this transition is analogous to a gas-liquid phase transition for which the density difference between the liquid and the gas is the order parameter. In the variational calculation this occurs at the transition number $`n`$1.3 while the numerical solution of the NLSE gives $`1.2`$. In region (III), we observe two local minima with different energies, a higher-density stable point and a lower-density metastable point. The metastable solution disappears in the beak at the boundary between region (III) and (IV). In regions (III) and (IV) the three-body repulsion stabilized a dense solution against the collapse induced by the two-body attraction. The qualitative features of the variational solution is clearly verified by the numerical solution of the NLSE, as shown by the dashed curve.

In Fig. 2, considering several values of $`g_3`$ (0, 0.012, 0.015, 0.0183 and 0.02), using exact numerical calculations, we present the evolution of some relevant physical quantities $`E,`$ $`\mu ,`$ $`\rho _c`$ and $`r^2`$ as functions of the reduced number of atoms $`n`$. For $`g_3=0`$, our calculation reproduces the result presented in Ref. , with the maximum number of atoms limited by $`n_{max}1.62`$ ($`n`$ is equal to $`|C_{nl}^{3D}|`$ of Ref. ). In the plot for the energy as a function of $`n`$ it is shown that for values of $`g_3>0.0183`$ the phase transition is absent. At $`g_30.0183`$ and $`n1.8`$, the stable, metastable and unstable solutions come to be the same. This corresponds to a critical point associated with a second order phase transition. At this point the derivatives of $`\mu ,`$ $`\rho _c`$ and $`r^2`$ as a function of $`n`$ all diverge.

As shown in the figure, for $`0<g_3<0.0183`$, the density $`\rho _c,`$ the chemical potential $`\mu `$ and the root-mean-squared radius $`r^2`$ present back bendings typical of a first order phase transition. For each $`g_3,`$ the transition point given by the crossing point in the $`E`$ versus $`n`$ corresponds to a Maxwell construction in the diagram of $`\mu `$ versus $`n`$. At this point an equilibrated condensate should undergo a phase transition from the branch extending to small $`n`$ to the branch extending to large $`n.`$ The system should never explore the back bending part of the diagram because as we have seen in figure 1 it is an unstable extremum of the energy. From this figure it is clear that the first branch is associated with large radii, small densities and positive chemical potentials while the second branch presents a more compact configuration with a smaller radius a larger density and a negative chemical potential. This justify the term gas for the first one and liquid for the second one. However we want to stress that both solutions are quantum fluids. With $`g_3=0.012`$ the gas phase happens for $`n<1.64`$ and the liquid phase for $`n>1.64`$. For $`g_3>0.0183`$ all the presented curves are well behaved and a single fluid phase is observed. We also checked that calculations with the variational expression of $`r^2`$, $`\rho _c`$ and $`\mu `$ are in good agreement with the ones depicted in Fig.2, following the same trend shown in Fig.1 for the energy.

Finally, in the lower frame of Fig. 3, we show the phase boundary separating the two phases in the plane defined by $`n`$ and $`g_3`$ and the critical point at $`n1.8`$ and $`g_30.0183`$. In the upper frame, we show the boundary of the forbidden region in the central density versus $`g_3`$ diagram.

To summarize, our calculation presents, at the mean-field level, the consequences of a repulsive three-body effective interaction for the Bose condensed wave-function, together with an attractive two-body interaction. A first-order liquid-gas phase-transition is observed for the condensed state as soon as a small repulsive effective three-body force is introduced. In dimensionless units the critical point is obtained when $`g_30.0183`$ and $`n1.8`$. The characterization of the two-phases through their energies, chemical potentials, central densities and radius were also given for several values of the three-body parameter $`g_3`$.

The results presented in this paper can be relevant to determine a possible clear signature of the presence of repulsive three-body interactions in Bose condensed atoms. It points to a new type of phase transition between two Bose fluids. Because of the condensation of the atoms in a single wave-function this transition may present very peculiar fluctuations and correlations properties. As a consequence, it may fall into a different universality class than the standard liquid-gas phase transition, which are strongly affected by many-body correlations. This question certainly deserves further studies.

Acknowledgments

We thank the organizers of the “International Workshop on Collective Excitations in Fermi and Bose Systems” (São Paulo, Brazil), C. Bertulani, L.F. Canto and M. Hussein, to provide the conditions for stimulating discussions and for the starting of this collaboration. AG, TF and LT also thank N. Akhmediev and M.P. Das for correspondence related to Ref., after conclusion of this work, which has pointed out their previous suggestion of relevance of three-body effects in BEC . This work was partially supported by Fundação de Amparo à Pesquisa do Estado de São Paulo and Conselho Nacional de Desenvolvimento Científico e Tecnológico.