$1/N$-expansion for the critical temperature of the Bose gas

TL;DR

This paper applies the large-N expansion to analyze the critical temperature of a three-dimensional Bose gas with short-range repulsion, providing numerical and analytical insights into its behavior across different interaction strengths.

Contribution

It introduces a full numerical calculation of the critical temperature using the 1/N-approximation for the Bose gas and estimates its strong-coupling behavior.

Findings

Critical temperature shift grows at small gas parameter

It reaches a maximum and then decreases, becoming negative

Provides analytical estimates for strong-coupling limit

Abstract

We revised the large-$N$ expansion for a three-dimensional Bose system with short-range repulsion in normal phase. Particularly, for the model potential that is characterised only by the $s$-wave scattering length $a$ the full numerical calculations of the critical temperature in the $1/N$-approximation as a function of the gas parameter $an^{1/3}$ are performed. Additionally to the well-known result in the dilute limit we estimated analytically the leading-order strong-coupling behavior of the Bose-Einstein condensation transition temperature. It is shown that the critical temperature shift of the non-ideal Bose gas grows at small $an^{1/3}$, reaches some maximal value and then falls down becoming negative.

Equations

$$\frac{T_c-T_0}{T_0}=ca{n}^{1/3}.$$

$$S=\int dx{\psi }_{\sigma }^{\ast }(x)\left\{{\partial }_{\tau }+{\hslash }^2{\nabla }^2/2m+\mu \right\}{\psi }_{\sigma }(x)$$

$$-\frac{1}{2N}\int dx\int d{x}^{\prime }\Phi (x-x^{\prime })|{\psi }_{\sigma }(x){|}^2|{\psi }_{{\sigma }^{\prime }}(x^{\prime }){|}^2,$$

$$S=\int dx{\psi }_{\sigma }^{\ast }(x)\left\{{\partial }_{\tau }+{\hslash }^2{\nabla }^2/2m+\mu -i\phi (x)\right\}{\psi }_{\sigma }(x)$$

$$-\frac{N}{2}\int dx\int d{x}^{\prime }{\Phi }^{-1}(x-x^{\prime })\phi (x)\phi (x^{\prime }),$$

$$S=\sum _P\left\{i{\omega }_p-{\epsilon }_p+\tilde{\mu }\right\}{\psi }_{\sigma P}^{\ast }{\psi }_{\sigma P}+\frac{N}{2}\beta V{n}^2{\nu }_0$$

$$-\frac{N}{2}\sum _K{\nu }_k^{-1}|{\phi }_K{|}^2-\frac{i}{\sqrt{\beta V}}\sum _{K,P}{\phi }_K{\psi }_{\sigma P}^{\ast }{\psi }_{\sigma P-K},$$

$$G_{\sigma }^{-1}(P)=i{\omega }_p-{\epsilon }_p+\tilde{\mu }-{\Sigma }_{\sigma }(P),$$

$$⟨|{\phi }_K{|}^2{⟩}^{-1}=N\left\{{\nu }_k^{-1}+\Pi (K)\right\},$$

$$\Sigma (P)=\frac{-1}{\beta V}\sum _K\Gamma (P+K,P)⟨|{\phi }_K{|}^2⟩G(P+K),$$

$$\Pi (K)=\frac{1}{\beta V}\sum _P\Gamma (P+K,P)G(P)G(P+K).$$

$$⟨|{\rho }_K{|}^2{⟩}^{-1}={\nu }_k+{\Pi }^{-1}(K),$$

$$\frac{\partial G^{-1}(P)}{\partial \tilde{\mu }}=\underset{K\to 0}{\lim }\Gamma (P+K,P),$$

$$\underset{K\to 0}{\lim }\Pi (K)=\frac{\partial n}{\partial \tilde{\mu }}.$$

$${\Sigma }^{(1)}(P)=-\frac{1}{\beta V}\sum _K\frac{{\nu }_k}{1+{\nu }_k\Pi (K)}G(P+K),$$

$$n=-\underset{\tau \to +0}{\lim }\frac{1}{\beta V}\sum _P{e}^{i{\omega }_p\tau }G(P).$$

$$n=n_0+\frac{1}{NV}\sum _{\mathbf{p}}\frac{\partial n(\beta {\epsilon }_p)}{\partial {\epsilon }_p}\Delta {\Sigma }_R^{(1)}({\tilde{\xi }}_p,p)+\dots ,$$

$$\frac{T_c-T_0}{T_0}=-\frac{2}{3}\frac{n-n_0}{n_0},$$

$$\frac{T_c-T_0}{T_0}=\frac{1}{N}{f}^{(1)}(a{n}^{1/3})+o(1/N),$$

$$f^{(1)}(a{n}^{1/3})=\frac{8\pi }{3[\zeta (3/2){]}^{4/3}}a{n}^{1/3}+\dots .$$

$$f^{(1)}(a{n}^{1/3})=-\frac{64}{45}\sqrt{\frac{a^{3}n}{\pi }}+\dots .$$

$${\Sigma }^{(1)}(P)$$

$$\frac{1}{V}\sum _{\mathbf{k}\ne 0}{\int }_{-\infty }^{\infty }\frac{d\omega }{\pi }\frac{{\nu }_k^{2}I(\omega ,k)}{{\left|1+{\nu }_k\Pi (\omega ,k)\right|}^2}$$

$$\times \frac{n(\beta \omega )-n(\beta {\tilde{\xi }}_{|\mathbf{k}+\mathbf{p}|})}{\omega -{\tilde{\xi }}_{|\mathbf{k}+\mathbf{p}|}+i{\omega }_p}.$$

$$\Delta {\Sigma }_R^{(1)}({\tilde{\xi }}_p,p)=$$

$$-\frac{4{\epsilon }_p}{3V}\sum _{\mathbf{k}\ne 0}\frac{n{\nu }_k^2{\epsilon }_k^2}{E_k(E_k+{\epsilon }_k{)}^3}+o(\sqrt{a^{3}n}).$$

$$I(\omega ,k)=\frac{\beta p_0^4}{16\pi k}\ln \left|\frac{1-e^{-[k/2p_0+\beta \omega /(2k/p_0){]}^2}}{1-e^{-[k/2p_0-\beta \omega /(2k/p_0){]}^2}}\right|.$$

$$R(\omega ,k)=\text{v.p.}{\int }_{-\infty }^{\infty }\frac{d{\omega }^{\prime }}{\pi }\frac{I({\omega }^{\prime },k)}{{\omega }^{\prime }-\omega }.$$

Full text

Bose system, critical temperature, $1/N$-expansion

I Introduction

Two decades ago it was clarified Gruter_Ceperley_Laloe that the transition temperature of a three-dimensional homogeneous weakly interacting system of hard-sphere bosons increases linearly with the increase of the $s$-wave scattering length

[TABLE]

Soon after its publication this result was confirmed Holzmann_Krauth except for the coefficient $c$. The value of this coefficient strongly depends on the specific calculation scheme and was lively debated Baym_etal_99 ; Baym_etal_01 ; Cruz_Pinto_Ramos ; Kneur_Pinto_Ramos ; Braaten_Radescu ; Wang ; Ledowski_Hasselmann_Kopietz ; Blaizot ; Morawetz_Mannel_Schreiber . More precise simulations with classical $\phi^{4}$ model performed in Refs. Kashurnikov_Prokofev_Svistunov ; Arnold_Moore are consistent to each other and give for $c\simeq 1.3$. The same value for the interaction-induced shift of the phase-transition temperature was also obtained in the later path-integral Monte Carlo simulations Nho_Landau . This discrepancy of different numerical approaches is discussed in detail in Ref. Pilati_Giorgini_Prokofev . In the leading order of the $1/N$-expansion, which effectively sums up particle-hole diagrams, the coefficient that determines the critical temperature shift is found to be $c=2.33$ Baym_Blaizot_Zinn-Justin . On the other hand, the summation of the particle-particle ladder diagrams Stoof gives a twice larger value of this coefficient. The calculations in the so-called fluctuation–exchange approximation Tsutsui_Kita which incorporates both particle-particle and particle-hole bubbles lead to $c=2.94$. The next-to-leading order in the $1/N$ term of $c$ lowers the result by $26\%$ (for a single component Bose system) Arnold_Tomasik which particularly signals a good convergence of the large-$N$ expansion. These calculations however are complicated by the fact that the first correction beyond the positive linear shift in $an^{1/3}$ to the critical temperature is nonanalytic function Holzmann_etal ; Arnold_Moore_Tomasik ; Cruz_etal and more a accurate value of $c$ can be obtained only through the five- Kleinert , and seven-loop Kastening calculations. An alternative way to improve the results for the transition temperature of a dilute Bose gas is to use various resummation procedures Yukalovs where very good agreement with Monte Carlo simulations is reached.

The goal of the present study is to explore the leading-order large-$N$ correction to the Bose-Einstein condensation temperature in the wide region of the interparticle interaction strength and compare the obtained results with recent Monte Carlo calculations Bronin ; Nguyen .

II Formulation

The imaginary-time action of the considered model reads

[TABLE]

where $x\equiv(\tau,{\bf r})$, $\int dx=\int_{0}^{\beta}d\tau\int_{V}d{\bf r}$, and the summation over repeated sort indices $\sigma,\sigma^{\prime}=1,2,\ldots,N$ is understood. As usual Popov , the path-integral is carried out over complex $\beta$-periodic ($1/\beta=T$ is the temperature) fields $\psi_{\sigma}(x)$. We also introduced the chemical potential $\mu$ that fixes the total equilibrium density $Nn$ of the system. The second term of action (II) describes interaction between particles with potential $\Phi(x)=\delta(\tau)\Phi({\bf r})$. For later consideration the explicit form of the two-body potential $\Phi({\bf r})$ is not important except it must have the Fourier transform defined. Introducing new real field $\varphi(x)$ and making use of the Hubbard-Stratonovich transformation we rewrite action $S$ in the equivalent form (up to an irrelevant constant term)

[TABLE]

where the inverse operator $\Phi^{-1}(x-x^{\prime})$ satisfies natural identity $\int dx^{\prime\prime}\Phi^{-1}(x-x^{\prime\prime})\Phi(x^{\prime\prime}-x^{\prime})=\delta(x-x^{\prime})$. The main advantage of such a decomposition lies in the rearranging of the perturbation theory in order to stress the leading role of particle-hole diagrams in the large-$N$ limit. Separating the uniform classical part of field $\varphi(x)=\varphi_{0}+\tilde{\varphi}(x)$ (with constraint $\int d{\bf r}\tilde{\varphi}(x)=0$) and using the steepest descend method for the grand canonical potential $\partial\Omega/\partial\varphi_{0}=0$ we obtain $i\varphi_{0}=n\int dx\Phi(x)$, where density $n$ of each sort of particles should be treated as a function of chemical potential. Finally, action (II) in momentum space is

[TABLE]

where capital letters denote four-momenta $P\equiv(\omega_{p},{\bf p})$, $K\equiv(\omega_{k},{\bf k}\neq 0)$ (here $\omega_{k},\omega_{p}$ are bosonic Matsubara frequencies) and $\varepsilon_{p}=\hbar^{2}p^{2}/2m$ is free-particle dispersion. We also used notations for the shifted chemical potential $\tilde{\mu}=\mu-n\nu_{0}$ and $\nu_{k}$ for the Fourier transform of the spherically symmetric inter-particle potential $\Phi(r)$. Dyson equations that determine both one-particle Green’s function $G_{\sigma}(P)=-\langle\psi_{\sigma P}\psi^{*}_{\sigma P}\rangle$ and correlator $\langle|\varphi_{K}|^{2}\rangle$ read

[TABLE]

[TABLE]

where self-energy $\Sigma_{\sigma}(P)$ and polarization operator $\Pi(K)$ are uniquely determined by vertex function $\Gamma_{\sigma}(P+K,P)$ (in the following we omit subscript $\sigma$ everywhere)

[TABLE]

[TABLE]

This formulation clearly simplifies the perturbative ($1/N$) expansion classifying diagrams by the number of $\varphi$-correlators.

The obtained effective theory constitutes Bose particles coupled to the field $\varphi_{K}$ that describes their own collective excitations. To reveal this connection explicitly it is enough to calculate the two-point density fluctuations function for each component of the Bose system

[TABLE]

which is related to the experimentally measured dynamical structure factor. Poles of $\langle|\rho_{K}|^{2}\rangle$ after analytical continuation determine the spectrum of collective modes.

The definition of the shifted chemical potential allows to rewrite the thermodynamic relation for the equilibrium density as follows $-\partial\tilde{\Omega}/\partial\tilde{\mu}=Nn$, where $\tilde{\Omega}=\Omega+\frac{N}{2}Vn^{2}\nu_{0}$ depends on temperature and $\tilde{\mu}$ only. It also naturally enters the Ward identities in the so-called static limit ($\omega_{k}=0,{\bf k}\rightarrow 0$)

[TABLE]

[TABLE]

Note that combining the above second identity (2.11) and the definition of $\tilde{\mu}$ we find a full agreement of the low-energy limit of the density-density correlation function (2.9) $\lim_{K\rightarrow 0}\langle|\rho_{K}|^{2}\rangle=\partial n/\partial\mu$ with the compressibility sum rule.

III Leading-order results

III.1 $1/N$-approximation

We restrict ourselves considering only the simplest approximation of order $1/N$. In this limit one neglects corrections to the vertex function, i.e., $\Gamma(P+K,P)=1$ and substitutes bare Green’s function $G^{-1}(P)=i\omega_{p}-\tilde{\xi}_{p}$, ($\tilde{\xi}_{p}=\varepsilon_{p}-\tilde{\mu}$) in Eqs. (2.5), (2.6) for the self-energy and for the polarization operator. Explicitly writing down the dependence of the self-energy $\Sigma(P)=\Sigma^{(1)}(P)/N+\ldots$ on $N$ in the leading order we find

[TABLE]

where $\Pi(K)=\frac{1}{\beta V}\sum_{P}G(P)G(P+K)$ is the polarization operator in the adopted approximation. In order to define the summation over $\omega_{k}$ correctly one should carefully single out the Hartree-Fock (actually, only the Fock) contribution while calculating $\Sigma^{(1)}(P)$ (see appendix A for details). After this the derivation of the renormalized one-particle spectrum becomes a typical routine connected with the analytical continuation $\Sigma(P)|_{i\omega_{p}=\omega+i0}=\Sigma_{R}(\omega,p)+i\Sigma_{I}(\omega,p)$ of the self-energy and finding poles of Green’s function. Up to the first order over the expansion parameter $1/N$ the impact of the interaction effects on the quasiparticle dispersion is given by $\xi^{*}_{p}=\tilde{\xi}_{p}+\Sigma^{(1)}_{R}(\tilde{\xi}_{p},p)/N$.

We perform the calculation of the Bose-Einstein condensation transition temperature $T_{c}$ in a standard way

[TABLE]

Obviously, in the transition point $\xi^{*}_{p}$ tends to zero in the long-length limit, therefore one has to take into account such a non-perturbative peculiarity of these calculations. Making use of spectral representation (via imaginary part $\Sigma_{I}(\omega,p)$ of the retarded function) for the self-energy and after straightforward summation over Matsubara frequencies we find in the accepted approximation

[TABLE]

where for convenience we denoted the shifted self-energy $\Delta\Sigma^{(1)}_{R}(\tilde{\xi}_{p},p)=\Sigma^{(1)}_{R}(\tilde{\xi}_{p},p)-\Sigma^{(1)}_{R}(-\tilde{\mu},0)$, introduced notations for the Bose distribution function $n(y)=(e^{y}-1)^{-1}$ and for density $n_{0}=(p_{0}/2\sqrt{\pi})^{3}\zeta(3/2)$ (here and after $p_{0}=\sqrt{2mT}/\hbar$, and $\zeta(s)=\sum_{n\geq 1}1/n^{s}$ is the Riemann zeta function) of the ideal gas at the critical temperature. Finally, we conclude this section with the equation that relates the shift of the critical temperature with the density jump Baym_etal_01 in three dimensions

[TABLE]

where $T_{0}$ is the Bose condensation temperature of non-interacting system.

Till now we were discussing the properties of Bose systems not specifying the form of the inter-particle interaction. The only requirement to potential was the existence of a well-defined Fourier transform. In the following sections the calculations of critical temperature will be performed for a model with a short-range repulsion, i.e., $\nu_{k}=4\pi\hbar^{2}a/m$ which is characterized only by one parameter, namely the $s$-wave scattering length $a$.

III.2 Analytical limits

We start with the discussion of the dilute limit which is well-studied due to the long history of investigations. In the present approach (see appendix A) this limit is reproduced in the so-called classical approximation, when all Bose distribution functions in the second term of Eq. (3.14) are replaced by their infrared asymptotes $n(y)\rightarrow 1/y$. After these simplifications, the function that determines the shift of the critical temperature

[TABLE]

can be calculated analytically in the limit of $an^{1/3}\rightarrow 0$ Baym_Blaizot_Zinn-Justin

[TABLE]

This linear dependence on $an^{1/3}$ appears only because the large-$N$ expansion even in the simplest approximation sums up infinite series of diagrams divergent near the phase transition point. Any finite order of the conventional perturbation theory necessarily leads to the incorrect behavior of the critical temperature shift (for instance, inclusion of the second-order terms only results in $f^{(1)}(an^{1/3})\propto\sqrt{an^{1/3}}$ Huang ; Bala_Srivastava_Pathak ).

Another highly non-trivial case, where the leading-order behavior of the function $f^{(1)}(an^{1/3})$ can be obtained analytically, is the strong-coupling limit $an^{1/3}\gg 1$ Pastukhov . To observe this possibility it is enough to analyze the integral over $\omega$ in Eq. (V.1). The first multiplier in the integrand is proportional to the imaginary part of the retarded density-density correlation function (2.9). To reveal the behavior of this correlator we calculated the polarization operator $\Pi(K)$ (see appendix B for details) and plot the sketch of the dynamical structure factor $S(\omega,k)=\frac{1}{\pi}(1-e^{-\beta\omega})^{-1}\Im\langle|\rho_{K}|^{2}\rangle|_{i\omega_{k}\rightarrow\omega+i0}$ of this model at fixed coupling constant and for some values of the wave-vector in Fig. 1.

Particularly, it is seen that in the strong-coupling limit the system possesses a well-defined spectrum of collective modes. It also easy to show that the damping rate in the phonon region of spectrum decreases exponentially as the gas parameter increases. Up to the leading order in $an^{1/3}$ this spectrum of collective excitations can be approximated by the Bogoliubov one. This excitingly simple fact allows to provide further calculations analytically and for the function $f^{(1)}(an^{1/3})$ we obtain (see appendix A) at large arguments

[TABLE]

First note that unlike to the dilute limit (3.17) the critical temperature in this case decreases. The second less noticeable fact is that this leading-order result in the limit $an^{1/3}\gg 1$ is fully determined by the effects of particle mass renormalization. On the level of our $1/N$-approximation this effective mass of a moving particle is nothing but the hydrodynamic mass of a single impurity atom immersed in the “phonon” field $\varphi(x)$. This situation with the function $f^{(1)}(an^{1/3})$ when it grows at small $an^{1/3}$ and rapidly falls down to negative values in the strong-coupling limit is consistent with experiments Reppy and with the results of simulations Bronin ; Nguyen and has a simple physical interpretation. For our model the enhancement of the interaction can be naively treated as an increase of the particle size. In the ultra-dilute Bose gas where the two-particle collision processes are less probable the presence of the short-range repulsion only leads to the reduction of the free volume of the system, that is why the Bose condensation temperature slightly increases in this limit. When interaction grows the Feynman mechanism Feynman also takes effect and leads to the effective increasing of a particle inertial mass. As a consequence, the competition of these two physical mechanisms of the critical temperature formation provides that $f^{(1)}(an^{1/3})$ is a non-monotonic function of $an^{1/3}$.

III.3 Numerical calculations

Although formula (3.14) with the self-energy given by (5.19) and (V.1) is applicative for the analytical calculations in diverse limits, the full numerical computations require another approach. From the practical point of view it is more convenient to perform primarily integration over wave-vector in the self-energy and after that to calculate the sum over the Matsubara frequency. Additionally, the zero-frequency term should be singled out and treated very carefully Baym_Blaizot_Zinn-Justin . The numerical estimation of the function $f^{(1)}(an^{1/3})$, which is the main result of the present paper, is depicted in Fig. 2.

The maximum of the critical temperature due to our $1/N$-calculations is reached at $an^{1/3}=0.125$ and after $an^{1/3}=0.345$ the shift of $T_{c}$ of the Bose system with the fully repulsive interparticle interaction is negative. For comparison we also provide in Fig. 3 the results of recent path-integral Bronin and quantum Nguyen Monte Carlo simulations.

IV Conclusions

In summary, by using $1/N$-expansion to the leading order we performed the full numerical calculations of the critical temperature for Bose system with short-range repulsion. Particularly it is shown by the direct computer computations as well as by the analytical estimations in various limits that the Bose-Einstein condensation transition temperature of this model is a non-monotonic function of the gas parameter. We have argued that for a repulsive inter-particle interaction with the large $s$-wave scattering length or for a dense Bose gas the critical temperature shift is negative even in the simplest $1/N$-approximation. Besides, it is demonstrated that the estimation for the transition temperature of hard-sphere bosons obtained by means of large-$N$ expansion coincides comparatively well with the essentially exact results of numerical methods.

Acknowledgements

Insightful discussions with Prof. Ivan Vakarchuk and Dr. Andrij Rovenchak are gratefully acknowledged. We thank Prof. Matthias Troyer and Dr. Sebastiano Pilati for providing us with the results of their quantum Monte Carlo simulations. This work was partly supported by Project FF-30F (No. 0116U001539) from the Ministry of Education and Science of Ukraine.

V Appendices

V.1 Self-energy

The summation over the Matsubara frequency $\omega_{k}$ in self-energy (3.12) is not well-defined. To deal with this problem we single out the leading-order ultraviolet term of the fraction and rewrite $\Sigma^{(1)}(P)$ in the equivalent form

[TABLE]

where the first term is the Fock correction. Now the second sum is convergent and by using of the spectral theorems can be written in terms of real and imaginary part of the retarded polarization operator $\Pi(K)|_{i\omega_{k}\rightarrow\omega+i0}=\Pi(\omega,k)=R(\omega,k)+iI(\omega,k)$

[TABLE]

In the strong-coupling limit, where the dynamical structure factor of the system has a sharp peak which corresponds to the quasiparticles with the Bogoliubov-like dispersion $E^{2}_{k}=\varepsilon^{2}_{k}+2n\nu_{k}\varepsilon_{k}$, one can integrate over $\omega$ in Eq. (V.1) and obtain in the leading order:

[TABLE]

V.2 Polarization operator

In order to reveal the structure of collective modes of the system and to calculate the self-energy we give some details of the polarization operator calculations. The imaginary part of the retarded polarization operator calculated at $T_{c}$ reads in the adopted approximation

[TABLE]

For the real part $R(\omega,k)$ we use well-known spectral theorem

[TABLE]

The spectrum of collective modes is determined by the equation $1+\nu_{k}R(E_{k},k)=0$. Using Eqs. (5.22), (5.23) it is easy to obtain the asymptotic formula for $R(\omega,k)$ Hore_Frankel and to argue that $E_{k}$ tends to the Bogoliubov spectrum when $an^{1/3}\gg 1$.