Bose–Einstein condensation temperature of weakly interacting atoms

Critical temperature is one of the main characteristics of systems experiencing Bose–Einstein condensation (BEC) phase transition. This temperature essentially depends on the system parameters and on the geometry of traps confining atomic clouds. In the present article, we give a survey of the BEC critical temperatures for weakly interacting Bose gases in typical systems, such as homogeneous Bose gas, harmonically trapped Bose gas, bosonic atoms confined by different power-law potentials, and Bose gas in an optical lattice.

We concentrate our attention on the systems for which it is possible to get analytical expressions for critical temperature. This is because having in hands an analytical, even maybe approximate, expression for T_c reveals the explicit role of the system parameters and makes it clear what are the optimal conditions for realizing Bose–Einstein condensation.

In the process of calculations, it is sometimes necessary to resort to nontrivial theoretical methods. One such an approach, allowing for relatively simple and accurate calculations, called self-similar approximation theory, is sketched in the appendix.

Throughout the paper, we use the system of units, where the Planck and Boltzmann constants are set to unity, ${{k}_{\text{B}}}=1$ and $\hbar =1$.

The ideal homogeneous Bose gas in three dimensions, as is well known, exhibits Bose–Einstein condensation at the critical temperature

$$\begin{eqnarray}&&{{T}_{0}}=\frac{2\pi}{m}{{\left[\frac{\rho}{\zeta (3/2)}\right]}^{2/3}},\end{eqnarray} \tag{ 1 }$$

where m is atomic mass, ρ is average particle density, and $\zeta (n)$ is the Riemann zeta function. The question, attracting for long time attention, is how this expression varies under switching on atomic interactions. This problem turned out to be highly nontrivial because the ideal Bose gas and the interacting Bose gas, even with asymptotically weak interactions, pertain to different classes of universality. The ideal Bose gas enjoys the Gaussian universality class, while the interacting three-dimensional Bose gas pertains to the $\mathcal{O}(2)$, or XY, universality class. Although the phase transition in both these cases is of continuous second order, but the physics in the vicinity of the phase transition is of different nature. Close to the phase transition, the physics in an interacting gas is governed by strong fluctuations. As a result, perturbation theory in powers of interaction strength becomes inapplicable, yielding infrared divergences.

One usually considers dilute Bose gas, characterized by the local interaction potential

$$\begin{eqnarray} &&\Phi \left(\mathbf{r}\right)=4\pi \;\frac{{{a}_{s}}}{m}\;\delta \left(\mathbf{r}\right),\end{eqnarray}$$

in which a_s is the s-wave scattering length. One keeps in mind repulsive interactions, with positive scattering length a_s  >  0, since a homogeneous system with a negative scattering length is unstable [1]. The interaction strength of dilute gas is conveniently described by the dimensionless gas parameter

$$\begin{eqnarray}&&\gamma \equiv {{\rho}^{1/3}}{{a}_{s}}.\end{eqnarray} \tag{ 2 }$$

One studies how the critical temperature T_c of an interacting Bose gas shifts from the temperature T₀ of the ideal gas, when switching on atomic interactions. The relative temperature shift is defined as

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}\equiv \frac{{{T}_{\text{c}}}-{{T}_{0}}}{{{T}_{0}}}.\end{eqnarray} \tag{ 3 }$$

For an asymptotically small gas parameter $\gamma \to 0$, the critical temperature shift behaves as [2, 3]

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}\simeq {{c}_{1}}\gamma +\left({{c}_{2}}+c_{2}^{\prime}\ln \gamma \right){{\gamma}^{2}}.\end{eqnarray} \tag{ 4 }$$

The coefficient $c_{2}^{\prime}$ can be found exactly using perturbation theory [3] that gives

$$\begin{eqnarray}&&c_{2}^{\prime}=-\;\frac{64\pi \zeta (1/2)}{3{{\left[\zeta (3/2)\right]}^{5/3}}}=19.7518.\end{eqnarray} \tag{ 5 }$$

While for c₁ and c₂ perturbation theory fails, and one needs more elaborate calculational methods.

There have been numerous attempts of calculating the nonperturbative coefficients c₁ and c₂, employing different techniques, such as Ursell operators and Green functions, renormalization group, the 1/N-expansion in the N-component field theory, and so on, as summarized in the review articles [4, 5]. The results ranged in wide intervals. Thus for c₁ one obtained the values between  −0.95 and 4.7 and for c₂ the values ranging from 4.9 to 101.4, as discussed in [4, 5]. Numerical calculations, using Monte Carlo simulations for three-dimensional lattice $\mathcal{O}(2)$ field theory, give ${{c}_{1}}=1.29\pm 0.05$ [6, 7] and ${{c}_{1}}=1.32\pm 0.02$ [8, 9], while ${{c}_{2}}=75.7\pm 0.4$ [3].

Optimized perturbation theory, advanced in [10, 11], has also been used for calculating the coefficients c₁ and c₂. The main idea of optimized perturbation theory is to define control functions making asymptotic perturbative sequences convergent. Control functions can be introduced in three ways: either by including them into an initial approximation, e.g. into the initial Lagrangian or Hamiltonian, or incorporating them into a sequence transformation, or including them into a change of variables [12]. Optimized perturbation theory has been used for calculating the coefficients c₁ and c₂ by two methods: including control functions into an initial approximation of Lagrangian [13–22] and by introducing them through the Kleinert [23] change of variables [24–26]. Both ways give results close to those of Monte Carlo simulations.

Determining the condensation temperature of a Bose gas can be reformulated as the problem of defining the critical temperature for a three-dimensional N-component field theory [3, 25, 26]. Different N correspond to different physical systems. Thus N  =  0 corresponds to dilute polymer solutions, N  =  1, to the Ising model, N  =  2, to magnetic XY models and to superfluids, and N  =  3 characterizes the Heisenberg model.

When perturbation theory is straightforwardly applied to the calculation of the coefficients c₁ and c₂, the following problem arises. Loop expansion yields asymptotic series in powers of the variable

$$\begin{eqnarray}&&x=\frac{{{\lambda}_{\text{eff}}}}{\sqrt{{{\mu}_{\text{eff}}}}}\;(N+2),\end{eqnarray}$$

in which N is the number of components, ${{\lambda}_{\text{eff}}}$ is an effective coupling parameter, and ${{\mu}_{\text{eff}}}$ is an effective chemical potential [26]. Then the coefficient, say c₁, is represented as an asymptotic series in powers of this variable,

$$\begin{eqnarray}&&{{c}_{1}}(x)\simeq \underset{n}{\sum}\,{{a}_{n}}{{x}^{n}}\quad \quad (x\to 0).\end{eqnarray} \tag{ 6 }$$

In the case of the seven-loop expansion [26], the coefficients a_n for different numbers of the field components N are listed in table 1.

However at T_c, the effective chemical potential tends to zero, ${{\mu}_{\text{eff}}}\to 0$, because of which the variable x tends to infinity, $x\to \infty$. Hence formally we need to find the limit

$$\begin{eqnarray}&&{{c}_{1}}=\underset{x\to \infty}{{\lim}}\,{{c}_{1}}(x).\end{eqnarray} \tag{ 7 }$$

Of course, this limit has no sense being applied directly to series (6). Before taking the limit, it is necessary to extrapolate the series for the asymptotically small variable to its arbitrary values, including asymptotically large values. Such an extrapolation is provided by optimized perturbation theory, as has been done in [13–22, 24–26].

Another method of extrapolation is based on self-similar approximation theory [27–33]. Employing the extrapolation by self-similar factor approximants [34–36], as described in the appendix, we find [37] the values of c₁ summarized in table 2. As is seen, these values are very close to the available Monte Carlo simulations for N  =  2 [6–9] and for N  =  1 and N  =  4 [38].

Table 2. Coefficient c₁ of the critical temperature, for a different number of components N, found by using self-similar factor approximants compared with the available Monte Carlo simulations.

N	c1	Monte Carlo
0	$0.77\pm 0.0$
1	$1.06\pm 0.05$	$1.09\pm 0.09$ [38]
2	$1.29\pm 0.07$	$1.29\pm 0.05$ [6]
		$1.32\pm 0.02$ [8]
3	$1.46\pm 0.08$
4	$1.60\pm 0.09$	$1.60\pm 0.10$ [38]

Note that for the formal limit $N\to \infty$, the coefficient c₁ is known exactly [39], being

$$\begin{eqnarray}&&\underset{N\to \infty}{{\lim}}\,{{c}_{1}}=\frac{8\pi}{3{{\left[\zeta (3/2)\right]}^{4/3}}}=2.328\,473.\end{eqnarray}$$

For the coefficient c₂, one finds the expression that we write here in a slightly different notation, as compared to that of Kastening [26],

$$\begin{eqnarray}&&{{c}_{2}}=\frac{32\pi \zeta (1/2)}{3{{\left[\zeta (3/2)\right]}^{5/3}}}\;\left\{\ln \;\frac{{{\left[\zeta (3/2)\right]}^{2/3}}}{128{{\pi}^{3}}}+\frac{\ln 2}{\sqrt{\pi}}\;\zeta \left(\frac{1}{2}\right)\right. \\ &&\left. \qquad +\frac{\sqrt{\pi}+0.270\,166}{\zeta (1/2)}\;-\;1\right\}\end{eqnarray}$$

$$\begin{eqnarray}&&+\,\,\frac{7}{4}\;c_{1}^{2}+\frac{8\zeta (1/2)}{{{\left[\zeta (3/2)\right]}^{1/3}}}\;\left(3{{d}_{2}}-2{{c}_{1}}\right).\end{eqnarray} \tag{ 8 }$$

In the seven-loop expansion, one has [26]

$$\begin{eqnarray}&&{{d}_{N}}(x)=\frac{1}{{{x}^{2}}}\left(\underset{n=0}{\overset{7}{\sum}}\,{{b}_{n}}{{x}^{n}}+b_{2}^{\prime}{{x}^{2}}\ln x\right).\end{eqnarray} \tag{ 9 }$$

Several first coefficients b_n can be written explicitly as

$$\begin{eqnarray}&&{{b}_{0}}=(N+2){{A}_{0}},\,\,\,\,\,\,\,\,\,\,{{b}_{1}}=\frac{N+2}{24\pi}\;{{A}_{0}},\end{eqnarray}$$

$$\begin{eqnarray}&&{{b}_{2}}=\frac{1-4\ln 6}{576{{\pi}^{2}}}\;{{A}_{0}}=-2.393\,297,\quad \quad \\ &&b_{2}^{\prime}=\frac{{{A}_{0}}}{288{{\pi}^{2}}}=0.776\,158,\end{eqnarray}$$

where

$$\begin{eqnarray}&&{{A}_{0}}=\frac{256{{\pi}^{3}}}{{{\left[\zeta (3/2)\right]}^{4/3}}}=2206.186\,117\,57.\end{eqnarray}$$

These and the higher-order coefficients b_n are given in table 3.

Table 3. Coefficients b_n of the asymptotic expansion for d_N(x), calculated in seven-loop perturbation theory [26] for a different number of components N.

N	0	1	2	3	4
b0	4412.37	6618.56	8824.74	11 030.9	13 237.1
b1	58.5209	87.7814	117.042	146.302	175.563
b2	−2.393 297	−2.393 297	−2.393 297	−2.393 297	−2.393 297
b3	$-1.966\,07\centerdot {{10}^{-2}}$	$-1.437\,53\centerdot {{10}^{-2}}$	$-1.173\,27\centerdot {{10}^{-2}}$	$-1.014\,70\centerdot {{10}^{-2}}$	$-9.089\,97\centerdot {{10}^{-3}}$
b4	$3.308\,62\centerdot {{10}^{-4}}$	$1.786\,49\centerdot {{10}^{-4}}$	$1.195\,96\centerdot {{10}^{-4}}$	$8.962\,15\centerdot {{10}^{-5}}$	$7.191\,20\centerdot {{10}^{-5}}$
b5	$-8.827\,44\centerdot {{10}^{-6}}$	$-3.513\,48\centerdot {{10}^{-6}}$	$-1.924\,34\centerdot {{10}^{-6}}$	$-1.246\,01\centerdot {{10}^{-6}}$	$-8.930\,33\centerdot {{10}^{-7}}$
b6	$2.997\,99\centerdot {{10}^{-7}}$	$8.756\,42\centerdot {{10}^{-8}}$	$3.908\,59\centerdot {{10}^{-8}}$	$2.179\,53\centerdot {{10}^{-8}}$	$1.391\,27\centerdot {{10}^{-8}}$
b7	$-1.196\,71\centerdot {{10}^{-8}}$	$-2.551\,54\centerdot {{10}^{-9}}$	$-9.244\,18\centerdot {{10}^{-10}}$	$-4.424\,57\centerdot {{10}^{-10}}$	$-2.5085\centerdot {{10}^{-10}}$

Again, one needs to extrapolate the asymptotic series (9) to finite values of the variable x and to find an effective limit

$$\begin{eqnarray}&&{{d}_{N}}=\underset{x\to \infty}{{\lim}}\,{{d}_{N}}(x).\end{eqnarray} \tag{ 10 }$$

In table 4, the results for d_N, obtained by Kastening [26], are presented, compared with the available Monte Carlo simulations for N  =  2 [9] and for N  =  1,4 [38]. Thus for Bose–Einstein condensation with N  =  2, the Monte Carlo simulations give ${{c}_{2}}=75.7\pm 0.4$ [9].

Table 4. Parameters d_N entering the coefficients c₂, found by Monte Carlo simulations and optimized perturbation theory, for a different number of components N.

N	dN (Monte Carlo)	dN (OPT) [26]
0		$0.636\pm 0.01$
1	$0.898\pm 0.004$ [38]	$0.893\pm 0.01$
2	$1.059\pm 0.001$ [9]	$1.060\pm 0.01$
3		$1.178\pm 0.01$
4	$1.255\pm 0.006$ [38]	$1.265\pm 0.01$

We may note that the large N limit for d_N is known exactly [26], being

$$\begin{eqnarray}&&\underset{N\to \infty}{{\lim}}\,{{d}_{N}}=\frac{1+2\ln 2}{288{{\pi}^{2}}}\;{{A}_{0}}=1.852\,1396.\end{eqnarray}$$

By Monte Carlo simulations, one can find the critical temperature not only for weak interactions, when the gas parameter is small, $\gamma \to 0$, but also for finite γ, as has been done by Pilati et al [40] for N  =  2, whose results are presented in table 5. In figure 1, we show the relative critical temperature of Bose–Einstein condensation ${{T}_{\text{c}}}/{{T}_{0}}$ as a function of the gas parameter γ, following from equation (4), as compared with the Monte Carlo results from table 5. As is clear, equation (4) is valid only for $\gamma \to 0$, as one could expect. In the simulations of Pilati et al [40], a Bose gas of hard spheres is considered, where a_s corresponds to the hard-sphere diameter. The hard-sphere system is often used as a reference system for realistic fluids, such as liquid helium [41, 42]. When a_s is much shorter than the mean interatomic distance a, the results for the local pseudopotential coincide with those for the hard-sphere system [43]. Moreover, the results for the local potential can be extended [44] to finite values of the ratio a_s/a up to a_s/a  =  0.65, where the fluid freezes [45]. Generally, it is possible to show that pseudopotentials can be employed for a rather accurate modeling of physical systems, including those with nonintegrable interaction potentials [46]. In the case of superfluid ⁴He at saturated vapor pressure, the effective hard-sphere diameter is ${{a}_{s}}=2.203~{\mathring{\rm A}}$, which corresponds to the gas parameter $\gamma =0.5944$ and to the relative critical temperature ${{T}_{\text{c}}}/{{T}_{0}}=0.7$.

Table 5. Relative critical temperature ${{T}_{\text{c}}}/{{T}_{0}}$ as a function of the gas parameter $\gamma$, obtained by Monte Carlo simulations [40] for the three-dimensional homogeneous Bose system.

$\gamma$	0.004 64	0.007 94	0.01	0.0215	0.0464	0.126	0.171	0.215	0.368
${{T}_{\text{c}}}/{{T}_{0}}$	1.0069	1.0091	1.0127	1.0214	1.0351	1.0624	1.0652	1.0627	1.0060

Zoom In Zoom Out Reset image size

Notice that the critical temperature, at asymptotically weak interactions increases with γ. This is because an ideal homogeneous Bose-condensed gas is unstable, while interactions stabilize it [47, 48]. At the same time, strong interactions destroy the condensate, because of which at higher values of γ the critical temperature diminishes.

Many finite quantum systems are well represented as being confined by effective potentials [49]. Most often, one considers a three-dimensional harmonic trapping potential

$$\begin{eqnarray}&&U\left(\mathbf{r}\right)=\frac{m}{2}\left(\omega _{x}^{2}{{x}^{2}}+\omega _{y}^{2}{{y}^{2}}+\omega _{z}^{2}{{z}^{2}}\right).\end{eqnarray} \tag{ 11 }$$

Strictly speaking, genuine Bose–Einstein condensation happens only in an infinite system. For a gas confined by the harmonic potential, this implies [50] the limits

$$\begin{eqnarray}&&N\to \infty,\quad \quad {{\omega}_{0}}\to 0,\quad \quad N\omega _{0}^{3}\to \text{const}.\end{eqnarray} \tag{ 12 }$$

Here N is the number of atoms in the trap and

$$\begin{eqnarray}&&{{\omega}_{0}}\equiv {{\left({{\omega}_{x}}{{\omega}_{y}}{{\omega}_{z}}\right)}^{1/3}}.\end{eqnarray} \tag{ 13 }$$

For finite N there occurs pseudocondensation or quasicondensation [51], which in what follows, for short, will also be called condensation.

The ideal Bose gas in a harmonic trap condenses at the critical temperature

$$\begin{eqnarray}&&{{T}_{0}}={{\left[\frac{N}{\zeta (3)}\right]}^{1/3}}{{\omega}_{0}}.\end{eqnarray} \tag{ 14 }$$

At this temperature, the thermal wavelength is

$$\begin{eqnarray}&&{{\lambda}_{0}}\equiv \sqrt{\frac{2\pi}{m{{T}_{0}}}}=\sqrt{\frac{2\pi}{m{{\omega}_{0}}}}{{\left[\frac{\zeta (3)}{N}\right]}^{1/6}}.\end{eqnarray} \tag{ 15 }$$

The ratio

$$\begin{eqnarray}&&\alpha \equiv \frac{{{a}_{s}}}{{{\lambda}_{0}}}\end{eqnarray} \tag{ 16 }$$

plays the role of a dimensionless coupling parameter. The critical temperature shift in terms of the asymptotically small parameter reads as [52]

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}\simeq {{\bar{c}}_{1}}\alpha +\left({{\bar{c}}_{2}}+\bar{c}_{2}^{\prime}\ln \alpha \right){{\alpha}^{2}}.\end{eqnarray} \tag{ 17 }$$

The first coefficient is known [50, 51] from perturbative calculations, and the coefficient $\bar{c}_{2}^{\prime}$ can also be calculated perturbatively [52],

$$\begin{eqnarray}&&{{\bar{c}}_{1}}=-3.426\,032,\bar{c}_{2}^{\prime}=-\;\frac{32\pi \zeta (2)}{3\zeta (3)}=-45.856\,623.\end{eqnarray} \tag{ 18 }$$

The coefficient ${{\bar{c}}_{2}}$ can be related to lattice simulations in three-dimensional $\mathcal{O}(2)$ field theory [52],

$$\begin{eqnarray}&&{{\bar{c}}_{2}}=21.4-\;\frac{16\pi \zeta (2)}{3\zeta (3)}\;\left[\ln \left(32{{\pi}^{3}}\right)+3.522\,272\right]\\ &&\qquad \ \,+12\zeta (2){{\left[\zeta \left(\frac{3}{2}\right)\right]}^{4/3}}{{d}_{2}},\end{eqnarray} \tag{ 19 }$$

with the same d₂ as in equation (9). Using the Monte Carlo result [9] for d₂, one has

$$\begin{eqnarray}&&{{\bar{c}}_{2}}=-155.0.\end{eqnarray}$$

We may notice the difference with the case of a homogeneous gas, where the linear term is positive, while for the trapped gas, according to equation (18), it is negative. This is because the ideal Bose-condensed gas in a three-dimensional harmonic trap is stable, so that switching on interactions immediately starts destroying the condensate. While, on the contrary, the ideal homogeneous Bose-condensed gas is unstable. Therefore interactions play the dual role, first stabilizing the system and then, when the interactions become sufficiently strong, they start depleting the condensate [5, 48].

The critical temperature shift due to repulsive interactions was studied experimentally for harmonically trapped ⁸⁷Rb atoms [53] and ³⁹K atoms [54, 55], varying the interaction strength by means of Feshbach resonance. The measurements were found to be in good agreement with the first coefficient ${{\bar{c}}_{1}}\approx -3.426$. However the next term was different, as compared with equation (17). Exploring the range $0.001<\alpha <0.06$, all experimental data points were fitted [54, 55] by the second-order polynomial

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}={{b}_{1}}\alpha +{{b}_{2}}{{\alpha}^{2}},\end{eqnarray} \tag{ 20 }$$

with the coefficients

$$\begin{eqnarray}&&{{b}_{1}}={{\bar{c}}_{1}}=-3.426,\quad \quad {{b}_{2}}=46\pm 5.\end{eqnarray}$$

The second-order positive term is due to atomic correlations beyond mean-field picture [56–60]. In order to show how expressions equations (17) and (20) are different, they are presented in figure 2.

Zoom In Zoom Out Reset image size

Traps, confining atoms, can be not only harmonic, but more generally, of power law, having in d-dimensional space the form

$$\begin{eqnarray}&&U\left(\mathbf{r}\right)=\underset{\alpha =1}{\overset{d}{\sum}}\,\frac{{{\omega}_{\alpha}}}{2}|\frac{{{r}_{\alpha}}}{{{l}_{\alpha}}}{{|}^{{{n}_{\alpha}}}}\quad \quad \left({{l}_{\alpha}}\equiv \frac{1}{\sqrt{m{{\omega}_{\alpha}}}}\right).\end{eqnarray} \tag{ 21 }$$

The characteristic trap frequency and length are defined as

$$\begin{eqnarray}&&{{\omega}_{0}}\equiv {{\left(\underset{\alpha =1}{\overset{d}{\prod}}\,{{\omega}_{\alpha}}\right)}^{1/d}}=\frac{1}{ml_{0}^{2}},\quad \quad {{l}_{0}}\equiv {{\left(\underset{\alpha =1}{\overset{d}{\prod}}\,{{l}_{\alpha}}\right)}^{1/d}}=\frac{1}{\sqrt{m{{\omega}_{0}}}},\end{eqnarray} \tag{ 22 }$$

respectively. For what follows, it is convenient to introduce the confining dimension [48, 61]

$$\begin{eqnarray}&&s\equiv \frac{d}{2}+\underset{\alpha =1}{\overset{d}{\sum}}\,\frac{1}{{{n}_{\alpha}}}.\end{eqnarray} \tag{ 23 }$$

Let us consider the ideal Bose gas trapped inside potential (21). In the semiclassical description, we have the density of states

$$\begin{eqnarray}&&\rho (\varepsilon )=\frac{{{\varepsilon}^{s-1}}}{{{\gamma}_{d}} \Gamma (s)}\quad \quad \left(s\geqslant \frac{1}{2}\right),\end{eqnarray} \tag{ 24 }$$

in which

$$\begin{eqnarray}&&{{\gamma}_{d}}\equiv \frac{{{\pi}^{d/2}}}{{{2}^{s}}}\underset{\alpha =1}{\overset{d}{\prod}}\,\frac{\omega _{\alpha}^{1/2+1/{{n}_{\alpha}}}}{ \Gamma \left(1+1/{{n}_{\alpha}}\right)}\quad \quad (d\geqslant 1).\end{eqnarray} \tag{ 25 }$$

In the standard semiclassical approximation, the energy variable ε varies between zero and infinity, as a result of which Bose condensation becomes impossible in some low-dimensional traps [62–64]. However, in quantum systems the energy varies not from zero, but from a finite value ${{\varepsilon}_{0}}$ corresponding to the lowest energy level of the quantum system. By the order of magnitude, it is possible to approximate the minimal energy as

$$\begin{eqnarray}&&{{\varepsilon}_{0}}\approx \frac{k_{\text{min}}^{2}}{2m}\approx \frac{1}{2ml_{0}^{2}}=\frac{{{\omega}_{0}}}{2}.\end{eqnarray}$$

Therefore the range of variation of ε has to be

$$\begin{eqnarray}&&{{\varepsilon}_{0}}\approx \frac{{{\omega}_{0}}}{2}\leqslant \varepsilon <\infty.\end{eqnarray} \tag{ 26 }$$

Taking this into account defines the modified semiclassical approximation [61]. Then the Bose–Einstein function is modified to the form

$$\begin{eqnarray}&&{{g}_{n}}(z)=\frac{1}{ \Gamma (n)}{\int}_{{{u}_{\text{min}}}}^{\infty}\frac{z{{u}^{n-1}}}{{{\text{e}}^{u}}-z}\;\text{d}u,\end{eqnarray} \tag{ 27 }$$

where $z={{\text{e}}^{\beta \mu}}$ is fugacity and the integration is not from zero, but from the finite lower limit

$$\begin{eqnarray}&&{{u}_{\text{min}}}\equiv \frac{{{\omega}_{0}}}{2T}.\end{eqnarray}$$

The value ${{\omega}_{0}}/{{T}_{0}}$, where T₀ is the critical temperature of trapped ideal gas, is assumed to be small.

Bose–Einstein condensation of ideal gas in a power-law trap occurs (see details in [48, 61]) at the temperature

$$\begin{eqnarray}&&{{T}_{0}}={{\left[\frac{{{\gamma}_{d}}N}{{{g}_{s}}(1)}\right]}^{1/s}}.\end{eqnarray} \tag{ 28 }$$

In particular, the condensation temperature formally exists in one-dimensional space, where for an anharmonic trap one has

$$\begin{eqnarray}&&{{T}_{0}}=\frac{\sqrt{\pi}(1-s) \Gamma (s)}{2 \Gamma (1+1/n)}\;N{{\omega}_{0}}\quad \quad (d=1,\;s<1,\;n>2).\end{eqnarray} \tag{ 29 }$$

For harmonic confinement, we come to the expression that coincides with that obtained in purely quantum consideration [65],

$$\begin{eqnarray}&&{{T}_{0}}=\frac{N{{\omega}_{0}}}{\ln (2N)}\quad \quad (s=d=1,\;n=2).\end{eqnarray} \tag{ 30 }$$

And for two- and three-dimensional harmonic traps, we have

$$\begin{eqnarray}&&{{T}_{0}}={{\omega}_{0}}{{\left[\frac{N}{\zeta (d)}\right]}^{1/d}}\quad \quad \left(s=d\geqslant 2,\;{{n}_{\alpha}}=2\right).\end{eqnarray} \tag{ 31 }$$

Recall that Bose–Einstein condensation, strictly speaking, assumes thermodynamic limit, which requires to study the system properties under a large number of atoms. Usually, in traps this number is really large. In the presence of confining potentials, thermodynamic limit is defined differently as compared to homogeneous systems.

The most general definition of thermodynamic limit, valid for any system, is as follows [47, 48]: extensive observable quantities A_N must increase together with the number of particles N, so that

$$\begin{eqnarray}&&N\to \infty,\quad \quad {{A}_{N}}\to \infty,\quad \quad \frac{{{A}_{N}}}{N}\to \text{const}.\end{eqnarray} \tag{ 32 }$$

As an observable, it is possible to take the system energy E_N, so that

$$\begin{eqnarray}&&N\to \infty,\quad \quad {{E}_{N}}\to \infty,\quad \quad \frac{{{E}_{N}}}{N}\to \text{const}.\end{eqnarray} \tag{ 33 }$$

For the considered case, the latter reduces to the limit

$$\begin{eqnarray}&&N\to \infty,\quad \quad {{\gamma}_{d}}\to 0,\quad \quad N{{\gamma}_{d}}\to \text{const}.\end{eqnarray}$$

In particular, for unipower potentials, when ${{n}_{\alpha}}=n$, we have

$$\begin{eqnarray}&&N\to \infty,\quad \quad {{\omega}_{0}}\to 0,\quad \quad N\omega _{0}^{s}\to \text{const}.\end{eqnarray} \tag{ 34 }$$

In this limit, the critical temperature (28), depending on the confining dimension s, behaves as

$$\begin{eqnarray}&&{{T}_{0}}\;\propto \;\frac{1}{{{N}^{1/s-1}}}\to 0\quad \quad (s<1),\end{eqnarray}$$

$$\begin{eqnarray}&&{{T}_{0}}\;\propto \;\frac{1}{\ln N}\to 0\quad \quad (s=1),\end{eqnarray}$$

$$\begin{eqnarray}&&{{T}_{0}}\;\propto \;\text{const}\quad \quad (s>1).\end{eqnarray} \tag{ 35 }$$

This tells us that a well defined condensation temperature exists only for s  >  1.

Moreover, a critical temperature may formally exist, but the system below this temperature looses its stability, which implies that, actually, such a system cannot exist. To check the stability of a Bose-condensed system, we need to study its compressibility

$$\begin{eqnarray}&&{{\kappa}_{T}}=\frac{\text{var}\left(\hat{N}\right)}{\rho TN}=\frac{1}{{{\rho}^{2}}}\left(\frac{\partial \rho}{\partial \mu}\right)=\frac{1}{\rho N}\left(\frac{\partial N}{\partial \mu}\right),\end{eqnarray} \tag{ 36 }$$

in which $\hat{N}$ is the number-of-particle operator and

$$\begin{eqnarray}&&\text{var}\left(\hat{N}\right)\equiv \langle {{\hat{N}}^{2}}\rangle -{{\langle \hat{N}\rangle}^{2}}.\end{eqnarray}$$

An equilibrium system is stable provided that its particle fluctuations are thermodynamically normal [47, 66, 67], such that

$$\begin{eqnarray}&&0\leqslant {{\kappa}_{T}}<\infty.\end{eqnarray} \tag{ 37 }$$

In other words, the stability condition is

$$\begin{eqnarray}&&0\leqslant \frac{\text{var}\left(\hat{N}\right)}{N}<\infty.\end{eqnarray} \tag{ 38 }$$

Calculating the compressibility, we find that, above the condensation phase transition, systems with the confining dimension s  <  1 are unstable, since

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;-{{N}^{1/s-1}}\to -\infty \quad \quad \left(s<1,\;T>{{T}_{0}}\right).\end{eqnarray} \tag{ 39 }$$

While for larger confining dimensions, systems are stable,

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;\text{const}\quad \quad \left(s\geqslant 1,\;T>{{T}_{0}}\right).\end{eqnarray} \tag{ 40 }$$

The situation is even more delicate below the transition temperature T₀, where the systems up to s  =  2 are unstable, with

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;-{{N}^{2/s-1}}\to -\infty \quad \quad (s<1),\end{eqnarray}$$

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;N\to \infty \quad \quad (s=1),\end{eqnarray}$$

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;{{N}^{2/s-1}}\to \infty \quad \quad (1<s<2),\end{eqnarray}$$

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;\ln N\to \infty \quad \quad (s=2).\end{eqnarray} \tag{ 41 }$$

And only for larger s  >  2, they are stable, because

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;\text{const}\quad \quad (s>2).\end{eqnarray} \tag{ 42 }$$

Remembering the definition of the confining dimension (23), we find the stability condition

$$\begin{eqnarray}&&\frac{d}{2}+\underset{\alpha =1}{\overset{d}{\sum}}\,\frac{1}{{{n}_{\alpha}}}>2.\end{eqnarray} \tag{ 43 }$$

For harmonic confinement, when ${{n}_{\alpha}}=2$, we have

$$\begin{eqnarray}&&s=d\quad \quad \left({{n}_{\alpha}}=2\right).\end{eqnarray}$$

Then we see that Bose-condensed systems in one- and two-dimensional harmonic traps are unstable, displaying thermodynamically anomalous particle fluctuations,

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;N\to \infty \quad \quad (s=d=1),\end{eqnarray}$$

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;\ln N\to \infty \quad \quad (s=d=2).\end{eqnarray} \tag{ 44 }$$

And only a three-dimensional harmonic trap can confine a stable Bose-condensed system, with a finite compressibility,

$$\begin{eqnarray}&&{{\kappa}_{T}}\;\propto \;\text{const}\quad \quad (s=d=3).\end{eqnarray} \tag{ 45 }$$

Although Bose-condensed gas is unstable in one- and two-dimensional harmonic traps, it can be stable for other powers of the confining potential. As an example, let us consider the unipower potentials, when ${{n}_{\alpha}}=n$. Then the confining dimension is

$$\begin{eqnarray}&&s=\left(\frac{1}{2}+\frac{1}{n}\right)d.\end{eqnarray}$$

According to condition (43), the ideal Bose-condensed gas is stable if the confining dimension s  >  2. Therefore, the gas is stable in different real-space dimensions d, provided that the potential powers n are limited from above:

$$\begin{eqnarray}&&n<\frac{2d}{4-d}.\end{eqnarray}$$

That is,

$$\begin{eqnarray}&&n<\frac{2}{3}\quad \quad (d=1),\end{eqnarray}$$

$$\begin{eqnarray}&&n<2\quad \quad (d=2),\end{eqnarray}$$

$$\begin{eqnarray}&&n<6\quad \quad (d=3).\end{eqnarray}$$

The existence of the upper potential power n for stability can be understood, if one remembers that the passage from a trapped gas to the uniform gas confined in a box of length L corresponds to the limits

$$\begin{eqnarray}&&n\to \infty,\quad \quad {{l}_{0}}\to \frac{L}{2},\quad \quad s\to \frac{d}{2}.\end{eqnarray}$$

Hence increasing n approaches the system to the uniform case that is known to be unstable [47, 48].

Generally, for stability of a system, it is also necessary that thermal fluctuations be thermodynamically normal [47]. Calculating specific heat, we find [48, 61] that it is positive and finite for all $s\geqslant 1/2$. From the definition of the confining dimension s it follows that it is always not smaller than 1/2. Hence thermal fluctuations are thermodynamically normal for the considered trapped systems. And stability is defined by the behavior of the isothermal compressibility.

Describing Bose gas in power-law traps, we have used the modified semiclassical approximation, accepting the modified function

$$\begin{eqnarray}&&{{g}_{s}}(z)=\frac{1}{ \Gamma (s)}{\int}_{\beta {{\varepsilon}_{0}}}^{\infty}\;\frac{z{{u}^{s-1}}}{{{\text{e}}^{u}}-z}\;\text{d}u,\end{eqnarray} \tag{ 46 }$$

with a finite lower integration limit, instead of the usual Bose–Einstein function

$$\begin{eqnarray}&&g_{s}^{(0)}(z)=\frac{1}{ \Gamma (s)}{\int}_{0}^{\infty}\frac{z{{u}^{s-1}}}{{{\text{e}}^{u}}-z}\;\text{d}u,\end{eqnarray} \tag{ 47 }$$

with the zero lower limit. In the modified variant, the transition temperature is given by equation (28), while in the standard case, taking into account that $g_{s}^{(0)}(1)=\zeta (s)$, it is

$$\begin{eqnarray}&&T_{0}^{(0)}={{\left[\frac{{{\gamma}_{d}}N}{\zeta (s)}\right]}^{1/s}}\quad \quad (s>1).\end{eqnarray} \tag{ 48 }$$

A finite positive value for the latter temperature exists only for s  >  1, since the Rieman zeta function $\zeta (s)$ is finite and positive for s  >  1, finite but negative for $0\leqslant s<1$, and infinite at s  =  1, where $\zeta (s)\to +\infty$, as $s\to 1+0$.

In order to study how the use of the modified function (46) influences the value of the transition temperature, let us consider the critical temperature shift under a large number of atoms,

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{0}}}{T_{0}^{(0)}}\equiv \frac{{{T}_{0}}-T_{0}^{(0)}}{T_{0}^{(0)}},\end{eqnarray}$$

where s  >  1, since finite temperature (48) exists only for this confining dimension. We find [48]

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{0}}}{T_{0}^{(0)}}=\frac{1}{s(s-1) \Gamma (s)\zeta (s)}{{\left(\frac{{{\varepsilon}_{0}}}{{{T}_{0}}}\right)}^{s-1}},\end{eqnarray} \tag{ 49 }$$

which shows that at large N one has

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{0}}}{T_{0}^{(0)}}\;\propto \;\frac{1}{{{N}^{1-1/s}}}\quad \quad (s>1).\end{eqnarray}$$

For harmonic traps, it follows

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{0}}}{T_{0}^{(0)}}=\frac{1}{4\sqrt{\zeta (2)}\;{{N}^{1/2}}}\quad \quad (s=d=2),\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{0}}}{T_{0}^{(0)}}=\frac{1}{48{{\left[\zeta (3)\right]}^{1/3}}\;{{N}^{2/3}}}\quad \quad (s=d=3).\end{eqnarray} \tag{ 50 }$$

The found critical temperature shifts represent finite-size corrections related to the quantum nature of trapped atomic systems. A finite quantum system possesses an energy spectrum with a nonzero energy ${{\varepsilon}_{0}}$ of the lowest level. This has been taken into account in modifying function (47) to (46). However, the quantum nature of trapped systems also plays the role in a different, and even more important, way, as is shown in the following section.

As far as a finite quantum system possesses a spectrum with a nonzero lowest energy level ${{\varepsilon}_{0}}$, then at the Bose–Einstein condensation temperature T_c, the chemical potential tends to ${{\varepsilon}_{0}}$, but not to zero. Therefore the fugacity $z=\exp (\beta \mu )$ tends to $\exp \left({{\varepsilon}_{0}}/{{T}_{\text{c}}}\right)$, but not to one.

In the previous section, we considered how the critical temperature varies when $g_{s}^{(0)}(1)$ is replaced by g_s(1). Now, our aim is to study the variation of the critical temperature under the replacement of $g_{s}^{(0)}(1)$ by $g_{s}^{(0)}\left({{\text{e}}^{\beta {{\varepsilon}_{0}}}}\right)$. Similarly to the previous section, we consider the confining dimensions s  >  1.

Let us denote by T_c the critical temperature calculated with the use of $g_{s}^{(0)}\left({{\text{e}}^{\beta {{\varepsilon}_{0}}}}\right)$ and by T₀ the condensation temperature corresponding to $g_{s}^{(0)}(1)$. For the relative critical temperature shift, we find [48]

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=-\;\frac{(s-1){{\varepsilon}_{0}}}{s(2-s) \Gamma (s)\zeta (s){{T}_{0}}}{{\left(\frac{{{T}_{0}}}{{{\varepsilon}_{0}}}\right)}^{2-s}}\quad \quad (1<s<2),\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=-\;\frac{{{\varepsilon}_{0}}}{2\zeta (2){{T}_{0}}}\;\ln \left(\frac{{{T}_{0}}}{{{\varepsilon}_{0}}}\right)\quad \quad (s=2),\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=-\,\;\frac{\zeta (s-1){{\varepsilon}_{0}}}{s\zeta (s){{T}_{0}}}\quad \quad (s>2).\end{eqnarray} \tag{ 51 }$$

Taking into account that ${{\varepsilon}_{0}}$ varies with increasing N in the same way as ${{\omega}_{0}}$, that is ${{\varepsilon}_{0}}\propto {{N}^{-1/s}}$, one finds

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}\;\propto \;-\;\frac{1}{{{N}^{1-1/s}}}\quad \quad (1<s<2),\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}\;\propto \;-\;\frac{\ln N}{\sqrt{N}}\quad \quad (s=2),\end{eqnarray}$$

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}\;\propto \;-\;\frac{1}{{{N}^{1/s}}}\quad \quad (s>2).\end{eqnarray} \tag{ 52 }$$

As is seen, the corrections here are of different sign as compared to the previous section. They are of the same order of magnitude for 1  <  s  <  2, but

$$\begin{eqnarray}&&|\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}|\gg \frac{ \Delta {{T}_{0}}}{T_{0}^{(0)}}\quad \quad (s\geqslant 2,\;N\to \infty )\end{eqnarray} \tag{ 53 }$$

for large N. Remembering that the ideal Bose-condensed gas is stable only for s  >  2, we come to the conclusion that, for stable systems, the quantum corrections to the critical temperature, found in the present section, are more important than those of the previous section.

For a three-dimensional harmonic trapping potential, one finds [68, 69].

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=-\;\frac{\zeta (2)\bar{\omega}}{2{{\left[\zeta (3)\right]}^{2/3}}{{\omega}_{0}}{{N}^{1/3}}}\quad \quad (s=d=3),\end{eqnarray} \tag{ 54 }$$

with the notation for the average trap frequency

$$\begin{eqnarray}&&\bar{\omega}\equiv \frac{1}{d}\underset{\alpha =1}{\overset{d}{\sum}}\,{{\omega}_{\alpha}}.\end{eqnarray}$$

This can be rewritten as

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=-0.727\,504\left(\frac{{\bar{\omega}}}{\omega}\right)\frac{1}{{{N}^{1/3}}}.\end{eqnarray} \tag{ 55 }$$

For anisotropic traps, $\bar{\omega}>{{\omega}_{0}}$. And in the case of a spherical harmonic trap, $\bar{\omega}={{\omega}_{0}}$.

The semiclassical approximation can also be used for calculating the critical temperature shift under switching on weak interactions of atoms trapped in a harmonic potential. Above the critical temperature, the spectrum of atoms interacting through the local potential

$$\begin{eqnarray} &&\Phi \left(\mathbf{r}\right)={{ \Phi }_{0}}\delta \left(\mathbf{r}\right)\quad \quad \left({{ \Phi }_{0}}\equiv 4\pi \;\frac{{{a}_{s}}}{m}\right),\end{eqnarray}$$

in the semiclassical approximation, is

$$\begin{eqnarray}&&{{\varepsilon}_{k}}\left(\mathbf{r}\right)=\frac{{{k}^{2}}}{2m}+U\left(\mathbf{r}\right)+2{{ \Phi }_{0}}\rho \left(\mathbf{r}\right)\quad \quad \left(T>{{T}_{\text{c}}}\right),\end{eqnarray} \tag{ 56 }$$

where $\rho \left(\mathbf{r}\right)$ is the local density of atoms. At the critical point, the chemical potential behaves as

$$\begin{eqnarray}&&\mu \to 2{{ \Phi }_{0}}\rho (0)\quad \quad \left(T\to {{T}_{\text{c}}}+0\right).\end{eqnarray} \tag{ 57 }$$

It is possible to show [48] that the interaction corrections can be small only for s  >  2 and d  >  2. Therefore in what follows, we take d  =  3 and assume that s  >  2, which yields

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=-{{c}_{3}}(s)\;\frac{{{a}_{s}}}{{{\lambda}_{0}}}\quad \quad (s>2),\end{eqnarray} \tag{ 58 }$$

with the thermal length

$$\begin{eqnarray}&&{{\lambda}_{0}}\equiv \sqrt{\frac{2\pi}{m{{T}_{0}}}}=\sqrt{2\pi}\;{{\left[\zeta (3)\right]}^{1/6}}\frac{{{l}_{0}}}{{{N}^{1/6}}},\end{eqnarray} \tag{ 59 }$$

and with the notation

$$\begin{eqnarray}&&{{c}_{3}}(s)\equiv \frac{4}{s\zeta (s)}\left[\zeta \left(\frac{3}{2}\right)\zeta (s-1)-{{S}_{3}}(s)\right],\\ &&{{S}_{3}}(s)\equiv \underset{m,n=1}{\overset{\infty}{\sum}}\,\frac{n}{{{(mn)}^{3/2}}{{(m+n)}^{s-3/2}}}.\end{eqnarray}$$

For a harmonic trap, when s  =  d  =  3, one gets [50].

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=-3.426\,032\;\frac{{{a}_{s}}}{{{\lambda}_{0}}}.\end{eqnarray} \tag{ 60 }$$

This reproduces the linear term in expansion (17).

In the similar way, it is possible to find the critical temperature shift for other types of interaction potentials, for instance for dipolar interactions [70–72], for which the shift

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}\;\propto \;\frac{{{a}_{D}}}{{{\lambda}_{0}}}\end{eqnarray}$$

is proportional to the dipolar length ${{a}_{D}}\equiv md_{0}^{2}$, with d₀ being dipolar moment.

Note that in the effective thermodynamic limit (34) we have ${{\gamma}_{d}}\propto \omega _{0}^{s}$ and ${{\omega}_{0}}\propto {{N}^{-1/s}}$. For a three-dimensional trap s  =  d  =  3, hence ${{\omega}_{0}}\propto {{N}^{-1/3}}$. Then ${{l}_{0}}\propto {{N}^{1/6}}$, while ${{\lambda}_{0}}\propto \text{const}$.

Quasi-uniform traps of box shape have recently become available [73]. If the trap has the shape of a box of linear length L and volume V  =  L³, then the atomic wave function has to satisfy the boundary conditions

$$\begin{eqnarray}&&\psi \left({{r}_{x}},{{r}_{y}},{{r}_{z}}\right)=0\quad \quad \left({{r}_{\alpha}}=0,L\right).\end{eqnarray} \tag{ 61 }$$

The finite size of the box makes the condensation temperature T_c different from the transition temperature T₀ of an infinite homogeneous system. In the case of ideal gas,

$$\begin{eqnarray}&&{{T}_{0}}=\frac{2\pi}{m}{{\left[\frac{\rho}{\zeta (3/2)}\right]}^{2/3}}\quad \quad (V\to \infty ).\end{eqnarray} \tag{ 62 }$$

For the ideal gas, the shift of the critical temperature, caused by the finite size of the box, reads as [74]

$$\begin{eqnarray}&&\frac{ \Delta {{T}_{\text{c}}}}{{{T}_{0}}}=0.351\,467\;\frac{\ln N}{{{N}^{1/3}}}.\end{eqnarray} \tag{ 63 }$$

However, if the box is strictly rectangular, and the Bose gas is ideal, then its condensation makes the system unstable, because of thermodynamically anomalous fluctuations, similarly to the homogeneous ideal Bose gas. This follows [48] from the number-of-particle variance

$$\begin{eqnarray}&&\text{var}\left(\hat{N}\right)=\frac{0.423}{2{{\pi}^{3}}}\;{{(mT)}^{2}}{{V}^{4/3}}\end{eqnarray} \tag{ 64 }$$

yielding the anomalous isothermal compressibility

$$\begin{eqnarray}&&{{\kappa}_{T}}=0.006\,82\;\frac{{{m}^{2}}T}{{{\rho}^{7/3}}}\;{{N}^{1/3}}.\end{eqnarray} \tag{ 65 }$$

The ideal Bose-condensed gas is unstable either in an infinite homogeneous system or in a finite box-shaped trap. Fortunately, real atomic systems always enjoy interactions that can stabilize the Bose-condensed gas.

Cold atoms can be loaded in different optical lattices created by laser beams [75–78]. The standard periodic potential, formed by the beams, reads as

$$\begin{eqnarray}&&{{V}_{L}}\left(\mathbf{r}\right)=\underset{\alpha =1}{\overset{d}{\sum}}\,{{V}_{\alpha}}{{\sin}^{2}}\left(k_{0}^{\alpha}{{r}_{\alpha}}\right),\end{eqnarray} \tag{ 66 }$$

where ${{\mathbf{k}}_{\mathbf{0}}}$ is the laser wave vector,

$$\begin{eqnarray}&&k_{0}^{\alpha}\equiv \frac{2\pi}{{{\lambda}_{\alpha}}}=\frac{\pi}{{{a}_{\alpha}}},\end{eqnarray}$$

${{\lambda}_{\alpha}}$ is a laser wavelength, and ${{a}_{\alpha}}$ is a lattice spacing in the α-direction. The lattice depth, or barrier, is characterized by the quantity

$$\begin{eqnarray}&&{{V}_{0}}\equiv \frac{1}{d}\underset{\alpha =1}{\overset{d}{\sum}}\,{{V}_{\alpha}}.\end{eqnarray} \tag{ 67 }$$

Atoms, subject to the action of the laser beams, have the recoil energy

$$\begin{eqnarray}&&{{E}_{R}}\equiv \frac{k_{0}^{2}}{2m}\quad \quad \left(k_{0}^{2}=\underset{\alpha =1}{\overset{d}{\sum}}\,{{\left(k_{0}^{\alpha}\right)}^{2}}\right).\end{eqnarray} \tag{ 68 }$$

The atomic field operator can be expanded over Wannier functions that can be chosen to be well localized [79],

$$\begin{eqnarray}&&\hat{\psi}\left(\mathbf{r}\right)=\underset{jn}{\sum}\,{{\hat{c}}_{jn}}{{w}_{n}}\left(\mathbf{r}-{{\mathbf{a}}_{j}}\right),\end{eqnarray} \tag{ 69 }$$

here $j=1,2,\ldots {{N}_{L}}$ enumerates the lattice sites and n is a band index. Substituting this expansion into the Hamiltonian, one usually considers only the lowest energy band. Then, keeping in mind local interactions, characterized by the scattering length a_s, yields the Hubbard Hamiltonian

$$\begin{eqnarray}&&\hat{H}=-J\underset{\langle ij\rangle}{\sum}\,\hat{c}_{i}^{\dagger}{{\hat{c}}_{j}}\;+\;\frac{U}{2}\underset{j}{\sum}\,\hat{c}_{j}^{\dagger}\hat{c}_{j}^{\dagger}{{\hat{c}}_{j}}{{\hat{c}}_{j}}\;+\;{{h}_{0}}\underset{j}{\sum}\,\hat{c}_{j}^{\dagger}{{\hat{c}}_{j}}.\end{eqnarray} \tag{ 70 }$$

The Hamiltonian parameters in three dimensions, resorting to the tight-binding approximation, can be presented [78] as follows. The tunneling parameter reads as

$$\begin{eqnarray}&&J=\frac{3}{4}\left({{\pi}^{2}}-4\right){{V}_{0}}\exp \left(-\;\frac{3{{\pi}^{2}}}{4}\;\sqrt{\frac{{{V}_{0}}}{{{E}_{R}}}}\right).\end{eqnarray} \tag{ 71 }$$

The on-site interaction is described by

$$\begin{eqnarray}&&U=\sqrt{\frac{8}{\pi}}\;{{k}_{0}}{{a}_{s}}{{E}_{R}}{{\left(\frac{{{V}_{0}}}{{{E}_{R}}}\right)}^{3/4}}.\end{eqnarray} \tag{ 72 }$$

And the last term contains

$$\begin{eqnarray}&&{{h}_{0}}=3{{E}_{R}}\;\sqrt{\frac{{{V}_{0}}}{{{E}_{R}}}}.\end{eqnarray} \tag{ 73 }$$

In the presence of Bose–Einstein condensate, global gauge symmetry is broken, which is the necessary and sufficient condition for the condensation [80–82]. The most convenient way for breaking the symmetry is through the Bogolubov shift that here is equivalent to the canonical transformation

$$\begin{eqnarray}&&{{\hat{c}}_{j}}=\sqrt{\nu {{n}_{0}}}\;+\;{{c}_{j}},\end{eqnarray} \tag{ 74 }$$

in which ν is the lattice filling factor

$$\begin{eqnarray}&&\nu \equiv \frac{N}{{{N}_{L}}}=\rho {{a}^{d}},\end{eqnarray} \tag{ 75 }$$

where ρ is average atomic density, a is a mean interatomic distance, and ${{n}_{0}}={{N}_{0}}/N$ is the condensate fraction.

The operator of uncondensed atoms satisfies [78] the properties

$$\begin{eqnarray}&&\langle {{c}_{j}}\rangle =0,\quad \quad \underset{j}{\sum}\,{{c}_{j}}=0.\end{eqnarray} \tag{ 76 }$$

In the Hartree–Fock–Bogolubov approximation, the critical temperature of Bose–Einstein condensation, in d dimensions, becomes [78, 83]

$$\begin{eqnarray}&&{{T}_{\text{c}}}=4\pi \;\frac{d-2}{d}{{\left[\Gamma \left(1+\frac{d}{2}\right)\right]}^{2/d}}J\nu.\end{eqnarray} \tag{ 77 }$$

In three dimensions, this reduces to

$$\begin{eqnarray}&&{{T}_{\text{c}}}=5.064\,222~J\nu.\end{eqnarray} \tag{ 78 }$$

In view of the tunneling parameter (71), we have

$$\begin{eqnarray}&&\frac{{{T}_{\text{c}}}}{\nu {{E}_{R}}}=22.294\;\frac{{{V}_{0}}}{{{E}_{R}}}\;\exp \left(-\;\frac{3{{\pi}^{2}}}{4}\;\sqrt{\frac{{{V}_{0}}}{{{E}_{R}}}}\right).\end{eqnarray} \tag{ 79 }$$

Recall that this expression is valid in the tight-binding approximation, when ${{E}_{R}}\ll {{V}_{0}}$.

The dependence of the critical temperature on the lattice and interaction parameters was studied by Monte Carlo simulations [84]. The linear variation of T_c with the filling factor, agreeing with equation (78), is confirmed for $\nu <6$.

Particle fluctuations in the lattice are thermodynamically normal, with the atomic variance

$$\begin{eqnarray}&&\text{var}\left(\hat{N}\right)=\frac{NT}{\nu U\left({{n}_{0}}+\sigma \right)},\end{eqnarray}$$

where

$$\begin{eqnarray}&&\sigma =\frac{1}{N}\underset{j}{\sum}\,\langle {{c}_{j}}{{c}_{j}}\rangle \end{eqnarray}$$

is the anomalous average. This gives the isothermal compressibility

$$\begin{eqnarray}&&{{\kappa}_{T}}=\frac{\text{var}\left(\hat{N}\right)}{\rho TN}=\frac{1}{\rho \nu \left({{n}_{0}}+\sigma \right)U}.\end{eqnarray} \tag{ 80 }$$

As is seen, if the interaction parameter U tends to zero, then the compressibility diverges, which means that the ideal Bose-condensed gas in a lattice is unstable, while repulsive interactions stabilize it.

Taking into account intersite atomic interactions results in the extended Hubbard model

$$\begin{eqnarray}&&\hat{H}=-J\underset{\langle ij\rangle}{\sum}\,\hat{c}_{i}^{\dagger}{{\hat{c}}_{j}}\;+\;\frac{U}{2}\underset{j}{\sum}\,\hat{c}_{j}^{\dagger}\hat{c}_{j}^{\dagger}{{\hat{c}}_{j}}{{\hat{c}}_{j}}\;\\ &&\qquad +\;\frac{1}{2}\underset{i\ne j}{\sum}\,{{U}_{ij}}\hat{c}_{i}^{\dagger}\hat{c}_{j}^{\dagger}{{\hat{c}}_{i}}{{\hat{c}}_{j}}\;+\;{{h}_{0}}\underset{j}{\sum}\,\hat{c}_{j}^{\dagger}{{\hat{c}}_{j}}.\end{eqnarray}$$

This model is usually studied by means of numerical calculations, such as density-matrix renormalization group [85, 86] and quantum Monte Carlo simulations [87–89]. In the presence of intersite interactions collective phonon excitations arise, which can make an optical lattice unstable although long living [90].

The critical temperature of Bose–Einstein condensation is one of the main characteristics of Bose systems. In the present brief review, we give a survey of critical temperatures for typical systems experiencing Bose–Einstein condensation. These are weakly interacting uniform Bose gases, trapped Bose gases confined in power-law potentials, and bosons in optical lattices. The emphasis is on the cases allowing for the derivation of explicit expressions for the critical temperature. Such explicit expressions, even being approximate, make it possible to better understand the dependence of the temperature on system parameters and to estimate optimal conditions for realizing Bose–Einstein condensation.

In the process of calculating a critical temperature, it is often necessary to employ elaborate mathematical methods. One such a very powerful method is based on self-similar approximation theory. Since this approach can be useful for calculating the critical temperature of different systems, it is sketched in the appendix.

In many cases, the considered system is rather complicated allowing only for the use of some kind of perturbation theory in powers of an asymptotically small parameter, while in real physical systems this parameter can be finite or even very large.

Suppose we are interested in a function f(x) of a real variable x, for which one can get an expansion in powers of this variable, obtaining

$$\begin{eqnarray}&&f(x)\simeq {{f}_{k}}(x)\quad \quad (x\to 0),\end{eqnarray}$$

with

$$\begin{eqnarray}&&{{f}_{k}}(x)={{f}_{0}}(x)\left(1+\underset{n=1}{\overset{k}{\sum}}\,{{a}_{n}}{{x}^{n}}\right),\end{eqnarray}$$

where f₀(x) is given. And assume that we need to know the value of the function for finite, or even large, variable x. That is, we need a method of extrapolating the function from asymptotically small x to the region of finite x. A very powerful and rather simple method is based on self-similar approximation theory [27–33]. Here we briefly mention the main ideas of the theory and its particular approach of extrapolation by self-similar factor approximants [34–36].

A finite series can be treated as a polynomial. By the fundamental theorem of algebra, a polynomial

$$\begin{eqnarray}&&{{P}_{k}}(z)=\underset{n=0}{\overset{k}{\sum}}\,{{a}_{n}}{{z}^{n}}\end{eqnarray}$$

over the field of complex numbers (generally, with a_n and z complex-valued) can be uniquely presented as the product

$$\begin{eqnarray}&&{{P}_{k}}(z)=\underset{j=1}{\overset{k}{\prod}}\,\left({{a}_{j}}+{{b}_{j}}z\right).\end{eqnarray}$$

Therefore, we can write

$$\begin{eqnarray}&&\left(1+\underset{n=1}{\overset{k}{\sum}}\,{{a}_{n}}{{x}^{n}}\right)=\underset{j=1}{\overset{k}{\prod}}\,\left(1+{{b}_{j}}x\right).\end{eqnarray}$$

The self-similar transformation of the linear function gives [34–36].

$$\begin{eqnarray}&&\left(1+{{b}_{j}}x\right)\to {{\left(1+{{A}_{j}}x\right)}^{{{n}_{j}}}}.\end{eqnarray}$$

Then the series f_k(x) transforms to the self-similar factor approximant

$$\begin{eqnarray}&&f_{k}^{\ast}(x)={{f}_{0}}(x)\underset{j=1}{\overset{{{N}_{k}}}{\prod}}\,{{\left(1+{{A}_{j}}x\right)}^{{{n}_{j}}}},\end{eqnarray}$$

where N_k will be given below and the parameters A_j and n_j are defined by the accuracy-through-order procedure, that is, by equating the like-order terms in the small-variable expansion, so that

$$\begin{eqnarray}&&f_{k}^{\ast}(x)\simeq {{f}_{k}}(x)\quad \quad (x\to 0).\end{eqnarray}$$

If the order k is even, then N_k  =  k/2. And the accuracy-through-order procedure yields k equations

$$\begin{eqnarray}&&\underset{j=1}{\overset{k/2}{\sum}}\,{{n}_{j}}A_{j}^{n}={{D}_{n}}(n=1,2,\ldots,k),\end{eqnarray}$$

in which

$$\begin{eqnarray}&&{{D}_{n}}\equiv \frac{{{(-1)}^{n-1}}}{(n-1)!}\;\underset{x\to 0}{{\lim}}\,\;\frac{{{\text{d}}^{n}}}{\text{d}{{x}^{n}}}\;\ln \left(1+\underset{m=1}{\overset{n}{\sum}}\,{{a}_{m}}{{x}^{m}}\right).\end{eqnarray}$$

For instance,

$$\begin{eqnarray}&&{{D}_{1}}={{a}_{1}},\quad \quad {{D}_{2}}=a_{1}^{2}-2{{a}_{2}},\quad \quad {{D}_{3}}=a_{1}^{3}-3{{a}_{1}}{{a}_{2}}+3{{a}_{3}}.\end{eqnarray}$$

In these equations, there are k/2 unknown A_j and k/2 unknown n_j, so that the total number of unknowns equals the number of equations.

If k is odd, then we may set N_k  =  (k  +  1)/2. However, then we again get k equations

$$\begin{eqnarray}&&\underset{j=1}{\overset{(k+1)/2}{\sum}}\,{{n}_{j}}A_{j}^{n}={{D}_{n}}\,\,\,\,\,\,\,\,\,\,(n=1,2,\ldots,k),\end{eqnarray}$$

but now with (k  +  1)/2 unknown A_j and (k  +  1)/2 unknown n_j, which makes k  +  1 unknowns. One of the parameters remains undefined, requiring to impose an additional condition. It is possible to resort to the scaling condition [36], agreeing to measure the parameters A_j in units of one of them, say A₁. This is equivalent to setting A₁  =  1. Thus

$$\begin{eqnarray}{{N}_{k}}=\left\{\begin{array}{*{35}{l}} k/2,~ & ~k=2,4,\ldots \\ (k+1)/2,~ & ~k=3,5,\ldots \end{array},\right. \end{eqnarray}$$

with the scaling condition A₁  =  1 for odd k.

The accuracy-through-order equations with respect to A_j are polynomial, the first equation being of first order, the second, of second order, and so on, with the last equation being of order k. Each polynomial equation of order n possesses n solutions. So that the total number of solutions is $1\times 2\times 3\cdots \times k=k!$. At the same time, from the form of the factor approximants it is evident that each of them is invariant with respect to the $k!$ permutations

$$\begin{eqnarray}&&{{A}_{i}}\to {{A}_{j}},\quad \quad {{n}_{i}}\to {{n}_{j}}.\end{eqnarray}$$

Therefore, the multiplicity of solutions for A_j is trivial, being related to the enumeration of the parameters. Up to this enumeration, the solutions are unique.

If the parameters A_j are given, then we get the linear algebraic equations with respect to n_j. The solutions for the latter have the form

$$\begin{eqnarray}&&{{n}_{j}}=\frac{{{U}_{j}}\left({{N}_{k}}\right)}{V\left({{N}_{k}}\right)}\quad \quad \left(\,j=1,2,\ldots,{{N}_{k}}\right).\end{eqnarray}$$

Here the nominator is the determinant

$$\begin{eqnarray}{{U}_{j}}\left({{N}_{k}}\right)=|\begin{array}{c} {{A}_{1}} & {{A}_{2}} & \ldots & {{D}_{1}} & \ldots & {{A}_{{{N}_{k}}}} \\ A_{1}^{2} & A_{2}^{2} & \ldots & {{D}_{2}} & \ldots & A_{{{N}_{k}}}^{2} \\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \\ A_{1}^{{{N}_{k}}} & A_{2}^{{{N}_{k}}} & \ldots & {{D}_{{{N}_{k}}}} & \ldots & A_{{{N}_{k}}}^{{{N}_{k}}} \end{array}|,\end{eqnarray}$$

with D_i in the jth column, and the denominator is the Vandermonde–Knuth determinant

$$\begin{eqnarray}V\left({{N}_{k}}\right)=|\begin{array}{c} {{A}_{1}} & {{A}_{2}} & \ldots & {{A}_{j}} & \ldots & {{A}_{{{N}_{k}}}} \\ A_{1}^{2} & A_{2}^{2} & \ldots & A_{j}^{2} & \ldots & A_{{{N}_{k}}}^{2} \\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \\ A_{1}^{{{N}_{k}}} & A_{2}^{{{N}_{k}}} & \ldots & A_{j}^{{{N}_{k}}} & \ldots & A_{{{N}_{k}}}^{{{N}_{k}}} \end{array}|.\end{eqnarray}$$

When k is odd, hence A₁  =  1, then V(N_k) is the standard Vandermonde determinant, while when k is even, this determinant by the relation

$$\begin{eqnarray}&&V\left({{N}_{k}}\right)=\left(\underset{j=1}{\overset{{{N}_{k}}}{\prod}}\,{{A}_{j}}\right)\bar{V}\left({{N}_{k}}\right)\end{eqnarray}$$

is connected with the Vandermonde determinant

$$\begin{eqnarray}\bar{V}\left({{N}_{k}}\right)=|\begin{array}{c} 1 & 1 & \ldots & 1 & \ldots & 1 \\ {{A}_{1}} & {{A}_{2}} & \ldots & {{A}_{j}} & \ldots & {{A}_{{{N}_{k}}}} \\ A_{1}^{2} & A_{2}^{2} & \ldots & A_{j}^{2} & \ldots & A_{{{N}_{k}}}^{2} \\ \ldots & \ldots & \ldots & \ldots & \ldots & \ldots \\ A_{1}^{{{N}_{k}}-1} & A_{2}^{{{N}_{k}}-1} & \ldots & A_{j}^{{{N}_{k}}-1} & \ldots & A_{{{N}_{k}}}^{{{N}_{k}}-1} \end{array}|=\underset{1\leqslant i<j\leqslant {{N}_{k}}}{\prod}\,\left({{A}_{j}}-{{A}_{i}}\right).\end{eqnarray}$$

The beauty of the factor approximants is that they, in a finite order, can exactly reconstruct a large class of functions from their asymptotic expansions. This class includes rational, irrational, as well as transcendental functions [34–36]. For example, exactly reproducible are all functions of the type

$$\begin{eqnarray}&&R(x)=\underset{j=1}{\overset{M}{\prod}}\,P_{{{m}_{j}}}^{{{\alpha}_{j}}}(x),\end{eqnarray}$$

where P_m(x) are polynomials and ${{\alpha}_{j}}$ are complex-valued numbers. To exactly reconstruct this function, one needs a factor approximant of the order

$$\begin{eqnarray}&&k=\underset{j=1}{\overset{M}{\sum}}\,{{m}_{j}}.\end{eqnarray}$$

Exactly reproducible are transcendental functions that can be defined as limits of polynomials [36]. For instance, as is easy to check, since

$$\begin{eqnarray}&&\underset{a\to 0}{{\lim}}\,P_{1}^{1/a}(x)={{\text{e}}^{x\,\,\,\,\,\,\,\,\,}}\left({{P}_{1}}(x)=1+ax\right),\end{eqnarray}$$

the exponential function is exactly reproduced in any order $k\geqslant 2$. Keeping in mind such limits, the class of exactly reproducible functions can be denoted as

$$\begin{eqnarray}&&\mathcal{R}=\left\{R(x),\;\lim R(x)\right\}.\end{eqnarray}$$

Self-similar factor approximants extrapolate asymptotic series to finite values of the variables, and even to the variables tending to infinity. Thus, assuming that

$$\begin{eqnarray}&&{{f}_{0}}(x)\simeq A{{x}^{\alpha}}\quad \quad (x\to \infty ),\end{eqnarray}$$

we get the large-variable behavior of the factor approximant as

$$\begin{eqnarray}&&f_{k}^{\ast}(x)\simeq {{B}_{k}}x_{k}^{\beta}\quad \quad (x\to \infty ),\end{eqnarray}$$

with the amplitude

$$\begin{eqnarray}&&{{B}_{k}}=A\underset{j=1}{\overset{{{N}_{k}}}{\prod}}\,A_{j}^{{{n}_{j}}}\end{eqnarray}$$

and the power

$$\begin{eqnarray}&&{{\beta}_{k}}=\alpha +\underset{j=1}{\overset{{{N}_{k}}}{\sum}}\,{{n}_{j}}.\end{eqnarray}$$

When the large-variable behavior of the sought function is known, say

$$\begin{eqnarray}&&f(x)\simeq B{{x}^{\beta}}\quad \quad (x\to \infty ),\end{eqnarray}$$

then it is possible to require that ${{\beta}_{k}}$ be equal to β,

$$\begin{eqnarray}&&\alpha +\underset{j=1}{\overset{{{N}_{k}}}{\sum}}\,{{n}_{j}}=\beta.\end{eqnarray}$$

If it is known that the limit of the sought function is finite, hence $\beta =0$, then one has the condition

$$\begin{eqnarray}&&\alpha +\underset{j=1}{\overset{{{N}_{k}}}{\sum}}\,{{n}_{j}}=0\,\,\,\,\,\,\,\,(\beta =0).\end{eqnarray}$$

Thus, for a finite series f_k(x), one gets a sequence $\left\{\,f_{n}^{\ast}(x)\right\}$ of self-similar factor approximants, with $n=1,2,\ldots k$. If the sequence converges numerically, then as the final answer one can accept the expression $\left[\,f_{k}^{\ast}(x)+f_{k-1}^{\ast}(x)\right]/2$, with the error bar  ±$|\,f_{k}^{\ast}(x)-f_{k-1}^{\ast}(x)|/2$.

The convergence of the sequence $\left\{\,f_{n}^{\ast}(x)\right\}$ can be accelerated in the following way. Being based on the last three factor approximants, one constructs a quadratic spline

$$\begin{eqnarray}&&q(x,t)=a(x)+b(x)t+c(x){{t}^{2}},\end{eqnarray}$$

whose coefficients are determined from the conditions

$$\begin{eqnarray}&&q(x,0)=f_{k-2}^{\ast}(x),\,\,\,\,\,\,q(x,1)=f_{k-1}^{\ast}(x),\,\,\,\,q(x,2)=f_{k}^{\ast}(x).\end{eqnarray}$$

The factor approximant for the spline is

$$\begin{eqnarray}&&{{q}^{\ast}}(x,t)=a(x){{\left[1+A(x)t\right]}^{n(x)}}.\end{eqnarray}$$

As the final answer, one can accept

$$\begin{eqnarray}&&{{q}^{\ast}}(x)=\frac{1}{2}\left[{{q}^{\ast}}(x,2)+{{q}^{\ast}}(x,3)\right],\end{eqnarray}$$

with the error bar  ±$|{{q}^{\ast}}(x,2)-{{q}^{\ast}}(x,3)|/2$.

The found q^*(x) extrapolates the initial asymptotic series to finite values of the variable x, including the limit $x\to \infty$.

This method was used [37] for calculating the limit (7) characterizing the critical-temperature shift for the uniform Bose gas in section 2. The results for c₁ are listed in table 2.