Identities and Many-Body Approaches in Bose-Einstein Condensates

TL;DR

This paper reviews exact relations in Bose-Einstein condensates and discusses many-body approximation methods that satisfy these relations, bridging theoretical exact identities and practical approximation approaches.

Contribution

It introduces key exact relations in BECs and discusses approximation methods that preserve these relations, connecting fundamental theory with practical modeling.

Findings

Identifies key exact relations like Hugenholtz--Pines and Nepomnyashchii identities.

Reviews many-body approximation methods that satisfy these exact relations.

Serves as a bridge between theoretical identities and approximation techniques.

Abstract

This paper discusses exact relations in Bose--Einstein condensates (BECs), starting from basic properties of an ideal Bose gas. In particular, focused on are the Hugenholtz--Pines relation, Nepomnyashchii--Nepomnyashchii identity, and identities for the density response function. After introducing these exact relations, a few approaches of many-body approximations are discussed, which satisfy the exact relations in BECs. This paper will serve as a bridge between theories on exact relations and those on approximations in BECs.

Tables

Table 1. Table 1: Critical temperature shift Δ ​ T c = c 1 ​ a ​ n 1 / 3 Δ subscript 𝑇 c subscript 𝑐 1 𝑎 superscript 𝑛 1 3 \Delta T_{\rm c}=c_{1}an^{1/3} , evaluated by numerical calculations: Monte-Carlo simulation  Arnold2001 ; Kashrnikov2001 , random-phase approximation (RPA), simplified RPA (s-RPA), and many-body T 𝑇 T -matrix theory (MBT)  Watabe2013 .

	Monte-Carlo Arnold2001 ; Kashrnikov2001	RPA Watabe2013	s-RPA Watabe2013	MBT  Watabe2013
$c_1$	$1.3$	$1.1$	$2.1$	$3.9$

Equations

$$\zeta \equiv \frac{1}{h^3}\int d\mathbf{p}\int d\mathbf{r}\exp \left(-\beta {ϵ}_{\mathbf{p}}\right)=\frac{V}{{\lambda }_T^3}.$$

$$S=k_{\mathrm{B}}N\left[\frac{5}{2}+\ln {\left(\frac{d}{{\lambda }_T}\right)}^3\right],$$

$$\Psi ({\mathbf{r}}_1,{\mathbf{r}}_2;t)=\pm \Psi ({\mathbf{r}}_2,{\mathbf{r}}_1;t),$$

$$f({ϵ}_i)=\frac{1}{\exp [\beta ({ϵ}_i-\mu )]-1},$$

$$\rho (ϵ)\equiv \sum _{\mathbf{p}}\delta (ϵ-{ϵ}_{\mathbf{p}})=\frac{V{m}^{3/2}}{\sqrt{2}{\pi }^2{\hslash }^3}\sqrt{ϵ},$$

$$N=\sum _{\mathbf{p}}f({ϵ}_{\mathbf{p}})={\int }_0^{\infty }dϵ\rho (ϵ)f(ϵ)=\frac{V}{{\lambda }_T^3}{\mathcal{G}}_{3/2}(z),$$

$${\mathcal{G}}_{\alpha }(z)\equiv \frac{1}{\Gamma (\alpha )}{\int }_0^{\infty }dx\frac{x^{\alpha -1}}{\exp (x)z^{-1}-1}=\sum _{n=1}^{\infty }\frac{z^n}{n^{\alpha }}$$

$${\rho }_{\mathrm{sd}}=\frac{N}{V}{\lambda }_T^3={\mathcal{G}}_{3/2}(z\to 1)=\zeta (3/2).$$

$$N=f({ϵ}_{\mathbf{p}=0})+\sum _{\mathbf{p}\ne 0}f({ϵ}_{\mathbf{p}})\equiv N_0+N^{\prime }.$$

$$N^{\prime }=$$

$$\mu =-\frac{k_{\mathrm{B}}T}{n_{0}V},$$

$$k_{\mathrm{B}}{T}_{\mathrm{c}}^0\equiv \frac{2\pi {\hslash }^2}{m}{\left[\frac{N}{V}\frac{1}{\zeta (3/2)}\right]}^{3/2}.$$

$${\hat{H}}_0=$$

$$G_0(\mathbf{r},\tau ;{\mathbf{r}}^{\prime },{\tau }^{\prime })\equiv$$

$$\equiv$$

$$\hat{\Psi }(\mathbf{r},\tau )=\frac{1}{\sqrt{V}}\sum _{\mathbf{p}}\exp \left(\frac{i\mathbf{p}\cdot \mathbf{r}}{\hslash }\right){\hat{a}}_{\mathbf{p}}(\tau ).$$

$$G_0(\mathbf{p},i{\omega }_n)=$$

$$=$$

$$f({ϵ}_{\mathbf{p}})=$$

$$=$$

$$\hat{\Psi }(\mathbf{r})=$$

$$\to$$

$$\left[\frac{{\hat{a}}_{\mathbf{p}=0}}{\sqrt{V}},\frac{{\hat{a}}_{\mathbf{p}=0}^{†}}{\sqrt{V}}\right]=\frac{1}{V},$$

$$\underset{|{\mathbf{r}}_1-{\mathbf{r}}_2|\to \infty }{\lim }{\rho }_{(1)}({\mathbf{r}}_1,{\mathbf{r}}_2)=\frac{N_0}{V}.$$

$$\hat{H}={\hat{H}}_0+\frac{U}{2}\int d\mathbf{r}{\hat{\Psi }}^{†}(\mathbf{r}){\hat{\Psi }}^{†}(\mathbf{r})\hat{\Psi }(\mathbf{r})\hat{\Psi }(\mathbf{r}).$$

$$\frac{4\pi a}{m}=$$

$$\mathbf{G}(p)={\mathbf{G}}_0(p)+{\mathbf{G}}_0(p)\mathbf{\Sigma }(p)\mathbf{G}(p),$$

$$\mathbf{G}(p)=\left(\begin{array}{cc}{G}_{11}(p)& G_{12}(p)\\ G_{21}(p)& G_{22}(p)\end{array}\right),\mathbf{\Sigma }(p)=\left(\begin{array}{cc}{\Sigma }_{11}(p)& {\Sigma }_{12}(p)\\ {\Sigma }_{21}(p)& {\Sigma }_{22}(p)\end{array}\right),$$

$$G_{11}(p)=$$

$$G_{12}(p)=$$

$$D(p)\equiv$$

$$\times \left[{ϵ}_{\mathbf{p}}-\mu +\hslash S(p)-\hslash {\Sigma }_{12}(p)\right],$$

$$A(p)\equiv$$

Full text

I Introduction

Bose–Einstein condensation is a dramatic phenomenon, where quantum mechanics emerge in a macroscopic scale. Bose–Einstein condensates (BECs) provide rich variety of topics in condensed matter physics, such as superfluidity, phase transition, symmetry breaking and so on AndersonBook ; LeggettBook ; PitaevskiiBook . These topics are strongly related to other condensed matter physics such as superconductors as well as magnetisms AndersonBook ; LeggettBook . The BEC is thus a cornucopia of basic concepts of condensed matter physics.

During a long history of the study of BECs, beautiful identities have been shown Hugenholtz1959 ; Gavoret1964 ; Nepomnyashchii1975 ; Nepomnyashchii1978 ; GriffinBook . These are related to the gapless excitation Hugenholtz1959 , the infrared divergence of the longitudinal susceptibility Nepomnyashchii1975 ; Nepomnyashchii1978 , which is strongly related to the magnetism Vaks1967 ; Patashinskii1973 , and the correspondence between sound speed of single-particle excitation and that of the collective density excitation in a BEC Gavoret1964 . On the other hand, from a view point of a more practical approximation side, not only the difficult problem still remains for the dilemma of the conserving-gapless approximation  Hohenberg1965 , but the treatment of the infrared divergence is also difficult Shi1998 . We furthermore still often encounter gaps between known exact properties and approximation frameworks that do not satisfy them.

This paper may serve as a bridge that connects these gaps, by starting to introduce basic properties of an ideal Bose gas. We then introduce exact relations for an interacting condensed Bose gas, such as the Hugenholtz–Pines relation Hugenholtz1959 , Nepomnyashchii–Nepomnyashchii identity Nepomnyashchii1975 ; Nepomnyashchii1978 , as well as the exact relations with respect to the density response function Gavoret1964 . Finally, we discuss many-body approaches satisfying these identities Watabe2013 ; Watabe2014 ; Watabe2018 . Development of approximation frameworks for BEC is not completed, in the sense that there is not a practical approximation satisfying all the exact relations in BECs. We hope that this paper helps to guide readers towards theories on BECs, and is useful for developing them.

II Quantum Statistical Mechanics of Bosons

An ideal classical gas is well described by the canonical ensemble. At very low temperature, however, the treatment of a classical gas is no longer valid. Its temperature scale is described by using one of the typical length scales of a nonzero temperature system — the thermal de Broglie length $\lambda_{T}\equiv[2\pi\hbar^{2}/(mk_{\rm B}T)]^{1/2}$, where $T$ is the temperature, and $m$ is an atomic mass of a gas LeggettBook . When this thermal de Broglie length has a same or longer length scale of the mean interatomic distance $d$, i.e., $\lambda_{T}\gtrsim d$, the treatment of the classical gas becomes failure, and the quantum statistical mechanics is needed.

Indeed, the partition function $Z$ of a three-dimensional ideal classical gas is given by $Z=\zeta^{N}/N!$ with the Boltzmann counting factor $N!$, where

[TABLE]

Here, $\beta$ is the inverse temperature $\beta=1/(k_{\rm B}T)$, $\epsilon_{\bf p}\equiv{\bf p}^{2}/(2m)$ is a kinetic energy of a particle in a uniform system, $V$ is the volume of the system. The entropy $S=-\partial F/\partial T$ with the Helmholtz free energy $F=-\beta^{-1}\ln Z$ is given by

[TABLE]

where the mean interatomic distance $d\equiv(V/N)^{1/3}$ GreinerBook . Now that we may find that an unphysical negative entropy is shown in the case $\lambda_{T}\gtrsim d$, we should not consider the classical gas treatment in this region. Instead, we consider the quantum statistical mechanics for identical particles at very low temperature $k_{\rm B}T\lesssim 2\pi\hbar^{2}/(md^{2})$.

This statistics may be found by the symmetry of the multiparticle wavefunction. The two-particle wave function for two identical particles at a time $t$ is given by $\Psi({\bf r}_{1},{\bf r}_{2};t)$. Since the two particles are identical, we have a relation $|\Psi({\bf r}_{1},{\bf r}_{2};t)|=|\Psi({\bf r}_{2},{\bf r}_{1};t)|$. This gives a symmetry relation LeggettBook

[TABLE]

where the upper sign is for boson with integer spin, and the lower sign is for fermion with half-integer spin. In contrast to the fermion case, where $\Psi({\bf r},{\bf r};t)=0$ for ${\bf r}\equiv{\bf r}_{1}={\bf r}_{2}$, which provides the anti-bunching effect, the bosons may give $\Psi({\bf r},{\bf r};t)\neq 0$, which provides the bunching effect for bosons. The fermion obeys the Pauli exclusions principle, where a single quantum state cannot be occupied by two or more fermions. On the other hand, a single quantum state can be occupied by any numbers of bosons.

In the grand canonical formulation, the Bose distribution function — the mean occupation number of a quantum state $i$ for boson — is given by PitaevskiiBook

[TABLE]

where $\epsilon_{i}$ is an energy of a single-particle state $i$, and $\mu$ is the chemical potential that satisfies $\epsilon_{0}>\mu$. Here, $\epsilon_{0}$ is the lowest single-particle energy. Using a density of states for a three-dimensional uniform system

[TABLE]

the number of particles without assuming the BEC is given by

[TABLE]

where $z$ is the fugacity $z\equiv\exp(\beta\mu)$, and

[TABLE]

with the Gamma function $\Gamma(\alpha)$.

By introducing the phase space density $\rho_{\rm sd}\equiv N\lambda_{T}^{3}/V$, we have the relation $\rho_{\rm sd}={\mathcal{G}}_{3/2}(z)$. Since the relation $\epsilon_{{\bf p}=0}=0>\mu$ holds, which leads $0<z<1$, we have the following relation in the limit $z\to 1$, given by

[TABLE]

With lowering temperature or increasing density, however, we should have a high phase space density with $\rho_{\rm sd}=N\lambda_{T}^{3}/V>\zeta(3/2)$. In this case, the relation $N=\lambda_{T}^{-3}V{\mathcal{G}}_{3/2}(z)$ is no longer valid, and we need to reconsider the relation

[TABLE]

The second term with $\mu\to 0$, i.e., $z\to 1$, gives

[TABLE]

which is intensive, where $N^{\prime}$ increases proportionally to $V$. The first term $f(\epsilon_{{\bf p}=0})$ with $\mu\neq 0$ may be negligible compared with the intensive second term $N^{\prime}$ in the thermodynamic limit $V\to\infty$. However, if the chemical potential approaches to zero in the limit $V\to\infty$ with the relation

[TABLE]

the number of particles occupying the lowest energy state $N_{0}=f(\epsilon_{{\bf p}=0})$ may become intensive in the thermodynamic limit $V,N\to\infty$ with $N/V={\rm const.}$, which provides $f(\epsilon_{{\bf p}=0})\simeq k_{\rm B}T/|\mu|=n_{0}V$. The coefficient $n_{0}$ is here found to be the condensate density given by $n_{0}=N_{0}/V$. This is the Bose–Einstein condensation, where the lowest energy state is macroscopically occupied, and its critical temperature for an ideal Bose gas $T_{\rm c}^{0}$ is given by

[TABLE]

III Bogoliubov Prescription

We now consider the non-zero temperature Green’s function formalism. We start with a non-interacting Bose gas, where the Hamiltonian is given by

[TABLE]

Here, $\hat{\Psi}({\bf r})$ and $\hat{\Psi}^{\dagger}({\bf r})$ are an annihilation and creation operator of a bosonic atom at a spatial position ${\bf r}$, satisfying a commutation relation $[\Psi({\bf r}),\Psi^{\dagger}({\bf r}^{\prime})]=\delta({\bf r}-{\bf r}^{\prime})$. The imaginary time Green’s function is a very powerful tool for studying physics at nonzero-temperatures, which is defined as

[TABLE]

where $\tau$ is the imaginary time. The real space representation of the annihilation operator $\hat{\Psi}({\bf r},\tau)$ and its momentum space $\hat{a}_{\bf p}(\tau)$ are related through

[TABLE]

The non-interacting Green’s function in the momentum and Matsubara frequency spaces is given by StoofBook

[TABLE]

where $T_{\tau}$ is the time ordering operator, where the operator at smaller imaginary time is ordered at the right, and $\omega_{n}=2\pi n/(\hbar\beta)$ with $n\in\mathbb{Z}$ is the bosonic Matsubara frequency. This non-zero temperature formalism can be mapped to the absolute zero temperature case with the analytic continuation $i\omega_{n}\to\omega+i\eta$ with an infinitesimally small number $\eta$. The pole of the Green’s function $G_{0}^{-1}(p)=0$ provides the dispersion relation of the excitation, where the ideal gas case provides $\hbar\omega=\epsilon_{\bf p}-\mu$. By taking the sum of the Matsubara frequency and using the counter integration, we obtain the Bose distribution function for an ideal Bose gas, given by StoofBook

[TABLE]

One of the treatments for the Bose–Einstein condensation is the Bogoliubov prescription, where the field operator of the zero-momentum state is replaced by a $c$-number, given by

[TABLE]

Here, $\Phi_{0}$ is the condensate wavefunction, which is related to the condensate density $n_{0}$ and the phase $\varphi_{0}$ through $\Phi_{0}=\sqrt{n_{0}}\exp(i\varphi_{0})$.

Since the field operators of the zero-momentum part satisfy AbrikosovBook ; FetterBook

[TABLE]

the error of the replacement of the operator with the $c$-number may be given by ${\mathcal{O}}(1/V)$, which is expected to be negligibly small in the thermodynamic limit $V\to\infty$. In the BEC phase, the density matrix $\rho_{(1)}({\bf r}_{1},{\bf r}_{2})=N_{0}/V-G_{0}({\bf r}_{1},\tau;{\bf r}_{2},\tau+\eta)$ has the off-diagonal long-range order, which can be given by PitaevskiiBook

[TABLE]

IV identities of interacting condensed Bose gas

We now consider an interacting Bose gas, where the Hamiltonian $\hat{H}$ is given by

[TABLE]

The interaction strength $U$ is related to an $s$-wave scattering length $a$ through a relation

[TABLE]

where $p_{\rm c}$ is a cutoff momentum. In the BEC phase, we apply the Bogoliubov prescription to the field operator. In the Green’s function formalism, interaction effect is included through the self-energy ${\bf\Sigma}(p)$, and the Green’s function in the BEC phase is given by the Dyson–Beliaev equation, where AbrikosovBook ; FetterBook

[TABLE]

with

[TABLE]

and ${\bf G}_{0}(p)={\rm diag}[G_{0}(p),G_{0}(-p)]$.

We here note the lowest contribution of ${\bf G}(p)$. By using the Dyson–Beliaev equation, we have the infinite series of the Green’s function ${\bf G}(p)={\bf G}_{0}(p)+{\bf G}_{0}(p){\bf\Sigma}(p){\bf G}_{0}(p)+{\bf G}_{0}(p){\bf\Sigma}(p){\bf G}_{0}(p){\bf\Sigma}(p){\bf G}_{0}(p)+\cdots$, and the lowest contribution of the normal Green’s function $G_{11}(p)$ is found to be $G_{11}(p)\simeq G_{0}(p)+\cdots$, which means that the first order contribution is the non-interacting Green’s function, which gives $G_{11}(p)={\mathcal{O}}(U^{0})$. On the other hand, the lowest contribution of the anomalous Green’s function is given by $G_{12}(p)\simeq G_{0}(p)\Sigma_{12}(p)G_{0}(-p)+\cdots$. In the lowest contributions of the self-energy $\Sigma_{12}(p)=U\Phi_{0}^{2}$, we can find that the anomalous Green’s function has $G_{12}(p)={\mathcal{O}}(U^{1})$. The order of the contribution of the interaction is different between the normal and anomalous Green’s functions.

The Dyson–Beliaev equation gives the full-Green’s function ${\bf G}(p)=\hbar/[{\bf G}_{0}^{-1}(p)-\hbar{\bf\Sigma}(p)]$, which provides AbrikosovBook ; FetterBook

[TABLE]

where

[TABLE]

with

[TABLE]

Here, we used $\Sigma_{11}(p)=\Sigma_{22}(-p)$ and $\Sigma_{12}(p)=\Sigma_{21}(-p)$. The same relation also holds for the Green’s function $G_{ij}(p)$.

The Hugenholtz–Pines relation Hohenberg1965

[TABLE]

ensures that an excitation of the single-particle excitation in a BEC is gapless. Indeed, (34) gives the denominator of (29) and (30) being zero, i.e., $D(p=0)=0$. This relation can be proven by various ways. One of the key ideas of the proof comes from the fact that the energy of the system is gauge invariant Hugenholtz1959 . A key idea of another proof comes from the linear response of the bosonic field operator with respect to an infinitesimally small symmetry breaking external field Hohenberg1965 . Using the Hugenholtz–Pines relation, the single-particle Green’s function in the low-energy limit is given by Gavoret1964

[TABLE]

where $c$ is the sound speed of the single particle excitation in the low-energy limit, and $n\equiv N/V$ is the total number density of particles.

The sound speed $c$ of the single particle excitation in (35) is equal to the thermodynamic sound speed, given by Gavoret1964 ; GriffinBook

[TABLE]

It indicates that the sound speed of the single-particle excitation is equal to that of the density mode given by the density-density correlation function, Gavoret1964 ; GriffinBook

[TABLE]

Indeed, the pole of this correlation function is originally comes from that of the single-particle Green’s function Gavoret1964 ; GriffinBook . These structures are specific to the BEC system, and important points are (i) the single-particle Green’s function is involved to the density correlation function because of the existence of the condensate, and (ii) the single-particle Green’s function has a pole giving the phonon dispersion relation, whose sound speed is the same as the thermodynamic sound speed.

The Nepomnyashchii–Nepomnyashchii identity Nepomnyashchii1975 ; Nepomnyashchii1978

[TABLE]

is another identity, which is caused by the non-vanishing weak infrared divergence of the higher order correlation functions. This identity provides the weak infrared divergence of the longitudinal susceptibility Nepomnyashchii1983 ; GriffinBook , which is also related to the fact that the main contribution of the excitation in the low-energy limit is the phase fluctuation.

In this paper, the details of the proofs for these exact properties are not shown. However, some useful ideas for them will be summarized.

For the Hugenholtz–Pines theorem, it is important that the system energy is gauge invariant, and the normal and anomalous self-energies in the limit $p\to 0$ is generated from the system energy $E$ per volume by eliminating the condensate wave-function Hugenholtz1959 ; Nepomnyashchii1978 . In the BEC phase, the system energy is given by the sum of the connected diagrams constructed by $G_{0}(p)$ as well as $\Phi_{0}^{(*)}$. Since the system energy is gauge invariant, it does not depend on the choice of the phase of the field operator as well as the condensate wave function. Indeed, the normal Green’s function $G_{0}(p)=-\langle T_{\tau}\hat{a}_{\bf p}(\tau)\hat{a}_{\bf p}^{\dagger}(0)\rangle$, which has the same numbers of $\hat{a}_{\bf p}$ and $\hat{a}_{\bf p}^{\dagger}$, is unchanged by the gauge transformation $(\hat{a}_{\bf p},\hat{a}_{\bf p}^{\dagger})\to(\exp(i\varphi)\hat{a}_{\bf p},\exp(-i\varphi)\hat{a}_{\bf p}^{\dagger})$. In order to keep the system energy gauge invariant, a certain contribution to the system energy should have the same number of $\Phi_{0}$ and $\Phi_{0}^{*}$, which provides the system energy is also unchanged by the gauge transformation $(\Phi_{0},\Phi_{0}^{*})\to(\exp(i\varphi)\Phi_{0},\exp(-i\varphi)\Phi_{0}^{*})$. If a condensate wave function $\Phi_{0}$ or $\Phi_{0}^{*}$ is eliminated from the energy functional, a single vertex point of $p=0$ is generated. This idea can be easily extended to generating multi-points vertex functions. Suppose that a certain contribution to the energy functional $E^{(s)}$ has $\Phi_{0}$ and $\Phi_{0}^{*}$, whose numbers are $s$, respectively. The normal self-energy is generated by eliminating a single $\Phi_{0}$ from $s$ possibilities of $\Phi_{0}$ in $E^{(s)}$, and by eliminating a single $\Phi_{0}^{*}$ from $s$ possibilities of $\Phi_{0}^{*}$ in $E^{(s)}$. On the other hand, the anomalous self-energy is generated by eliminating a single $\Phi_{0}^{*}$ from $s$ possibilities of $\Phi_{0}^{*}$ in $E^{(s)}$, and then eliminating a single $\Phi_{0}^{*}$ from the remaining $s-1$ possibilities of $\Phi_{0}^{*}$. As a result, we have relations

[TABLE]

which leads

[TABLE]

The right hand side, on the other hand, is given by Hugenholtz1959 .

[TABLE]

By employing the Bogoliubov prescription, we may have a term $-\mu_{0}|\Phi_{0}|^{2}=-\mu_{0}n_{0}$ satisfying $\mu=\mu_{0}$ in the original energy $\bar{E}(T,\mu,\mu_{0})$. However, the contribution $-\mu_{0}n_{0}$ is not generated from the energy functional constructed by the Green’s function, which instead gives $E(T,\mu,n_{0})=\bar{E}(T,\mu,\mu_{0})+\mu_{0}n_{0}$. This is a kind of the Legendre transformation, and we find $dE=-SdT-pdV-n^{\prime}d\mu+\mu_{0}dn_{0}$, which gives $\mu=\mu_{0}=\partial E/\partial n_{0}$. After collecting all the possible contributions with respect to $s$, we have the Hugenholtz–Pines relation

[TABLE]

Since the phase of the condensate wave function is often taken as to be $\varphi_{0}=0$, Eq. (42) can be reduced to Eq. (34), where (34) is found to be a specific representation for $\Phi_{0}=\sqrt{n_{0}}$.

From (39), we can generate the relations Nepomnyashchii1978

[TABLE]

Since the phase $\varphi_{0}$ may be taken arbitrary, we take the average of self-energies for the phase of the condensate wave-function, i.e., $\langle\Sigma_{ij}(0)\rangle=(2\pi)^{-1}\int_{0}^{2\pi}d\varphi_{0}\Sigma_{ij}(0)$. It is clearly seen that the anomalous self-energy is gauge dependent, and we can find $\langle\Sigma_{12}(0)\rangle=0$ Watabe2014 . On the other hand, the normal self-energy is gauge independent, and $\langle\Sigma_{11}(0)\rangle\neq 0$. It is also the case for the chemical potential. This is one of the simplest ways for understanding the Nepomnyashchii–Nepomnyashchii identity $\Sigma_{12}(0)=0$.

This identity is also strongly related to the weak infrared divergence of the longitudinal susceptibility, where we here introduce the longitudinal and transverse susceptibilities Weichman1988 ; Giorgini1992 ,

[TABLE]

Here, $\hat{a}_{\perp,{\bf p}}=(\hat{a}_{\bf p}-\hat{a}_{-{\bf p}}^{\dagger})/(2i)$, and $\hat{a}_{\parallel,{\bf p}}=(\hat{a}_{\bf p}+\hat{a}_{-{\bf p}}^{\dagger})/2$ are transverse and longitudinal operators, respectively. The transverse fluctuation is consistent with the phase fluctuation. Although the longitudinal operator may be often referred to as the density fluctuation or the amplitude fluctuation, those are not exactly the same and the careful treatment of those operators are needed, because the longitudinal operator $\hat{a}_{\parallel}=(\hat{a}+\hat{a}^{\dagger})/2$ has the gauge dependence. The amplitude mode, known as the Higgs mode, may be detected by the scalar susceptibility, not the longitudinal susceptibility Podolsky2011 . In the low-energy regime, we have relations

[TABLE]

If the Nepomnyashchii–Nepomnyashchii identity $\Sigma_{12}(0)=0$ holds, the infrared divergence of the longitudinal susceptibility emerges.

The Popov’s hydrodynamic theory is known as an absolute-zero temperature approximation that satisfies both the Hugenholtz–Pines relation and the Nepomnyashchii–Nepomnyashchii identity PopovBook ; Popov1972 ; Popov1979 . In this theory, the bosonic field operator $\hat{\Psi}(x)\equiv\sqrt{n_{0}+\hat{\pi}(x)}e^{i\hat{\varphi}(x)}$ is expanded by hydrodynamic operators $\hat{\pi}$ and $\hat{\varphi}$ PopovBook ; Popov1972 ; Popov1979 ; Dupuis2011 , and the Green’s function ${\bf G}$ is approximately given by the sum of correlation functions $G_{\hat{\pi}\hat{\pi}},G_{\hat{\pi}\hat{\varphi}},G_{\hat{\varphi}\hat{\pi}},G_{\hat{\varphi}\hat{\varphi}}$ as well as the convolution of $G_{\hat{\varphi}\hat{\varphi}}$. By inversely solving the Dyson–Beliaev equation (27) for ${\bf\Sigma}$ with the given ${\bf G}$ in the Popov’s hydrodynamic theory, we can find the self-energy satisfies both the Hugenholtz–Pines relation and the Nepomnyashchii–Nepomnyashchii identity Popov1979 . In this formalism, the leading order contribution of the longitudinal susceptibility is the convolution of $G_{\hat{\varphi}\hat{\varphi}}\propto 1/(\hbar^{2}\omega^{2}-c^{2}|{\bf p}|^{2})$ Popov1979 . This convolution in the low-energy regime shows the weak infrared divergence, given by GriffinBook

[TABLE]

As a result, the anomalous self-energy vanishes as $\Sigma_{12}\simeq 1/(4\chi_{\parallel})\propto|{\bf p}|$ in the limit ${\bf p}\to 0$ at $T\neq 0$.

Appearance of the convolution of the transverse-field (or phase) correlation in the longitudinal susceptibility is originated from the fact that the excitation in the low-energy limit is exhausted by the transverse field (or the phase of BEC), and the amplitude of the order-parameter may not change in this limit Patashinskii1973 . Let $\Phi_{0}=\sqrt{n_{0}}$ be the original field, and in the low-energy excitation with the fluctuation, the field may be approximately given by

[TABLE]

If the amplitude of the order parameter does not change in the low-energy limit, i.e., $|\Phi_{0}^{\prime}|^{2}=|\Phi_{0}|^{2}$, we have a condition

[TABLE]

As a result, the leading order of the longitudinal correlation function $\langle\delta\Phi_{\parallel}\delta\Phi_{\parallel}\rangle$ is given by the transverse fluctuations, i.e, $\langle\delta\Phi_{\perp}^{2}\delta\Phi_{\perp}^{2}\rangle$, which provides the convolution of the transverse correlation function $\langle\delta\Phi_{\perp}\delta\Phi_{\perp}\rangle$. This story is consistent with the infrared divergence of the longitudinal susceptibility and that of the Popov’s hydrodynamic theory.

The original derivation of the Nepomnyashchii–Nepomnyashchii identity Nepomnyashchii1975 ; Nepomnyashchii1978 is different from those approximation theories. The self-energies — two-point vertices — can be generated from the three-point vertices, where two of three vertex points are connected to two single-particle Green’s functions, which are also connected to an interaction line as well as a single condensate line  Nepomnyashchii1978 . Although almost all infrared divergences are canceled out each other, some contributions remains, which is the same contribution as that in (47) Nepomnyashchii1978 ; Watabe2014 . On the other hand, as discussed in the topic of the Hugenholtz–Pines identity, the vertex functions in the zero-energy limit is generated from the system energy by eliminating the condensate wave-functions. As a result, three-point vertices in the low-energy limit can be related to the two-point vertices (i.e., self-energies), and one can obtain an equation with respect to the anomalous self-energy. By solving this equation, one can find the anomalous self-energy vanishes in the zero-momentum and zero-energy limits, which is caused by the inverse of the weak-infrared divergence shown in (47) Nepomnyashchii1978 ; Watabe2014 .

As related to the fact that the self-energy in the zero-energy limit can be generated from the energy functional, the self-energy with the small but nonzero momentum and that with the small but nonzero energy is also generated by the similar way Gavoret1964 . The self-energy contribution in the low-energy regime is thus related to the thermodynamic quantities. Indeed, we have a relation Gavoret1964

[TABLE]

where we used relations $n^{\prime}=-\partial E^{\prime}/\partial\mu$, and $\partial n^{\prime}/\partial\mu=dn/d\mu=n/(mc^{2})$. Here, $c$ is the thermodynamic sound speed. The thermodynamic functional is constructed by the non-interacting single-particle Green’s function $G_{0}$, the condensate wavefunction $\Phi_{0}$ and $\Phi_{0}^{*}$, as well as the interaction parameter $U$. If one increase the energy of $G_{0}$ with the infinitesimally small value, i.e., $G_{0}(p)=\hbar/[\hbar(\omega+\delta\omega)-\epsilon_{\bf p}+\mu]$, we may regard it as the infinitesimally small increment of the chemical potential, which provides the relation $\hbar^{-1}\partial G_{0}/\partial\omega=\partial G_{0}/\partial\mu$ Gavoret1964 . The second derivative of the self-energy with respect to the frequency can be thus related to the second derivative of the energy functional with respect to the chemical potential. Using (34) and (51), one can show that the sound speed of the single-particle excitation in the low-energy limit is given by the thermodynamic sound speed, as discussed in (36) and (37).

The density response function in the BEC phase is constructed by two contributions, one-particle reducible (1PR) part, and one-particle irreducible (1PI) part Gavoret1964 . The 1PR part is specific to the BEC phase, given by

[TABLE]

where ${\bf\Upsilon}$ and ${\bf\Upsilon}^{\dagger}$ are the density vertices that include the single-particle Green’s function into the density response function and vanish above the critical temperature. This vertex is known to show the zero-frequency density vertex identity Nepomnyashchii1978

[TABLE]

and the 1PR part in the low energy regime is given by Gavoret1964

[TABLE]

This implies that $\lim_{{\bf p}\to 0}\chi^{{\rm 1PR}}(0,{\bf p})=0$ Gavoret1964 . On the other hand, the 1PI part in the low-energy limit exactly shows Gavoret1964

[TABLE]

Since the density response function is the sum of the 1PR and 1PI parts, i.e., $\chi=\chi^{\rm 1PR}+\chi^{\rm 1PI}$, the compressibility zero-frequency sum-rule $\chi(0)=-n/(mc^{2})$ is exhausted by the 1PI part.

V Many-body treatment

This section focuses on many-body approaches at non-zero temperatures, such as the random-phase approximation, and the many-body $T$-matrix theory Watabe2013 ; Watabe2014 ; Watabe2018 . For simplicity, we take the convention $V=\hbar=k_{\rm B}=1$ in this section. It may be convenient to construct building blocks for many-body contributions from the single-particle Green’s function in the Hartree–Fock–Bogoliubov–Popov approximation (Shohno model) GriffinBook , given by

[TABLE]

where $\xi_{\bf p}=\epsilon_{\bf p}+Un_{0}$. Here, $\sigma_{1,2,3}$ are the Pauli matrices. This gives the Bogoliubov excitation — the gapless phonon excitation in the low-energy limit —, which captures properties of the exact single-particle Green’s function. One of the building blocks is the correlation function, given by Watabe2013 ; Watabe2014 ; Watabe2018

[TABLE]

where $\otimes$ is the Kronecker product. Since the normal and anomalous Green’s function has the opposite sign in the low-energy regime, i.e., $g_{ij}(p)\simeq(-1)^{i+j}mc_{0}^{2}/(\omega^{2}-c_{0}^{2}{\bf p}^{2})$, where the Bogoliubov phonon sound speed is $c_{0}=\sqrt{Un_{0}/m}$, the infrared divergence of this correlation function $\Pi$ has a simple structure, which can be extracted by the following matrix Watabe2013 ; Watabe2014 ; Watabe2018

[TABLE]

where

[TABLE]

with $\Pi_{14}(p)=-T\sum\limits_{q}g_{12}(p+q)g_{12}(-q)$. The function $\Pi_{14}$ shows the same structure of the infrared divergence in (47), where $c$ is replaced by $c_{0}$. The correlation function $\Pi^{\rm R}\equiv\Pi-\Pi^{\rm IR}$ converges to a finite value in the low-energy limit.

Using this correlation function $\Pi$, the four point vertex $\Gamma$ and the regular part of the density response function $\chi_{\rm R}$ are respectively given by Watabe2013 ; Watabe2014 ; Watabe2018

[TABLE]

where $\langle f_{0}|=(0,1,1,0)$ and $|f_{0}\rangle=(0,1,1,0)^{\rm T}$. This regular part does not show the infrared divergence and converges to Watabe2013 ; Watabe2014 ; Watabe2018

[TABLE]

with

[TABLE]

The random-phase approximation (RPA) includes the density fluctuation into the effective interaction, where the effective interaction $U_{\rm eff}$ is given by Watabe2013 ; Watabe2018

[TABLE]

Using this interaction, we consider the self-energies as

[TABLE]

where $n_{0}^{\prime}\equiv-T\sum_{p}g_{11}(p)e^{i\omega_{n}\delta}$. Since $\chi_{\rm R}(0)$ does not show the infrared divergence, we find $U_{\rm eff}(0)\neq 0$. This random-phase approximation provides $\Sigma_{12}(0)=n_{0}U_{\rm eff}(0)\neq 0$, which does not satisfy the Nepomnyashchii–Nepomnyashchii identity.

The variant of this RPA, which satisfies the Nepomnyashchii–Nepomnyashchii identity, is the simplified RPA (s-RPA), where we include the lowest order of the contribution Watabe2013 ; Watabe2018

[TABLE]

The effective interaction in this s-RPA is given by

[TABLE]

Self-energies are also given by replacing $U_{\rm eff}$ with $U_{\rm eff}^{0}$ in (65) and (66). Since $\chi_{\rm R}^{0}$ shows the infrared divergence, the effective interaction $U_{\rm eff}^{0}$ converges to zeros $U_{\rm eff}^{0}(0)=0$. As a result, the anomalous self-energy satisfies the Nepomnyashchii–Nepomnyashchii identity $\Sigma_{12}(0)=n_{0}U_{\rm eff}^{0}(0)=0$.

The many-body $T$-matrix (MBT) approximation is given by Watabe2013

[TABLE]

where $\Gamma_{11}$ is the $(1,1)$-element of the four-point vertex $\Gamma$ that is connected to the well-known ladder diagrams above $T_{\rm c}$.

The many-body effect is important for considering the critical temperature shift from that of an ideal Bose gas. The critical temperature of the Bose gas in a harmonic trap is known to decrease by the repulsive interaction Giorgini1996 . This fact can be imagined based on the local density approximation. In this approximation, the critical temperature is a monotonic increase function of the peak-density at the center of a harmonic trap. In the presence of the repulsive interaction, this density becomes lower, because the repulsive interaction tends to exclude the other bosons. This scenario cannot applied to a uniform Bose system, in the case where the system volume and the number of particles are fixed. The critical temperature in the uniform system may be shifted by a repulsive interaction in the competition of two effects; the depletion from BEC and the suppression of the fluctuation Kashrnikov2001 . The depletion may make the critical temperature lower, because the particles are excluded from the condensate by the repulsive interaction. The suppression of the fluctuation may make the critical temperature higher. In the first place, an ideal Bose gas has a strong density fluctuation, and its compressibility diverges. On the other hand, by the repulsive interaction, the density fluctuation is suppressed, and the compressibility in the repulsively interacting condensed Bose system is finite with satisfying the zero-frequency compressibility sum-rule. In the weakly interacting regime, effect of the suppression of the fluctuation is dominant, and the critical temperature increases by the weak repulsive interaction Gruter1997 ; Reppy2000 . In the strongly interacting regime, effect of the depletion is dominant, and the critical temperature decreases by the strong repulsive interaction Gruter1997 ; Reppy2000 .

The critical temperature shift in the weak interaction regime is often discussed by the following formula, given by Baym2001

[TABLE]

where $\Delta T_{\rm c}\equiv T_{\rm c}-T_{\rm c}^{0}$, and $\Delta n_{\rm c}(T_{\rm c}^{0})$ is the difference between the critical density of an interacting Bose gas and that of an ideal Bose gas both at the critical temperature of an ideal Bose gas $T_{\rm c}^{0}$. This difference is given by

[TABLE]

Here, we used the Hugenholtz–Pines relation $\mu=\Sigma_{11}(0)$ at the critical temperature, and also used the fact that the chemical potential is zeros at $T_{\rm c}^{0}$ in an ideal Bose gas. If we take the mean-field approximation, such as the Hartree–Fock–Bogoliubov–Popov approximation (Shohno-model), the self-energy $\Sigma_{11}$ does not include the frequency and momentum-dependence. As a result, the right hand side of (72) vanishes, and the mean-field type approximation cannot describe the critical temperature shift Baym2001 . The critical temperatures in those approximations are the same as that of an ideal Bose gas. In order to study the critical temperature shift by the interaction, the frequency or momentum dependent self-energy is needed. The critical temperature shift in the dilute limit is characterized by the following equation Baym2001 ; Andersen2004

[TABLE]

The values $c_{1}$ are summarized in table I. The Monte-Carlo simulation gives $c_{1}\simeq 1.3$ Arnold2001 ; Kashrnikov2001 . The closet value in our many-body approximation is $c_{1}\simeq 1.1$, which is given by the RPA Watabe2013 . Many other values of $c_{1}$ are summarized in the review Andersen2004 .

The density response function can be constructed by the RPA formalism, where the 1PI and 1PR parts are respectively given by Watabe2018

[TABLE]

where $\chi_{\rm R}$ is given in (61), and the density vertices are given by Watabe2018

[TABLE]

with three point vertices ${\bf\gamma}(p)=\Gamma(p)\Pi(p)|f_{0}\rangle$ and ${\bf\gamma}^{\dagger}(p)=\langle f_{0}|\Pi(p)\Gamma(p)$. Here, we introduced the condensate Green’s function $G_{1/2}\equiv\sqrt{-n_{0}}(1,1)^{\rm T}$, $G_{1/2}^{\dagger}\equiv\sqrt{-n_{0}}(1,1)$, ${\mathcal{G}}_{1/2}=\sqrt{-n_{0}}\eta_{\rm g}$, and ${\mathcal{G}}_{1/2}^{\dagger}=\sqrt{-n_{0}}\eta_{\rm g}^{\dagger}$. The matrix $\hat{T}$ and $\eta_{\rm g}$ are respectively given by Watabe2018

[TABLE]

and $\eta_{\rm g}^{\dagger}$ is the transpose of $\eta_{\rm g}$. The vertex coefficient $A(p)$ may be given by Watabe2018

[TABLE]

This is one of the way to include infinite series of the selected diagrams into the density response function in the RPA manner for BEC. In this case, we find that this density vertex converges a non-zero finite value ${\bf\Upsilon}(0)\neq 0$, which does not satisfy the zero-frequency density vertex identity ${\bf\Upsilon}(0)=0$. In order to satisfy the identity ${\bf\Upsilon}(0)=0$, we employ the simple replacement of the vertex coefficient $A(p)$ in (79) by Watabe2018

[TABLE]

Since the one-loop approximation $\chi_{\rm R}^{0}$ shows the infrared divergence, $A(p)$ vanishes in the low-energy limit, which provides the zero-frequency density vertex identity ${\bf\Upsilon}(0)=0$ Watabe2018 .

If the 1PI part is constructed by the infrared divergent one-loop approximation, i.e., $\chi^{\rm 1PI}(p)=\chi_{\rm R}^{0}(p)/[1-U\chi_{\rm R}^{0}(p)]$, the 1PI part in the low-energy limit is given by $\chi^{\rm 1PI}(0)=-1/U$. Using the result of the sum-rule $\chi^{\rm 1PI}(0)=-1/U=-n/(mc^{2})$, the sound speed is given by $c=\sqrt{Un/m}$, which is unphysical temperature-independent sound speed. The approximation (74) rather well reproduces the temperature-dependent sound speed $c=\sqrt{-n/[m\chi^{\rm 1PI}(0)]}$, where Watabe2018

[TABLE]

The self-energy in the simplified-RPA reproduces the Nepomnyashchii–Nepomnyashchii identity Watabe2013 . However, the critical temperature shift $c_{1}$ in this approximation is over-estimated as shown in table I. It will be useful to control the self-energy contributions of the Green’s function so as to satisfy the Nepomnyashchii–Nepomnyashchii identity as well as to reproduce the reasonable value of the critical temperature shift. One of the idea is to construct the Green’s function similar to that of the Popov’s hydrodynamic theory developed in the case at $T=0$ Popov1979 ; Watabe2014 .

In general, the interaction effect can be included to the Green’s function by the Dyson–Beliaev equation (27), with the irreducible self-energy ${\bf\Sigma}$. Another way is to use the non-interacting Green’s function ${\bf G}_{0}$ with the reducible self-energy ${\bf\Sigma}^{\prime}$, given by

[TABLE]

where

[TABLE]

The hybrid version Watabe2014

[TABLE]

is also possible, where $\tilde{\bf G}^{-1}={\bf G}_{0}^{-1}-{\bf\Sigma}_{\rm a}$ and $\tilde{\bf\Sigma}=(1-{\bf\Sigma}_{\rm b}\tilde{\bf G})^{-1}{\bf\Sigma}_{\rm b}$ with ${\bf\Sigma}={\bf\Sigma}_{\rm a}+{\bf\Sigma}_{\rm b}$.

The approximation theory satisfying the Nepomnyashchii–Nepomnyashchii identity is constructed by separating the self-energy into two-parts; one is the regular part $\tilde{\bf\Sigma}^{\rm R}$, where the bare infrared divergent contribution is not included, and the other is the infrared divergent part $\tilde{\bf\Sigma}^{\rm IR}$ Watabe2014 . We include the regular part into the Green’s function by the Dyson–Beliaev type equation, i.e., ${\bf\Sigma}_{\rm a}=\tilde{\bf\Sigma}^{\rm R}$. For example, the regular part of the self-energy in the MBT approximation may be given by Watabe2014

[TABLE]

where $\Gamma_{11}^{\rm R}$ is the $(1,1)$-element of the four-point vertex

[TABLE]

where $\Pi^{\rm R}=\Pi-\Pi^{\rm IR}$. On the other hand, the infrared divergent part ${\bf\Sigma}^{\rm IR}$ is included into $\tilde{\bf\Sigma}$ within the first order, i.e., Watabe2014

[TABLE]

In particular, we take the infrared divergent part as Watabe2014

[TABLE]

The Green’s function in this prescription reproduces the weak infrared divergent longitudinal susceptibility Watabe2014

[TABLE]

Given the approximated Green’s function in the form ${\bf G}=\tilde{\bf G}+\tilde{\bf G}\tilde{\bf\Sigma}\tilde{\bf G}$, we obtain the self-energy ${\bf\Sigma}$ by solving the Dyson–Beliaev equation (27), given in the form Watabe2014

[TABLE]

Since $\tilde{\bf\Sigma}={\bf\Sigma}^{\rm IR}$ shows the infrared divergence, we have ${\bf\Sigma}(0)={\bf G}_{0}^{-1}(0)=\mu$. It indicates that the present prescription can reproduce the Hugenholtz–Pines relation, as well as the Nepomnyashchii–Nepomnyashchii identity. This formulation is consistent with the Popov’s hydrodynamic theory at $T=0$, where the Green’s function is given by the correlation functions of the hydrodynamic variables.

VI Summary

In this paper, starting from the discussion of properties of an ideal Bose gas, exact relations in an interacting condensed Bose gas is discussed. In particular, we focused on the Hugenholtz–Pines relation, Nepomnyashchii–Nepomnyashchii identity, and the density response function. The Hugenholtz–Pines relation is related to the gapless excitation of the single-particle excitation. The Nepomnyashchii–Nepomnyashchii identity gives a weak infrared divergence of the longitudinal susceptibility. The zero-frequency compressibility sum-rule is exhausted by the one-particle irreducible part of the density response function. On the other hand, the one-particle reducible part of the density response function, which is specific to the BEC, has the same pole of the single-particle Green’s function, whose sound speed is equal to that of the thermodynamic sound mode.

In the last part of this paper, we discussed the approaches of many-body approximations satisfying these exact relations. We introduced two random phase approximations and the many-body $T$-matrix theory. We discussed the critical temperature shift by the many-body effect in these approximations, and the systematic approach to reproduce the Hugenholtz–Pines relation, Nepomnyashchii–Nepomnyashchii identity, as well as the weak-infrared divergence of the longitudinal susceptibility. We also discussed the random phase approximation for studying the density response function with satisfying the identity. There is not a many-body approximation satisfying all the exact relations in a Bose–Einstein condensate (BEC). This paper will be useful for developing beyond mean-field theory of an interacting condensed Bose system consistent with exact relations in BECs.

Acknowledgements.

The author thanks Y. Kato, Y. Ohashi, M. Ueda, and T. Nikuni for valuable discussions and comments on BEC and superfluidity. A series of the author’s study on BECs are supported by JSPS KAKENHI Grant No. (249416, JP16K17774, JP18K03499). He also gratefully acknowledges the support by Nanotech Career-up Alliance (Nanotech CUPAL), Japan Science and Technology Agency (JST).