Bose-Einstein condensation in two-dimensional traps

TL;DR

This paper develops an analytical method to accurately determine the critical temperature and condensate fraction for Bose-Einstein condensation in two-dimensional traps, addressing divergence issues and validating results with numerical calculations.

Contribution

It introduces an analytical continuation approach to resolve divergence problems in 2D BEC studies and provides improved expressions for critical parameters.

Findings

Analytical expressions for critical temperature and condensate fraction in 2D traps.

Good agreement between analytical results and numerical calculations.

Validation of grand canonical ensemble calculations for 2D BEC.

Abstract

In two-dimensional traps, since the theoretical study of Bose-Einstein condensation (BEC) will encounter the problem of divergence, the actual contribution of the divergent terms is often estimated in some indirect ways with the accuracy to the leading order. In this paper, by using an analytical continuation method to solve the divergence problem, we obtain the analytical expressions of critical temperature and condensate fraction for Bose gases in a two-dimensional anisotropic box and harmonic trap, respectively. They are consistent with or better than previous studies. Then, we further consider the nonvanishing chemical potential, and obtain the expressions of chemical potential and more precise condensate fraction. These results agree with the numerical calculation well, especially for the case of harmonic traps. The comparison between the grand canonical and canonical ensembles shows that our calculation in the grand canonical ensemble is reliable.

Equations

$$⟨N⟩=\sum _i\frac{1}{z^{-1}{e}^{\beta E_i}-1}=\sum _{\ell =1}^{\infty }{z}^{\ell }\sum _i{e}^{-\ell \beta E_i}=\sum _{\ell =1}^{\infty }K\left(\ell \frac{{\hslash }^2}{2m}\beta \right)z^{\ell },$$

$$K\left(t\right)=\sum _{i=0}^{\infty }{e}^{-{\lambda }_{i}t}.$$

$$K\left(t\right)\approx \frac{1}{{\left(4\pi t\right)}^{d/2}}\sum _{k=0,\frac{1}{2},1,\cdots }^{\infty }{B}_k{t}^k,\left(t\to 0\right)$$

$$K^{\prime }\left(t\right)=\sum _{i=1}^{\infty }{e}^{-\left({\lambda }_i-{\lambda }_0\right)t},$$

$$B_0^{\prime }=B_0,B_{1/2}^{\prime }=B_{1/2},B_1^{\prime }=B_1+{\lambda }_0{B}_0-4\pi ,\cdots$$

$$\lambda \left(n_x,n_y\right)={\pi }^2\left(\frac{n_x^2}{L_x^2}+\frac{n_y^2}{L_y^2}\right)-{\pi }^2\left(\frac{1}{L_x^2}+\frac{1}{L_y^2}\right).\left(n_x,n_y=1,2,3,\cdots \right)$$

$$B_0^{\prime }=S=L_x{L}_y,B_{1/2}^{\prime }=-\sqrt{\pi }\left(L_x+L_y\right),B_1^{\prime }={\pi }^2\left(\frac{L_y}{L_x}+\frac{L_x}{L_y}\right)-3\pi ,\cdots .$$

$$N_e=\sum _{\ell =1}^{\infty }{K}^{\prime }\left(\ell \frac{{\hslash }^2}{2m}\beta \right)z^{\ell }=\sum _{k=0,\frac{1}{2},1,\cdots }^{\infty }\frac{B_k^{\prime }}{{\left(4\pi \right)}^k}{\lambda }^{2k-2}{g}_{1-k}\left(z\right),$$

$$g_{\sigma }\left(z\right)=\frac{1}{\Gamma (\sigma )}{\int }_0^{\infty }\frac{x^{\sigma -1}}{z^{-1}{e}^x-1}dx=\sum _{k=1}^{\infty }\frac{z^k}{k^{\sigma }}$$

$$g_{\sigma}\left(e^{\beta\mu}\right)\approx\left\{\begin{array}[c]{lll}\zeta\left(\sigma\right),&\left(\sigma\geq\frac{3}{2}\right)&\\ -\ln\left(-\beta\mu\right),&\left(\sigma=1\right)&\\ \Gamma\left(1-\sigma\right)\frac{1}{\left(-\beta\mu\right)^{1-\sigma}},\text{ }&\left(\sigma\leq\frac{1}{2}\right)&\text{ \ \ \ }\left(\mu\rightarrow 0\right)\end{array}\right.$$

$$n_e{\lambda }^2\approx -\ln \left(-\beta \mu \right)+\sum _{k=\frac{1}{2},1,\cdots }^{\infty }\Gamma \left(k\right)\frac{B_k^{\prime }}{{\left(4\pi \right)}^{k}S}{\lambda }^{2k}\frac{1}{{\left(-\beta \mu \right)}^k},$$

$$\Gamma \left(\xi \right)={\int }_0^{\infty }{x}^{\xi -1+s}{e}^{-x}dx.\left(s\to 0\right)$$

$$n_e{\lambda }^2$$

$$=-\ln \left(-\beta \mu \right)+\frac{2\pi {\hslash }^2}{mS\left(-\mu \right)}{\int }_0^{\infty }dx{x}^s{e}^{-x}{K}^{\prime }\left(\frac{{\hslash }^2}{2m\left(-\mu \right)}x\right)-\Gamma \left(s\right).$$

$$n_e{\lambda }^2=-\ln \left(-\beta \mu \right)+\frac{2\pi {\hslash }^2}{mS}\Gamma \left(1+s\right){\left(-\mu \right)}^s{\sum }^{\prime }\frac{1}{{\left[E\left(n_x,n_y\right)-\mu \right]}^{1+s}}-\Gamma \left(s\right),$$

$${\sum }^{\prime }\frac{1}{{\left[E\left(n_x,n_y\right)\right]}^{1+s}}={\left(\frac{2mS}{{\pi }^2{\hslash }^2}\right)}^{1+s}\sum _{\left(n_x,n_y\right)\ne \left(1,1\right)}^{\infty }\frac{1}{{\left[\left({\chi }^{-1}{n}_x^2+\chi n_y^2\right)-\left(\chi +{\chi }^{-1}\right)\right]}^{1+s}}$$

$$\displaystyle=\left(\frac{2mS}{\pi^{2}\hbar^{2}}\right)^{1+s}\sum_{p=0}^{\infty}\left(\begin{array}[c]{c}p+s\\ p\end{array}\right)\left[\left(\chi+\chi^{-1}\right)^{p}E_{2}\left(1+s+p;\chi^{-1},\chi\right)-\left(\chi+\chi^{-1}\right)^{-1-s}\right],$$

$$E_2\left(\sigma ;a_1,a_2\right)=\sum _{n_1,n_2=1}^{\infty }\frac{1}{{\left(a_1{n}_1^2+a_2{n}_2^2\right)}^{\sigma }}$$

$${\sum }^{\prime }\frac{1}{{\left[E\left(n_x,n_y\right)\right]}^{1+s}}\approx {\left(\frac{2mS}{{\pi }^2{\hslash }^2}\right)}^{1+s}\frac{\pi }{4}\left(\frac{1}{s}+{\Omega }_2\right),$$

$${\Omega }_2$$

$$+\frac{4}{\pi }\sum _{p=1}^{\infty }\left[{\left(\chi +{\chi }^{-1}\right)}^p{E}_2\left(1+p;{\chi }^{-1},\chi \right)-{\left(\chi +{\chi }^{-1}\right)}^{-1}\right]$$

$$\eta \left(\tau \right)=e^{\frac{\pi i\tau }{12}}\prod _{n=1}^{\infty }\left(1-e^{2n\pi i\tau }\right)$$

$$\Gamma \left(s-n\right)\approx \frac{{\left(-1\right)}^n}{n!}\left[\frac{1}{s}+\psi \left(n+1\right)\right],\left(n=0,1,2,\cdots \right)$$

$$n_e{\lambda }^2=\ln \frac{2mS}{{\pi }^2{\hslash }^2\beta }+{\Omega }_2=\ln \frac{N}{n{\lambda }^2}+\Omega ,$$

$$\Omega =\ln \frac{4}{\pi }+{\Omega }_2$$

$$T_c=\frac{2\pi {\hslash }^2}{m{k}_B}\frac{n}{W\left(N{e}^{\Omega }\right)},$$

$$T_b=T_c\left(N=1000,\chi =1\right)\approx 0.228\frac{2\pi {\hslash }^{2}n}{m{k}_B}$$

$$\frac{N_0^{\left(0\right)}}{N}=1-\frac{1}{n{\lambda }^2}\left(\ln \frac{N}{n{\lambda }^2}+\Omega \right),$$

$$N_0^{\left(1\right)}=\frac{1}{e^{-\beta \mu }-1}\approx \frac{1}{-\beta \mu },$$

$$g_1\left(e^{\beta \mu }\right)\approx -\ln \left(-\beta \mu \right)+\frac{-\beta \mu }{2},$$

$${\sum }^{\prime }\frac{1}{{\left[E\left(n_x,n_y\right)-\mu \right]}^{1+s}}\approx {\sum }^{\prime }\left\{\frac{1}{{\left[E\left(n_x,n_y\right)\right]}^{1+s}}-\frac{\left(1+s\right)\left(-\mu \right)}{{\left[E\left(n_x,n_y\right)\right]}^{2+s}}\right\}.$$

$${\Omega }_3=\sum _{p=0}^{\infty }\left(p+1\right){\left(\chi +{\chi }^{-1}\right)}^p\left[E_2\left(p+2;{\chi }^{-1},\chi \right)-{\left(\chi +{\chi }^{-1}\right)}^{-p-2}\right],$$

Full text

1 Introduction

In recent years, BEC in two-dimensional systems attracts much research. First, the BEC of cold atoms in (quasi)two-dimensional traps has been realized in experiments [1, 2, 3]. Then, more interestingly, the BEC of various bosonic quasiparticles in many-body systems has been widely investigated, such as excitons [4], magnons [5, 6, 7], cavity photons [8, 9, 10], and exciton-polaritons [11, 12, 13]. Many experiments of quasiparticles are realized in two-dimensional traps.

In two dimensions, the realization of BEC is mainly in a box or harmonic trap. In the thermodynamic limit, these two cases have a remarkable difference: As the temperature descends, an ideal Bose gas in a two-dimensional harmonic trap will undergo the BEC phase transition, but in two-dimensional infinite space there is no phase transition. In finite systems, however, their difference becomes small since genuine phase transition cannot occur in either case. In both cases, at low enough temperature, a large fraction of particles will fall into the ground state, so the condensation can still occur. This kind of condensation phenomenon can be observed in experiments.

Unfortunately, there is an obstacle in the theoretical interpretation of the influence of trapping potentials or boundaries on the critical temperature of BEC for ideal Bose gases (We will still use the word ’critical temperature’ in this paper though there is no genuine phase transition in a finite system). In the thermodynamic limit, the critical temperature is determined by the condition that the excited-state population $N_{e}$ is equal to the total particle number $N$ when the chemical potential $\mu=0$. In a finite system, this condition can still be used as an approximate method. However, for trapped gases, the expression of $N_{e}$ is usually divergent at $\mu=0$. This problem is not too serious for a two-dimensional harmonic trap since the leading term is convergent. By neglecting all the other divergent terms, one can obtain the zero-order critical temperature, which is actually the result in the thermodynamic limit and is widely used in the literature [14, 15, 16, 17, 18]. In a two-dimensional box, the problem is particularly serious since all terms of $N_{e}$ are divergent at $\mu=0$. Then even the zero-order result cannot be obtained. In the literature, the critical temperature is determined by, for example, setting a given condensate fraction [19] or numerical calculation [20]. To obtain more precise results, the finite-size effect has been studies for many years, some approximate results of critical temperature and condensate fraction are also presented, often based on the analysis of the nonvanishing ground-state energy in a finite system and only including the leading correction [21, 22, 23]. A systematic method for studying the influence of potentials and boundaries is still lacking.

In this paper, we will use an analytical continuation method to deal with the divergence problem at $\mu=0$, which is based on the heat kernel expansion and $\zeta$-function regularization [24]. First, we will show that the divergence can be removed by a general treatment, and the analytical expressions for critical temperature and condensate fraction for ideal Bose gases in a two-dimensional anisotropic box or harmonic trap are presented, respectively. These results are consistent with or better than the previous studies. Then, more precisely, $\mu=0$ does not exactly hold below the transition point in a finite system, but the divergence problem makes it difficult to solve the chemical potential. We will show that our method is applicable to this problem, and we will give the analytic expressions of the chemical potential and the more precise condensate fraction, respectively. These results agree with the numerical calculation well, especially for the harmonic traps. In addition, to check the influence of the fluctuation in the grand canonical ensemble, we compare the condensate fraction in the grand canonical and canonical ensembles. The comparison indicates that the difference between these two ensembles is very small for particle number $N\sim 10^{3}$.

The paper is organized as follows. In section 2, we discuss the BEC of an ideal Bose gas in a two-dimensional rectangle box. The analytical expressions of the critical temperature, the condensate fraction, and the chemical potential are obtained. In section 3, we discuss the Bose gas in a two-dimensional anisotropic harmonic trap. The first-order correction to the critical temperature, and the analytical expressions of condensate fraction and chemical potential are obtained. They agree with the numerical results very well. In section 4, we give a comparison between the grand canonical and canonical ensembles to show the influence of fluctuation in the grand canonical ensemble. The conclusion and some discussion are presented in section 5. A kind of the Epstein $\zeta$-function is used in our calculation, so we give its asymptotic expansion in Appendix A.

2 Two-dimensional rectangle box

The main tool used in this paper is the heat kernel expansion. In the grand canonical ensemble, the average particle number of an ideal Bose gas can be expanded as

[TABLE]

where $z=e^{\beta\mu}$ is the fugacity, $\beta=1/\left(k_{B}T\right)$ with $k_{B}$ denoting the Boltzmann constant, $\left\{E_{i}\right\}$ is the single-particle energy spectrum, which is proportional to the spectrum $\left\{\lambda_{i}\right\}$ of the Laplacian operator $D=-\nabla^{2}+\left(2m/\hbar^{2}\right)V\left(\mathbf{x}\right)$, $E_{i}=\left(\hbar^{2}/2m\right)\lambda_{i}$, and $K\left(t\right)$ denotes the global heat kernel of the operator $D$ [25, 26, 27]

[TABLE]

For small $t$, the heat kernel expansion of $K\left(t\right)$ has the asymptotic form [25, 26, 27]

[TABLE]

where $d$ is the spatial dimension and $B_{k}$ $\left(k=0,1/2,1,\cdots\right)$ are the heat kernel coefficients. Thus, eq. (1) expresses the average particle number of the Bose gas as a series of global heat kernels.

In the thermodynamic limit, the critical temperature of BEC is determined by the condition that the excited-state population $N_{e}$ equals the total particle number $N$ at $\mu=0$. In a finite system, although genuine phase transitions cannot occur, we can expect to obtain the critical temperature by the same condition as an approximation.

The excited-state population is easy to find from eq. (1) by excluding the ground-state contribution. Furthermore, the transition occurring at $\mu=0$ means that the ground-state energy should be zero, so we need to shift the energy spectrum so that the ground-state energy vanishes. In other words, we will replace the heat kernel eq. (2) by

[TABLE]

in which the ground-state contribution is excluded. Therefore, for the two-dimensional case, the corresponding heat kernel coefficients change to

[TABLE]

In the following, we will consider a Bose gas in a two-dimensional rectangle box of length sides $L_{x}$ and $L_{y}$ with Dirichlet boundary conditions. The shifted spectrum is

[TABLE]

According to eq. (5) and the usual heat kernel coefficients [28], the heat kernel coefficients for $K^{\prime}\left(t\right)$ are

[TABLE]

Replacing the $K\left(t\right)$ in eq. (1) by $K^{\prime}\left(t\right)$, we can obtain the excited-state population as

[TABLE]

where

[TABLE]

is the Bose-Einstein integral, and $\lambda=\sqrt{2\pi\beta}\hbar/\sqrt{m}$ is the mean thermal wavelength. In eq. (8) we have replaced $\left\langle N_{e}\right\rangle$ by $N_{e}$ for simplicity.

In eq. (8), the heat kernel coefficient $B_{k}^{\prime}$ has a dimension of $\left[L^{2-2k}\right]$. If we denote the characteristic length scale of the system as $\bar{L}$, $B_{k}^{\prime}$ will be roughly proportional to $\bar{L}^{2-2k}$, just as in eq. (7). Therefore, eq. (8) is in fact a series of $\lambda/\bar{L}$.

2.1 Critical temperature $T_{c}$

The critical temperature of BEC is determined by $N_{e}=N$ at $\mu=0$. In eq. (8), $N_{e}$ is expressed as a series of a small parameter $\lambda/\bar{L}$, so usually the higher-order terms are just small corrections. However, when $\mu\rightarrow 0$, since the asymptotic behavior of the Bose-Einstein integral is

[TABLE]

where $\zeta\left(\sigma\right)=\sum_{n=1}^{\infty}n^{-\sigma}$ is the Riemann zeta function, every term in eq. (8) is divergent, and the divergence becomes more severe in the higher orders. As a result, it will not work to truncate this series at any finite order. To overcome this divergence problem, we will use an analytical continuation method with the help of the heat kernel expansion and $\zeta$-function regularization [24], in which all the terms in the series are considered.

First, substituting the leading term in the asymptotic expansion of the Bose-Einstein integral eq. (10) into eq. (8) gives

[TABLE]

where $n_{e}=N_{e}/S$ is the number density of excited-state particles. We hope to express the divergent sum in eq. (11) by the heat kernel. For this purpose, introduce a regularization parameter $s$ which will be set to [math] at the end of the calculation in the gamma function

[TABLE]

Eq. (11) becomes

[TABLE]

In the last line we have replaced the divergent series by the heat kernel $K^{\prime}\left(t\right)$ according to the heat kernel expansion.

Then, by the definition of heat kernel eq. (4), we can perform the integral in eq. (13),

[TABLE]

where the prime on the sum $\sum\nolimits^{\prime}$ denotes that the ground state is excluded. Since the transition occurs at $\mu=0$, by neglecting the chemical potential $\mu$ in the denominator, the sum in eq. (14) becomes

[TABLE]

where we have introduced a shape factor $\chi=L_{x}/L_{y}$, and \left(\begin{array}[c]{c}n\\ k\end{array}\right)=\frac{n!}{k!\left(n-k\right)!} is the binomial coefficient,

[TABLE]

is the Epstein $\zeta$-function. By use of eq. (72) in Appendix A, when $s\rightarrow 0$, eq. (17) is divergent and its asymptotic form is

[TABLE]

where

[TABLE]

is a parameter only related to the shape factor $\chi$, $\gamma\approx 0.5772$ is the Euler constant, $\psi\left(z\right)=\Gamma^{\prime}\left(z\right)/\Gamma\left(z\right)$ is the digamma function, and

[TABLE]

is the Dedekind $\eta$-function. Since for $s\rightarrow 0$,

[TABLE]

the divergent term of $s$ from eq. (19) and that from the term with $\Gamma\left(s\right)$ are exactly canceled.

Finally, eq. (14) becomes

[TABLE]

where we have introduced

[TABLE]

for simplicity. In eq. (23), all of the divergent terms of $\mu$ are also canceled, and the final result is fully analytical, so the critical temperature is

[TABLE]

where $W\left(z\right)$ is the Lambert $W$ function, satisfying $z=W\left(ze^{z}\right)$.

Eq. (25) gives the influence of the particle number and the shape of box on the critical temperature. In fig. 1 we plot the relation between critical temperature and $\chi$ at fixed density of particles. It shows that the anisotropy lowers the critical temperature. In this and the following figures, the temperature is rescaled to $T/T_{b}$, where

[TABLE]

is the critical temperature for $N=1000$ in a square box.

There are many studies on the BEC in cavities, most of them concentrate on the three-dimensional cases [29, 30]. For two-dimensional boxes, in [21], the authors give a relation between the critical temperature and particle number, which is similar to eq. (23) but with $\Omega=0$. In ref. [20], the authors discuss the property of an ideal Bose gas in a square box in both the grand canonical ensemble and canonical ensemble in details. Their research is based on numerical calculation, and obtain an expression of critical temperature by fitting the numerical solution. By taking $\chi=1$ so that $\Omega=-1.0468$ in eq. (25), our result will go back to the square box case. The relation between the critical temperature and particle number given by eq. (25) and refs. [20] and [21] are shown in fig. 2. Our result agrees with the numerical calculation in ref. [20] quite well.

2.2 Condensate fraction and chemical potential

In the above discussion, the chemical potential $\mu$ is assumed to be zero at the transition point. It implies that $\mu=0$ holds for $T<T_{c}$ just like in the thermodynamic limit case. Under this assumption, the condensate fraction can be directly obtained from eq. (23):

[TABLE]

which will be called the zero-order condensate fraction in this paper.

The chemical potential $\mu$ cannot be exactly zero at $T<T_{c}$ in a finite system, but because of the divergence problem, directly solving $\mu$ is difficult, especially near the transition point. When $T\ll T_{c}$, $\mu$ can be approximate to $-k_{B}T/N_{0}^{\left(0\right)}$, but this approximation is invalid for $T\sim T_{c}$ since $N_{0}^{\left(0\right)}=0$ at the transition point.

The discussion in the above section provides a way to avoid the divergence, so we can solve $\mu$ by the similar way. Specifically, accurate to $\mu^{1}$, we will add three more terms in eq. (14) to obtain the expression of total particle number: the contribution from the ground-state particles

[TABLE]

the next-to-leading term in the asymptotic expansion of the Bose-Einstein integral in the first term

[TABLE]

and the first-order contribution of $\mu$ in the sum of energy spectrum

[TABLE]

In the right-hand side of this equation, the first sum has been given in eq. (19); the second sum is analytical at $s=0$, so we can directly set $s=0$ in it. By introducing a parameter only related to $\chi$,

[TABLE]

we can express the asymptotic expansion of eq. (30) at $s\rightarrow 0$ as

[TABLE]

Thus, eq. (14) with the additional terms becomes

[TABLE]

where the divergent terms of $s$ have also been canceled. The term $-1/2$ in the parentheses in the second term can be neglected, which means that the contribution from the second term of $g_{1}\left(e^{\beta\mu}\right)$ in eq. (29) is much smaller than that from the second term in the right-hand side of eq. (30). After neglecting this small term, we can solve the chemical potential as

[TABLE]

where

[TABLE]

is the chemical potential at the transition point.

In fig. 3 we plot the relation between the chemical potential and temperature given by eq. (34) for different $N$. The result for $\mu$ in the literature is rare, and we include the numerical results in the figure for comparison.

The first-order condensate fraction $N_{0}^{\left(1\right)}/N$ in eq. (28) is straightforward from eq. (34). In fig. 4 we show the relation between the condensate fraction and temperature for different $N$. We can find that the zero-order condensate fraction $N_{0}^{\left(0\right)}/N$ vanishes at the transition point as expected.

3 Two-dimensional anisotropic harmonic trap

The harmonic trap is the most commonly used trap in BEC experiments and also in the theoretical research. In fact, the thermodynamic properties of Bose gases in two-dimensional harmonic traps can be exactly obtained [31, 32]. On the other hand, due to the divergence problem at the transition point, the critical temperature of BEC in a two-dimensional harmonic trap is often approximately regarded as the thermodynamic-limit value [14, 15, 16, 17, 18]. In the following we will remove the divergence, and give the analytical forms of the critical temperature, the condensate fraction, and the chemical potential.

Consider an ideal Bose gas trapped in an anisotropic harmonic potential

[TABLE]

The single-particle energy spectrum has the form

[TABLE]

where $\omega_{0}=\sqrt{\omega_{x}\omega_{y}}$ and

[TABLE]

where we have introduced $\kappa=\omega_{x}/\omega_{y}$ for convenience, and the ground-state energy has been shifted to [math]. Consequently, the exact solution and the asymptotic expansion of the global heat kernel are

[TABLE]

where $\sum\nolimits^{\prime}$ still represents that the ground state is excluded in the sum, and the expansion coefficients are

[TABLE]

In such a trap, the excited-state population of an ideal Bose gas is

[TABLE]

3.1 Critical temperature $T_{c}$

To determine the critical temperature, we need to know the value of eq. (41) at $\mu=0$. However, under this condition, except the first term of eq. (41), all the other ones are divergent. This divergence can also be removed by the method used in last section.

First, substituting the leading term of the asymptotic expansion of the Bose-Einstein integral eq. (10) into eq. (41) and replacing the gamma function by eq. (12), we have

[TABLE]

where

[TABLE]

The integral in the first term becomes a sum over the spectrum,

[TABLE]

For simplicity, we assume that $\kappa$ is an integer. For $\mu=0$, the sum then becomes

[TABLE]

where $n=\kappa n_{x}+n_{y}$, and $\left[x\right]$ denotes the greatest integer not exceeding $x$. Thus,

[TABLE]

where

[TABLE]

is the Hurwitz $\zeta$-function. Eq. (46) is divergent at $s\rightarrow 0$, but the divergent term is exactly canceled by another divergent term coming from $\Gamma\left(s\right)$ in the last term in eq. (43). Asymptotically expanding eq. (43) at $s\rightarrow 0$ and dropping the term proportional to $\mu$, we have

[TABLE]

where

[TABLE]

is a parameter only related to $\kappa$. Then eq. (42) becomes

[TABLE]

In this equation, both of the divergent terms of $s$ and $\mu$ are canceled, so the critical temperature can be obtained analytically by setting $N_{e}=N$. Compared with the thermodynamic-limit result, the second term in the right-hand side in eq. (50) is an extra correction. When the correction is small, the critical temperature is approximately

[TABLE]

where

[TABLE]

is the critical temperature in the thermodynamic limit. The leading term in the correction to the critical temperature is proportional to $\ln N/\sqrt{N}$, which is consistent with the leading term of the quantum correction given in ref. [23].

In fig. 5, we plot the critical temperatures eqs. (51) and (52) for different $\kappa$. It shows that our result is lower than the thermodynamic-limit value ($T_{c}<T_{0}$), and the anisotropy increases the difference between them. In this and the following figures, the temperature is rescaled to $T/T_{h}$, where $T_{h}=T_{c}\left(N=1000,\kappa=1\right)$.

3.2 Condensate fraction and chemical potential

Under the assumption $\mu=0$, the zero-order condensate fraction is easy to obtain from eq. (50),

[TABLE]

For the isotropic case, i.e. $\kappa=1$, neglecting the higher-order contribution in the third term, the zero-order condensate fraction can be expressed as

[TABLE]

In Ref. [22], the author gives an approximate result of the condensate fraction in an isotropic harmonic trap, which has the similar form as eq. (54) but the coefficient of the third term is twice as large as our result. The comparison with the numerical calculation confirms that eq. (54) is much more precise (see fig. 7).

In a finite system, the chemical potential $\mu$ is not exactly zero below the transition point. To find the analysis form of $\mu$, we need to add three terms in eq. (42) to give an equation of $N$: the ground-state particle, the next-to-leading term of the Bose-Einstein integral, and the first-order correction of $\mu$ in eq. (44). Thus eq. (42) becomes

[TABLE]

The sum of the spectrum is approximately

[TABLE]

The first term has been calculated in eq. (46), and the second term is also divergent at $s\rightarrow 0$:

[TABLE]

where

[TABLE]

is only related to $\kappa$. However, the term with $\Gamma\left(s-1\right)$ in eq. (55) is proportional to $\mu$ and should be included in this approximation. Easy to check that the divergent term coming from the gamma function and that in eq. (57) are exactly canceled. Therefore all the divergent terms of $s$ are canceled in eq. (55):

[TABLE]

By using eqs. (53) and (52), it can be rewritten as

[TABLE]

Neglecting the higher-order terms, we solve the chemical potential as

[TABLE]

where

[TABLE]

is the chemical potential at the transition point.

In fig. 6 we plot the relation between the chemical potential and temperature for different $N$. For $T<T_{c}$, eq. (61) agrees with the numerical solution quite good.

From eq. (61), the first-order condensate fraction $N_{0}^{\left(1\right)}/N$ is straightforward according to eq. (28). In fig. 7 we plot the relation between the condensate fraction and temperature for different $N$. At the critical temperature, the zero-order condensate fraction vanishes, but the first-order one matches the numerical solution very well.

4 Comparison with the canonical ensemble

In the above sections, our discussion on BEC is in the grand canonical ensemble. However, in a finite system, the fluctuation of particle number in the grand canonical ensemble may be non-negligible. For investigating the influence of fluctuation, we will consider the behavior of Bose gases in the canonical ensemble and compare the result with the grand canonical ensemble.

There are many studies on the similarities and differences between different ensembles for finite systems [14, 20, 33, 34]. In this section, we will take the two-dimensional harmonic trap as an example to show the difference between these two ensembles.

In the canonical ensemble, the partition function of a $N$-particle system is

[TABLE]

where $E_{k}^{\left(S\right)}$ is the total energy of the $k$-th system in the ensemble. However, the constraint of fixed particle number makes the exact analytical form of partition function for a quantum system hard to obtain, even for ideal gases. One method is to express the partition function by a complex integral of the grand partition function as

[TABLE]

where the integral path is a loop surrounding the original point. However, although this integral can be approximately evaluated by the saddle point method for large $N$, the exact integral can hardly be performed.

To give a direct comparison between different ensembles, we need the exact partition function. For not very large $N$, this can be achieved by use of the recursion relation [20, 33]

[TABLE]

where

[TABLE]

is the partition function for a single particle at the temperature $T/k$.

We will take the condensate fraction as an example to compare with that in the grand canonical ensemble. The average particle number in a state with energy $E_{i}$ in the canonical ensemble can be expressed as [20, 33]

[TABLE]

Combined with eq. (65), it will give the particle number in the ground state and the condensate fraction.

In fig. 8 we plot the numerical solutions of condensate fraction in the grand canonical and canonical ensembles for different $N$ in a two-dimensional harmonic trap. It is clear that for $N\sim 10^{3}$, the difference between these two ensembles is very small.

5 Conclusion and discussion

In the above, by using an analytical continuation method to solve the divergence problem in BEC, we discuss the low-temperature behavior of ideal Bose gases in the two-dimensional anisotropic box and harmonic trap, respectively. We show that the influence of boundaries and external potentials can be dealt with by a general treatment. We obtain the critical temperature, the condensate fraction and the chemical potential for Bose gases in these two kinds of traps, respectively. The results are consistent with or better than the corresponding studies in the literature, and they agree with the numerical calculation well. To check the influence of fluctuation in the canonical ensemble, we compare the condensate fraction in the grand canonical and canonical ensembles. The result shows that for about $N\sim 10^{3}$, the difference between these two ensembles is negligible.

Although some previous studies also discussed the corrections to critical temperature and condensate fraction in finite systems, our method is not an order-of-magnitude estimate, so we can obtain more precise results, including not only the leading correction. Besides, our method provides a general treatment to the problem of BEC in finite systems. As long as the heat kernel expansion is known, the critical temperature and the thermodynamic quantities of the Bose gas can be calculated.

The grand potential of a finite system also contains divergent terms at $\mu\rightarrow 0$, and this problem can also be solved by similar treatment. The analytical expressions of the grand potential and other thermodynamic quantities below the transition point can be obtained as well. However, the divergence problem is often not serious in the grand potential. For the two cases considered in this paper, the divergence appears from the third term of the grand potential. Therefore, our method will give the corrections to the third terms. Such corrections are usually negligible, so their expressions are not presented in this paper.

The advantage of our method is to remove the divergence at the transition point, so the magnitude of the correction tightly depends on the specific nature of the systems. For example, for the critical temperature, it gives the second-order correction in the case of three-dimensional harmonic traps, which is usually negligible [24]. In a two-dimensional harmonic trap, the correction is first-order and is expected to be observed in experiments. In a two-dimensional box, since no phase transition exists in the thermodynamic limit, the correction is zero-order and its influence is significant.

Recently, many experimental studies on BEC are performed in two-dimensional traps, especially the BEC of quasiparticles, such as excitons in graphene and surface exciton-polaritons. We hope that more precise experiments at this field will test our results.

The author is very indebted to Prof. Wu-Sheng Dai for his help. The author is grateful to an anonymous referee for helpful comments and suggestions which greatly improved this paper. This work is supported in part by NSF of China, under Project No. 11575125.

Appendix A: Asymptotic expansion of the Epstein $\zeta$-function

$E_{2}\left(\sigma;a_{1},a_{2}\right)$

According to ref. [35], the Epstein $\zeta$-function

[TABLE]

has a singularity $\sigma=1$, so we need the asymptotic expansion of Epstein $\zeta$-function around $\sigma=1$.

Around $\sigma=1$, only the second term in eq. (68) is divergent, which is

[TABLE]

The first and third terms in eq. (68) is convergent, so substituting $\sigma=1$ into them gives

[TABLE]

and

[TABLE]

Therefore, around $\sigma=1$, we have

[TABLE]