Non-Saturation in Interacting Bose Gas Integral

Question

I am independently working through some problems on Bose-Einstein condensation. In particular, I am trying to show that—in the Hartree-Fock mean-field approximation—for a Bose gas with contact interactions, the first order non-saturation effect is given by:

$$ \frac{N'}{N_c^0} = 1 + \alpha \frac{\mu_0}{k_bT} $$

where $N_c^0 = \zeta(3)\left(\frac{k_bT}{\hbar \bar{\omega}} \right)^3$, with $\zeta$ the Riemann zeta function and $\bar{\omega} = (\omega_x\omega_y\omega_z)^{1/3}$ where $\omega_i$ are the trap frequencies, and $\alpha = \zeta(2)/\zeta(3)$.

From what I understand, this result is obtained by evaluating the integral

$$\int \frac{g_{3/2}(e^{(\mu - V(\textbf{r})/k_bT})}{\lambda^3} d\textbf{r}$$

with the effective interaction potential $V_\text{eff}(\textbf{r}) = |V(\textbf{r}) - \mu_0| + \mu_0$, where $\mu_0 = gn_0(\textbf{r})$ is the chemical potential due to contact interactions and $V(\textbf{r})$ is the external trapping potential $\sum (1/2) m\omega_i^2r_i^2$.

While I am able to obtain $N_c^0 = \zeta(3)\left(\frac{k_bT}{\hbar \bar{\omega}} \right)^3$ by setting $\mu = 0$ in the same integral, I am unable to reproduce the non-saturation result. Even when I use the relation,

$$\int_0^x d\epsilon \sqrt{\epsilon}e^{\epsilon-x} + \int_x^\infty d\epsilon\sqrt{\epsilon}e^{x-\epsilon} \approx \frac{\sqrt{2}}{2}(1+x)$$

for $x\ll1$, I get stuck and am unable to reproduce to analytically evaluate the integral.

Any guidance would be appreciated.

Answer 1

I don’t know what solution the author of the problem had in mind, but the very first formula is valid for an ideal Bose gas.

The local concentration of an ideal Bose gas in an external potential in three dimensions is equal to $$ n(\vec{r}) = \frac1{(2\pi\hbar)^3}\int d^3\vec{p}\ \frac1{e^{\beta\frac{\vec{p}^2}{2m}-\beta\mu(\vec{r})}-1},\quad \beta = \frac1{k_bT}=\frac1\theta. $$ The gas equilibrium condition in an external potential has the form $\mu(\vec{r}) + V(\vec{r}) = \mu(0) + V(0) = \mu_0$. Therefore the number of particles of Bose gas is equal to $$ N(\mu_0,T) = \int d^3\vec{r} n(\vec{r}) = \int\!\!\int d^3\vec{r} d^3\vec{p}\ \frac1{e^{\beta\frac{\vec{p}^2}{2m}+\beta V(\vec{r})-\beta\mu_0}-1} $$ Here $V(\vec{r}) = \frac{m}2\sum_{i=1}^3\omega_i^2x_i^2$, therefore the following change of variables in the last integral $$ p_i = \sqrt{2m\theta}\ \xi_i, \quad x_i = \sqrt{\frac{2\theta}{m\omega_i^2}}\ \xi_{i+3},\quad i = 1,2,3 $$ gives $$ N(\mu_0,T) = \left(\frac{\theta}{\pi\hbar\omega_0}\right)^3\int d^6\vec{\xi}\ \frac1{e^{\ \vec{\xi}^2-\beta\mu_0}-1},\quad \omega_0 = (\omega_1\omega_2\omega_3)^\frac13 $$ The integrand in the last integral has symmetry with respect to rotations of the vector $\vec{\xi}$, so we further have $$ N(\mu_0,T) = \left(\frac{\theta}{\pi\hbar\omega_0}\right)^3 S_6\int_0^\infty d\xi\ \xi^5\frac1{e^{\ \xi^2-\beta\mu_0}-1},\quad S_6 = \pi^3. $$ Next, change of the integration variable $\xi = \sqrt{t}$ and representing the fraction as the sum of a geometric progression leads to $$ N(\mu_0,T) = \left(\frac{\theta}{\hbar\omega_0}\right)^3\frac12\int_0^\infty dt\ t^2\frac1{e^{t-\beta\mu_0}-1} = \left(\frac{\theta}{\hbar\omega_0}\right)^3\frac12\sum_{n=1}^\infty \int_0^\infty dt\ t^2e^{-tn}e^{\beta\mu_0 n} = $$ $$ = \left(\frac{\theta}{\hbar\omega_0}\right)^3 \sum_{n=1}^\infty \frac1{n^3} e^{\beta\mu_0 n}\tag{1} $$ In the case when $\mu_0 = 0$, the last equality implies $$ N_c = N(0, T) = \left(\frac{\theta}{\hbar\omega_0}\right)^3\sum_{n=1}^\infty \frac1{n^3} = \left(\frac{k_b T}{\hbar\omega_0}\right)^3\zeta(3). $$ In the case when $0< -\beta\mu_0 \ll 1$, the expansion of the exponential in (1) leads to $$ N' = \left(\frac{\theta}{\hbar\omega_0}\right)^3\left(\sum_{n=1}\frac1{n^3} + \beta\mu_0\sum_{n=1}^\infty\frac1{n^2}+\ldots\right)\approx N_c\left(1 + \frac{\mu_0}{k_bT}\frac{\zeta(2)}{\zeta(3)}\right), $$ and this is exactly the formula that needs to be derived.

If using the equilibrium condition $\mu(\vec{r}) + V(\vec{r}) = \mu(0) + V(0) = \mu_0$ is undesirable for some reason, I am sure that the same result can be obtained by considering the energy levels of a quantum oscillator and approximating the sums by integrals.