Beliaev technique for a weakly interacting Bose gas

# Beliaev technique for a weakly interacting Bose gas

B Capogrosso-Sansone, S Giorgini, S Pilati, L Pollet, N Prokof’ev, B Svistunov and M Troyer Department of Physics, University of Massachusetts, Amherst, MA 01003, USA Institute for Theoretical Atomic, Molecular and Optical Physics, Harvard-Smithsonian Center of Astrophysics, Cambridge, MA, 02138, USA Dipartimento di Fisica, Università di Trento and CNR-INFM BEC Center, I-38050 Povo, Trento, Italy Theoretische Physik, ETH Zurich, CH-8093 Zurich, Switzerland Physics Department, Harvard University, Cambridge, MA 02138, USA Russian Research Center “Kurchatov Institute”, 123182 Moscow, Russia
###### Abstract

Aiming for simplicity of explicit equations and at the same time controllable accuracy of the theory we present results for all thermodynamic quantities and correlation functions for the weakly interacting Bose gas at short-to-intermediate distances obtained within an improved version of Beliaev’s diagrammatic technique. With a small symmetry breaking term Beliaev’s diagrammatic technique becomes regular in the infrared limit. Up to higher-order terms (for which we present order-of-magnitude estimates), the partition function and entropy of the system formally correspond to those of a non-interacting bosonic (pseudo-)Hamiltonian with a temperature dependent Bogoliubov-type dispersion relation. Away from the fluctuation region, this approach provides the most accurate—in fact, the best possible within the Bogoliubov-type pseudo-Hamiltonian framework—description of the system with controlled accuracy. It produces accurate answers for the off-diagonal correlation functions up to distances where the behaviour of correlators is controlled by generic hydrodynamic relations, and thus can be accurately extrapolated to arbitrarily large scales. In the fluctuation region, the non-perturbative contributions are given by universal (for all weakly interacting U(1) systems) constants and scaling functions, which can be obtained separately—by simulating classical U(1) models—and then used to extend the description of the weakly interacting Bose gas to the fluctuation region. The theory works in all spatial dimensions and we explicitly check its validity against first-principle Monte Carlo simulations for various thermodynamic properties and the single-particle density matrix.

###### pacs:
05.30.Jp, 03.75.Hh, 67.10.-j

## 1 Introduction

For nearly half a century the theory of the weakly interacting Bose gases (WIBG) remained in the realm of purely theoretical investigations [1, 2, 4, 3, 5, 6, 7] (for a recent review, see [8]) providing insight into the nature of superfluid states of matter but not directly relating to existing experimental systems. The situation changed with the realization of Bose Einstein condensation in cold atomic gases [9, 10, 11]. The typical values of the interaction parameter for alkali atoms are very small , where is the number density of the gas and is the -wave scattering length. Since 1995 both the theory of WIBG and experiments have progressed with a close relationship between the two [12].

Until recently, existing mean-field and variational treatments such as the Bogoliubov zero-temperature approximation and the finite-temperature quasi-condensate theory [13] were capable of describing the data within the relatively large experimental uncertainties. However, the improvements in detection techniques, the studies of optical lattice systems in which the effective gas parameter can be made arbitrarily large, and the need for reliable thermometry, e.g., through precise entropy matching, all indicate the need for a more accurate and controllable theory. It is highly desirable, similarly to the well-known corrections to the ground state energy [14] and condensate density [1], to have a theory which accounts for leading corrections to all thermodynamic properties at finite temperature, works in any dimension, and provides estimates for omitted higher order terms.

The accuracy of the description plays an important role in this paper. We systematically estimate the errors of all approximations, which is crucial for comparing seemingly different schemes proposed for Bose gases. Indeed, alternative theories can be equivalent to each other within the level of their accuracy, but a definitive conclusion can only be drawn on the basis of analyzing their systematic errors. For example, within the same accuracy gap-less and gap-full finite-temperature schemes may often be regarded as equivalent without any preference for one of them, as long as the gap is on the order of or smaller than the terms in the chemical potential which the theory ignores. Also, in the fluctuation region all perturbative schemes are equally inaccurate in terms of their treatment of the order parameter, and it makes no sense whatsoever to distinguish between them by the type of the transition they predict.

In this work, we provide a rigorous framework for obtaining a consistent description of all thermodynamic finite- properties of the WIBG in one, two, and three dimensions. It is based on Beliaev’s regularized diagrammatic technique [2] which allows us to calculate all relevant thermodynamic functions of the system. The diagrammatic technique also leads to an accurate description of correlation functions up to sufficiently large distances where Popov’s hydrodynamic description takes over [5]. The key results of our calculations are summarized in section 2, which is accessible to the general reader.

It turns out that all final results for thermodynamic functions perfectly fit into the pseudo-Hamiltonian picture, i.e. we prove that the same answers would be obtained if they were calculated in the standard way with Bogoliubov’s non-interacting quasi-particle picture (with a self-consistently defined temperature dependent spectrum).

In addition, we find an explicit expression for the pressure and give a more accurate expression for the chemical potential. We obtain the equations of state for all basic thermodynamic quantities in an integral parametric form using temperature, effective chemical potential, and the interaction strength as parameters. These results include sub-leading corrections 111The only exceptions are the quantities, like entropy and heat capacity, that vanish as . At low enough temperatures, these quantities lose their ideal-gas leading terms, getting nothing instead—as opposed to energy, chemical potential, and pressure, having straight-forward mean-field contributions that start to lead at low temperatures.. The analysis of systematic errors shows that the accuracy of our results cannot be further improved within the paradigm of non-interacting quasiparticles with a (temperature-dependent) Bogoliubov spectrum. Moreover, there exists a two-parametric continuum of possibilities for recasting the form of the final answers with the same accuracy due to the freedom of modifying the definitions of approximate notions of “effective interaction” and “effective chemical potential”. Using this freedom, one can trade the higher order terms in the explicit integral expressions for the thermodynamic quantities for those of the effective interaction vertex. Our analysis shows that the question of the energy and momentum dependence of the interaction vertex is more a matter of taste rather than accuracy. We show that even in 2D—where, in view of the vanishing -scattering amplitude, a low-energy cutoff is unavoidable and one could expect that the momentum-dependent effective-interaction vertex is crucially important in any theory taking care of sub-leading corrections—the effects associated with the energy and momentum dependence of the effective interaction can be naturally accounted for on equal footing with all the other effects of the same order. Within the systematic error bars of our pseudo-Hamiltonian approach, the effective interaction of the 2D WIBG can be safely treated as an energy and momentum independent constant.

Such analytical approaches are not supposed to work in the vicinity of the critical point . This so-called fluctuation region is characterized by scaling functions that are universal for all U(1) weakly-interacting systems. It has already been studied numerically with high accuracy [15, 16, 17], so that a combination of analytical and numerical data provides an accurate description at all temperatures.

Despite its success in treating the weakly interacting Bose gas, Beliaev’s original technique features a subtlety that has seriously questioned its reliability: the normal and anomalous self energies—the fundamental quantities of the theory—have been demonstrated to be essentially momentum-dependent, with the anomalous self-energy vanishing in the limit of zero momentum and frequency [18, 19, 20]. On the other hand the first-order diagrams in the theory of the weakly interacting gas with a contact interaction yield momentum independent self energies. In the literature, one can find a number of recipes of how to deal with this inherent infrared problem of Beliaev’s technique. In Popov’s approach [5], the diagrammatic technique applies only to the higher-momentum part of the bosonic field, while the lower-momentum part is treated separately within the hydrodynamic representation. Similar results can be achieved by working with a finite-size system and taking advantage of the bi-modality of correlation properties of weakly interacting superfluid bosonic systems (away from the fluctuation region) [21]. The bi-modality allows one to select such a system size that in terms of local thermodynamic properties the system is already well in the macroscopic limit while the effect of the infra-red renormalization of the self energies is still negligibly small. Another option of a controllable microscopic description of weakly interacting Bose gas is the renormalization-group treatment [22, 23]. Similar to the finite-size regularization, our approach is to exploit the bi-modality of correlation properties of the WIBG, see figure 1, by cutting off infra-red divergencies using small U(1) symmetry breaking terms. The amplitude of the symmetry-breaking terms is such that their effect on the local thermodynamic properties of the system is negligibly small while all infrared singularities in Beliaev’s technique are gone. Clearly, the solution to the infrared problem comes at the expense of suppressing long-range fluctuations of the phase of the order parameter (with opening of a gap in the spectrum at small momenta), and thus distorting the long-range behaviour of correlation functions. Nevertheless, this distortion is essentially irrelevant since it takes place at large distances where the (generic to all superfluids) behaviour of correlators is governed by hydrodynamic fluctuations of the phase of the order parameter, with a simple and transparent theoretical description [24]. Moreover, in all final answers one can explicitly take the limit of vanishing symmetry breaking terms.

In terms of the final results, we confine our analysis to the most typical (and most difficult in two and three dimensions) case of the dilute gas, when the size of the potential is much smaller than the distance between the particles:

 R0≪n−1/d. (1.1)

In principle, the theory works also at , and becomes even simpler, since in this case the weakness of interaction literally implies Born type of the potential (while in the opposite case the amplitude of potential can be arbitrarily large in two and three dimensions). However, only with the condition (1.1) the final answers become insensitive to the details of the interaction potential. The second condition which is assumed in our final results is

 λT≫R0, (1.2)

where is the de Broglie wavelength. Given the inequality (1.1), the condition (1.2) is satisfied automatically as long as one is interested in essentially quantum properties of the system, implying . An extension of the theory to the Boltzmann high-temperature regime , where answers become sensitive to the particular form of the interaction potential, goes beyond the scope of this paper.

The paper is organized as follows. In section 2 we first summarize the key results of our approach for thermodynamic quantities and correlation functions in a form easily accessible to the general reader. The next chapters systematically introduce the technique and derive the results. In section 3 we re-derive the diagrammatic technique for the phase with the broken U(1) symmetry by introducing explicit (arbitrary small) symmetry breaking terms in the Hamiltonian. We keep the discussion as general as possible to cover inhomogeneous, e.g. trapped, systems and states with non-zero average current. In section 3 we also derive the central expressions for the proper self energies and chemical potential. In section 4 we consider thermodynamics (including the superfluid density) of a homogeneous system in the superfluid phase. In particular, we introduce the notion of the renormalized chemical potential and show that and temperature form a natural pair of thermodynamic parameters, in terms of which all the other quantities are obtained from the single-particle momentum-space integrals. In section 5 we discuss, at the order-of-magnitude level, the structure of higher-order corrections to the spectrum of elementary excitations and thermodynamic quantities. This analysis, in particular, leads to the conclusion that within the pseudo-Hamiltonian approximation and Bogoliubov ansatz for the spectrum of elementary excitations our results cannot be further improved. In section 6 we show that the answers for the asymptotic long-range behaviour of off-diagonal correlation functions are readily obtained on top of Beliaev’s diagrammatic treatment for medium-range correlations, by employing the generic hydrodynamic description of a superfluid. The short section 7 deals with the normal region. In section 8 we show that in both normal and superfluid regimes (but away from the fluctuation region around the critical point) the thermodynamic relations can be associated with a pseudo-Hamiltonian in the following sense. The partition function and the entropy of the system turn out to correspond to those of a non-interacting bosonic Hamiltonian with temperature-dependent parameters and a temperature-dependent global energy shift. In section 9 we render the results for the fluctuation region, which previously were obtained in 2D and 3D, but not in 1D. In section 10 we finally compare our solutions with first-principle Monte Carlo data for WIBG both in continuous space and on a lattice.

## 2 Key results

Here we summarize the most important relations which we derive in the main part of the paper, assuming that the conditions (1.1) and (1.2) are met. The results of this section work for both continuous-space and lattice systems; in the latter case, the mass is understood as the effective mass of low-energy motion. The equilibrium state of the system is conveniently parametrized by two independent variables, the temperature and the effective chemical potential . The latter is always negative, so that .

### 2.1 Coupling constant

The interaction is characterized by an (effective) coupling constant . In one dimension is the zero-momentum Fourier component of the bare potential . In two and three dimensions, is naturally expressed in terms of the -wave scattering length ; the latter being defined as the radius of hard disk/sphere potential

 Uhard(r)={∞,      r

with low-energy scattering properties identical to those of the potential . A representational freedom allows one to use slightly different expressions for , without sacrificing the order of accuracy. Our particular choice is

 U=4πam     (d=3), (2.2)

in three dimensions, and

 U=4π/mln(2/a2mϵ0)−2γ     (d=2), (2.3)

in two dimensions, where is Euler’s constant. A particular expression for comes from fixing the specific form of an auxiliary function . The latter enters thermodynamic relations in such a way that the values of those remain insensitive—up to irrelevant higher-order corrections—to slightly changing by adopting a slightly different . Equations (2.2) and (2.3) correspond to

 Π(k)=−12ϵ(k)     (d=3), (2.4)
 Π(k)=−121ϵ(k)+ϵ0     (d=2), (2.5)

respectively, where is the particle dispersion relation. In one dimension . Within the representational freedom, different choices of the -pair are straightforwardly connected with each other by (3.49)-(3.50).

In the 2D case, the coupling constant , and, correspondingly, function slowly evolve with changing effective chemical potential and temperature:

 ϵ0≡ϵ0(~μ,T)={(e/2)|~μ|   (TTc). (2.6)

It is well within the representational freedom to multiply the r.h.s. of (2.6) by a constant of order unity, provided the same expression for is used in both (2.3) and (2.5). The only reason for introducing the order-unity factor is to cast the zero-temperature equations of state in the form (4.46)-(4.47) adopted previously by some authors.

In the case when the problem of finding the scattering length is not straightforward, one can resort directly to the integral equation (3.41) defining in terms of and . Numeric solution of this equation can be obtained within bold diagrammatic Monte Carlo technique [25], which gives the scattering length as a by-product of the calculation.

### 2.2 Thermodynamics in the superfluid region

Thermodynamic quantities are obtained in parametric form, as functions of the pair of parameters and . The main relation is for the number density:

 n(~μ,T)=|~μ|/U+n′, (2.7)
 n′=∑k[ϵ(k)−E(k)2E(k)−|~μ|Π(k)+ϵ(k)E(k)NE], (2.8)

where

 E(k)=√ϵ(k)[ϵ(k)+2|~μ|] (2.9)

is the Bogoliubov-type quasiparticle dispersion relation and

 NE=(eE/T−1)−1 (2.10)

is the Bose-Einstein occupation number. In equation (2.8) and below throughout the text where we sum over momenta, we set the system volume equal to unity thus omitting a trivial normalization factor.

The genuine chemical potential is then given by

 μ=nU−|~μ|. (2.11)

The pressure and entropy density are obtained as

 p=~μ2/2U+2n′|~μ|+(n′)2U −12∑k{E(k)−ϵ(k)−|~μ|−~μ2Π(k)+2Tln[1−e−E(k)/T]}, (2.12)

and

 s=∑k[E(k)TNE−ln(1−e−E(k)/T)], (2.13)

from which the energy density follows from the general relation .

At zero temperature the integrals can be evaluated analytically, as discussed in subsection 4.4.

In continuous space, where , the superfluid density is given compactly by Landau’s formula

 ns=n+1d∑kk2m∂N∂E. (2.14)

Due to the condition (1.1) the answer for a lattice differs only by a global factor equal to the ratio between the bare and effective masses.

### 2.3 Single-particle density matrix. Quasi-condensate

In the superfluid region, the single-particle density matrix can be parameterized as

 ρ(r)=~ρ(r)e−Λ(r), (2.15)

with

 ~ρ(r)=n−∑k(1−eikr)ϵE2[E−ϵ−|~μ|2+(E+|~μ|)NE], (2.16)
 Λ(r)=|~μ|n∑k(1−eikr)E−ϵ2E2[1+2NE]. (2.17)

By definition, is the condensate density. The function is finite in any dimensions. In 1D and at finite temperature in 2D, the function diverges at , consistently with the general fact of absence of condensate in those systems. Nevertheless, at distances at which the function saturates to its asymptotic value , the function is still much smaller than unity. This characteristic feature of the weakly interacting system allows one to speak of the quasi-condensate with the density given by . Up to the distances at which becomes , the correlation properties of condensed and quasi-condensed systems are essentially the same. Details of the general approach to long-range off-diagonal many-particle correlation functions are described in section 6.

### 2.4 Thermodynamics in the normal region

The density in the normal region is given by

 n=∑k[e~ϵ(k)/T−1]−1, (2.18)

with

 ~ϵ(k)=ϵ(k)+|~μ|. (2.19)

The pressure is obtained by

 p=n2U−T∑kln[1−e−~ϵ(k)/T]. (2.20)

For the entropy and energy densities we have the relations

 s=∑k[~ϵ(k)TN~ϵ−ln(1−e−~ϵ(k)/T)], (2.21)
 ε=Un2+∑kϵ(k)N~ϵ. (2.22)

### 2.5 Accuracy control and fluctuation region

A necessary condition for the above-outlined relations to apply is the smallness of the parameter

 γ0=⎧⎪ ⎪⎨⎪ ⎪⎩√na3(d=3),mU(d=2),√mU/n(d=1). (2.23)

At , the parameter directly controls the systematic error of the theory, and the condition is sufficient for the theory to be accurate. At higher temperatures, the condition (2.23) is only necessary, since there exists the so-called fluctuation region where the fluctuations of the order parameter are essentially non-linear and are not captured by our theory. The closeness to the fluctuation region is described by the dimensionless parameter

 x=μ−μ(d)c(T)(mdT2U2)14−d. (2.24)

In 2D and 3D systems, is the critical value of the chemical potential for a given temperature (for the explicit expressions, see section 9), in 1D systems, where there is no finite-temperature phase transition (see subsection 5.4 for more detail), . The theory applies as long as , getting progressively less accurate with decreasing . At , the theory fails to properly describe condensate and superfluid densities the values of which being defined by the fluctuating classical-field order parameter. The description of other thermodynamic quantities is better since the fluctuation contributions to them are of the higher order than the leading (and sometimes even sub-leading) terms.

The estimates for systematic errors for major thermodynamic quantities away from and within the fluctuation region are given in sections 5 (in the superfluid region) and 7.3 (in the normal region).

In the fluctuation region and its vicinity, the theory can be improved by incorporating accurate description of fluctuation contributions by—universal to all weakly interacting U(1) models and available in literature—dimensionless scaling functions. The description is outlined in section 9.

## 3 Beliaev diagrammatic technique revisited

### 3.1 The Hamiltonian and approach

The diagrammatic technique for bosons can be derived in terms of functional integrals over the classical complex-valued fields , propagating in the imaginary time from to (we set ), and subject to the -periodicity constraint with respect to the variable [5]. The classical-field grand-canonical Hamiltonian has a form

 H=H0+Hint+H1, (3.1)

where

 H0=12m∫|∇ψ|2dr, (3.2)

is the ideal-gas term, with the particle mass, and

 Hint=12∫U(r1−r2)|ψ(r1)|2|ψ(r2)|2dr1dr2, (3.3)

is the pairwise interaction term. The term contains linear and bilinear terms associated with the chemical potential, , external potential, , and the field which explicitly breaks the U(1) symmetry,

 H1=∫[V(r)−μ]|ψ|2dr−∫[η∗ψ+c.c.]dr. (3.4)

For a homogeneous system . Strictly speaking, at finite one cannot refer to as a chemical potential because the total number of particles is not conserved. However, we employ -terms solely for explicitly breaking the symmetry and stabilizing supercurrent states; at the end of the calculation we take the limit . We will therefore continue referring to as the chemical potential.

Standard finite-temperature treatments of the WIBG typically suffer from infra-red divergencies [18, 19, 20, 26], especially in lower dimensions. The ultra-violet divergence is removed by introducing the notion of the pseudopotential which allows one to express final answers in terms of the scattering length alone. Our treatment successfully and systematically deals with divergence issues, and features (at least) three important considerations:
- First, the derivation is based on the appropriate definition of the pseudopotential for low momenta [see the discussion of equation (3.44) below] in any dimension.
- Second, the U(1)-symmetry breaking field introduces a gap in the spectrum of Goldstone modes in a well-controlled way. This suppresses long-range phase fluctuations of the order parameter and removes infra-red divergencies in the behaviour of correlation functions and self-energies without modifying any of the physically important quantities in the relevant order of approximation. In lower dimensions, finite results in a finite value of the genuine condensate density which dramatically simplifies the description. In this respect, a small does essentially the same job as Berezinskii’s finite-size trick, cf. [21], or Popov’s special (hydrodynamic) treatment of long-wave parts of the fields [5]. In all final answers one can set .
- Third, the effects of both interaction and chemical/external potential are fundamentally non-perturbative separately below the critical temperature. In the absence of interactions a positive chemical potential immediately leads to the instability of the ideal Bose gas, so that there is no way of consistently introducing a positive —crucial for describing the system below the critical point—before switching on the interaction. This observation explains the idea behind omitting all terms in (3.4) from the non-interacting Hamiltonian, . At the end of the day, we will use and as two variables describing the thermodynamic ensemble.

### 3.2 Diagrammatic expansion

Since the non-interacting Hamiltonian corresponds to the ideal gas at chemical potential approaching zero from below, i.e. with no Bose-Einstein condensate, the corresponding (bare) Matsubara Green’s function reads

 G(0)(τ1,τ2,r1,r2)=∑ξ,kG(0)(ξ,k)eik⋅(r1−r2)−iξ(τ1−τ2), (3.5)
 G(0)(ξ,k)=[iξ−ϵ(k)]−1, (3.6)

where is the single particle dispersion relation. We typically assume a parabolic dispersion relation , but our final answers (excluding formulae for superfluid properties which explicitly invoke Galilean invariance) are valid for arbitrary with a parabolic dependence on momentum in the long-wave limit, e.g. for the tight-binding spectrum. We use the Matsubara imaginary-frequency representation: () with the temperature. For graphical representation of diagrammatic expansions and relations we introduce a set of objects in figure 2 depicting the single-particle propagator, the condensate, and various terms in the Hamiltonian.

The non-perturbative response to and is accounted for by Dyson summation. First, we consider diagrams that feature only one incoming or outgoing line—we call them tails. Given our starting point with no condensate in the non-perturbed system such particle-number-changing diagrams exist only due to the symmetry breaking field . The frequency of the line connecting to the rest of the diagram is zero, by frequency conservation. The Dyson summation of all tails attached to a given point replaces them with a single line which we denote as and , if the tails are incoming and outgoing, respectively, see figure 3. Below we write expressions in the frequency representation; if the frequency argument is not mentioned it is implied that its value is zero. We also adopt a convention of integration over repeated coordinate/momentum/frequency arguments. The Dyson equation for then reads:

 ψin(r)=−G(0)(r−r′)η(r′)+ G(0)(r−r′)[V(r′)−μ]ψin(r′)+G(0)(r−r′)Θin(r′). (3.7)

Here is the sum of all other diagrammatic elements attached to the first line which are not accounted for by the first two terms, i.e. excluding diagrams with the field and diagrams connected to the first solid line by the vertex. The subscript ‘in’ reminds that has an extra incoming particle line. Similarly,

 ψout(r)=−η∗(r′)G(0)(r′−r)+ ψout(r′)[V(r′)−μ]G(0)(r′−r)+Θout(r′)G(0)(r′−r). (3.8)

The fact that is a real even function of its argument implies that

 ψout=ψ∗in,Θout=Θ∗in. (3.9)

The diagrammatic expansion for the tail is identical to that for the condensate wave function defined as an anomalous average

 ψ0(r) = ⟨ψ(r)⟩ ≡ ψin(r). (3.10)

Correspondingly, the condensate density is defined as . From now on we will write for and for , and use the notions of condensate lines and tails on equal footing. Of course, in the limit of speaking of the condensate wave-function is meaningful only when the long-range fluctuations of the phase are small.

### 3.3 Normal and anomalous propagators

The next step is to introduce exact, or ‘bold’, particle propagators and work with skeleton diagrams. There are a couple of differences between our approach and the standard Beliaev technique. Our bare (thin-line) propagators have zero chemical potential and this, according to (3.7), immediately results in a constraint relating the condensate density to the chemical potential and [see (3.21) below]. Also, the chemical/external potential have to be explicitly introduced into the standard Beliaev-Dyson equations for the normal and anomalous Green’s functions, see e.g., [27, 28, 26]. These equations can be defined purely diagrammatically. To this end—proceeding in the frequency representation for the sake of definiteness—we introduce the normal Green’s function, , defined as the sum—with a global minus sign, which is a mere convention—of all diagrams which have an incoming -line with frequency to point and an outgoing -line (with the same frequency) from point . The anomalous Green’s function, , by definition, has an incoming -line with frequency to point and another incoming -line with frequency , by conservation of frequency, to point . The anomalous Green’s function is a counterpart of the function : instead of two incoming it has two outgoing -lines; one with frequency from point and another with frequency from point . The symmetry under exchanging the end points of the anomalous Green’s functions immediately implies the following relations

 Fin(ξ,r1,r2)=Fin(−ξ,r2,r1), (3.11)
 Fout(ξ,r1,r2)=Fout(−ξ,r2,r1). (3.12)

Since complex conjugation is equivalent to changing the sign of the Matsubara frequency and direction of propagation we also have

 [G(ξ,r1,r2)]∗=G(−ξ,r2,r1), (3.13)
 F∗in(ξ,r1,r2)=Fout(−ξ,r2,r1)=Fout(ξ,r1,r2). (3.14)

The physical meaning of , , and follows from the structure of the two-point correlation functions in the imaginary-time–coordinate representation:

 ⟨ψ(τ1,r1)ψ∗(τ2,r2)⟩=−G(τ1−τ2,r1,r2)+ψ0(r1)ψ∗0(r2), (3.15)
 ⟨ψ(τ1,r1)ψ(τ2,r2)⟩=Fin(τ1−τ2,r1,r2)+ψ0(r1)ψ0(r2). (3.16)

These relations can be readily checked by expanding the averages into diagrammatic series. The special case of (3.15), corresponding to and , relates the local density to the normal Green’s function and the condensate density:

 n(r)=−G(τ=−0,r,r)+n0(r). (3.17)

 G(ξ,r1,r2)=G(0)(ξ,r1,r2)+G(ξ,r1,r′)[V(r′)−μ]G(0)(ξ,r′,r2) +G(ξ,r1,r′)Σ11(ξ,r′,r′′)G(0)(ξ,r′′,r2) +Fin(ξ,r1,r′)Σ20(ξ,r′,r′′)G(0)(ξ,r′′,r2), (3.18)
 Fin(ξ,r1,r2)=Fin(ξ,r1,r′)[V(r′)−μ]G(0)(−ξ,r2,r′) +Fin(ξ,r1,r′)Σ11(ξ,r′,r′′)G(0)(−ξ,r2,r′′) +G(ξ,r1,r′)Σ02(ξ,r′,r′′)G(0)(−ξ,r2,r′′), (3.19)

with the standard definition of self-energies ’s as sums of diagrams which can not be cut through a single or line. Equations (3.18) and (3.19) are shown graphically in figure 5 for a homogeneous system in momentum representation. The complexity of the theoretical solution is in the evaluation of the and functions.

### 3.4 The chemical potential and the Hugenholtz-Pines relation

Since the bare Green’s function at zero frequency is identical to the inverse Laplacian operator one can cast (3.7) in the differential form

 −Δ2mψ0(r)+[V(r)−μ]ψ0(r)+Θin(r)=η(r). (3.20)

This equation reduces to the Gross-Pitaevskii equation at low enough temperature when the leading term in is . In the homogeneous case at , when and are coordinate-independent, equation (3.20) simplifies to

 μ=Θin/ψ0≡Θout/ψ∗0, (3.21)

generalizing relations (21.3)-(21.5), and , in [28] to finite temperatures.

There is also an exact relation between , , , and (see also [4]). Here (and only here!) we assume that all diagrams for ’s are in terms of and condensate lines. Let be the sum of diagrams contributing to with incoming and outgoing condensate lines. Then (for )

 Σ11(r,r′)ψ0(r′)=∞∑l=1lD(l)in, (3.22)

because each diagram, with incoming condensate lines, contributing to produces—upon integration over with the weight —a diagram contributing to , and there are such diagrams contributing to . An identical argument leads to

 Σ02(r,r′)ψ∗0(r′)=∞∑l=2(l−1)D(l)in. (3.23)

By subtracting (3.23) from (3.22) we obtain

 Σ11(r,r′)ψ0(r′)−Σ02(r,r′)ψ∗0(r′)=Θin(r). (3.24)

In the homogeneous case . Then, in the limit, equations(3.24) and (3.21) can be combined to yield the Hugenholtz-Pines relation [4] (and its finite-temperature version [3])

 μ=Σ11(0,0)−Σ02(0,0). (3.25)

### 3.5 Beliaev-Dyson equations in the presence of homogeneous superflow

In order to discuss the superfluid properties of the homogeneous system, we add a phase factor to the symmetry breaking field

 η(r)=η0eik0⋅r. (3.26)

which readily translates into phase of the condensate wave function

 ψ0(r)=√n0eik0⋅r, (3.27)

The only difference with the previous discussion is that now we have to associate a finite momentum carried by the condensate lines and modify the momentum conservation laws accordingly. Then (3.21) and (3.25) become

 μ=Θ/√n0+k20/2m−η0/√n0, (3.28)
 Σ11(k0)−Σ02(0)=Θ√n0=μ−k202m+η0√n0, (3.29)

where .

It is convenient (for transparency of expressions which follow) to combine frequency and momentum into a single “4-momentum” variable and to introduce an auxiliary momentum . The symmetry between the two ends of the anomalous Green’s functions and, equivalently, between the two ends of the anomalous self-energies is then expressed by (accounting for the momentum carried by the condensate lines)

 Fin/out(P)=Fin/out(P′), (3.30)
 Σ20/02(P)=Σ20/02(P′). (3.31)

[In a more comprehensive notation scheme one has to mention momenta of both incoming lines in .]

Complex conjugation of propagators and condensate lines changes the signs of their 4-momenta. This property can be used to prove the symmetry relation for the Green’s function

 G∗(P)|k0=G(−P)|−k0. (3.32)

Similar symmetry relations take place for the anomalous Green’s functions and all three self-energies.

In the momentum representation, inverting the direction of does not change the analytical expression for propagators. Hence, by inverting the direction of all the lines, including the condensate ones, it follows that

 Fin(P)=Fout(P)≡F(P), (3.33)
 Σ20(P)=Σ02(P)≡~Σ(P). (3.34)

(For the normal Green’s function, , and the self-energy, , inverting momentum directions results in the same series, and in this sense is trivial.)

We are ready to formulate the pair of Beliaev-Dyson equations in the momentum representation. Using symmetry properties into consideration and the shorthand notation we obtain

 G(P)=G(0)(P)+G(P)[Σ(P)−μ]G(0)(P) +F(P)~Σ(P)G(0)(P), (3.35) F(P)=F(P)[Σ(P′)−μ]G(0)(P′)+G(P)~Σ(P)G(0)(P′). (3.36)

The solution in terms of self-energies reads

 G(P)=iξ+ϵ(|2k0−k|)+Σ(P′)−μD(P), (3.37)
 F(P)=−~Σ(P)D(P), (3.38)

where

 D=~Σ2(P)−[ϵ(|2k0−k|)+Σ(P′)−μ+iξ][ϵ(k)+Σ(P)−μ−iξ]. (3.39)

With these relations at hand, one can calculate the current density induced by the phase gradient in the condensate wave-function, see subsection 4.5.

### 3.6 Low-density limit in 3D and 2D: Pseudo-potentials and scattering lengths

In two and three dimensions, the expansion in terms of the bare interaction potential can be (and in most realistic cases, is) non-perturbative. The system is regarded as weakly-interacting only because of the low density of particles. In one dimension, the physics is perturbative instead in the high-density limit for a given potential. Theoretically, dealing with the strong bare potential implies summation of an infinite sequence of ladder diagrams and produces an effective interaction in the form of the four-point vertex [2], , see figure 6. The analytical relation behind figure 6 reads

 Γ(P1,P2,Q)= U(q)+∑KU(q−k)G(P1+K)G(P2−K)Γ(P1,P2,K) (3.40) ≡U(q)+∑KΓ(P1+K,P2−K,Q−K)G(P1+K)G(P2−K)U(k).

When the bare-interaction lines are replaced with ’s, the rest of the series becomes perturbative (excluding the critical region). On the technical side, working with ’s is convenient as long as one does not intend to systematically take into account higher-order corrections, utilizing thus only the (simple and transparent) leading-order expression for . For the higher-order corrections, equation (3.40) involves three different length scales: the size of the potential, , the healing length of the Bose condensed system, and the de Broglie wavelength, while only the non-perturbative physics at the scale requires summation of the ladder diagrams. Hence, it makes sense to construct an object with a much simpler structure than that captures the non-perturbative physics (and thus coincides with the leading-order expression for ), and then systematically investigate the difference between this object and . To achieve this goal let us define the pseudo-pontential by the equation (see figure 7)

 U(q)=U(q)+∑kU(q−k)Π(k)U(k) (3.41) ≡U(q)+∑kU(q−k)Π(k)U(k),

where is such that when is much larger than the inverse healing length or thermal momentum the value of approaches that of the product summed over the frequency of the 4-vector . That is,

 Π(k)→−12ϵ(k)atk→∞. (3.42)

In this case, , and one can expand the difference in a perturbative series. As a result, we arrive at the diagrammatic technique where instead of thin dashed lines standing for the bare interaction potential we have bold dashed lines representing the pseudo-potential , with an additional requirement that whenever two (normal) Green’s function lines are sandwiched between two pseudo-potential lines, the former are supposed to be summed over the frequency difference (which is always possible in view of the frequency-independence of the pseudo-potential), and then has to be subtracted from the result of summation. Specifically, if and are the two external 4-momenta of the above mentioned diagrammatic element, then the following replacement is supposed to take place for internal propagator lines

 ∑ξ(K)G(P1+K)G(P2−K)→∑ξ(K)G(P1+K)G(P2−K)−Π(k), (3.43)

where is the frequency of the 4-momentum .

As long as we are not interested in the non-universal ultraviolet corrections, we can replace with , the systematic error introduced by the replacement being controlled by the following dimensional estimate

 U(q)=U(0)[1+O(q2R20)]. (3.44)

One may wonder how this estimate is reconciled with the momentum dependence of , which is of the order in 3D, and of the order in 2D. The solution is provided by (3.43) explaining that this dependence, which is both universal and perturbative, is taken into account by the second-order ladder-type diagram in .

In 3D, a natural choice for is

 Π(k)=−12ϵ(k)(d=3), (3.45)

in which case we have the usual [27, 28] (see also below)

 U(0)=4πam(d=3), (3.46)

with the -wave scattering length. In 2D one cannot use (3.45) because of the infra-red divergence of the integral in (3.41). A reasonable choice here is

 Π(k)=−121ϵ(k)+ϵ0(d=2), (3.47)

where is an arbitrary low-energy cutoff. The particular value of is not important since final answers are not sensitive to it. Within the systematic-error bars of the pseudo-Hamiltonian description discussed below, any choice of such that is equally reasonable in terms of accuracy, provided the temperature is not much larger than ; moreover, replacing with is also acceptable. (An optimal choice of in the regime will be discussed in subsection 7.2.) There is also no need to introduce in . We will assume that formally in in which case our final answers can be used as written in all spatial dimensions.

Given that (3.45) and (3.47) are not unique [there is a free parameter in (3.47)], it is instructive to explicitly relate two pseudo-potentials, and —corresponding to and , respectively—to each other. We notice that (3.41) implies

 U2(q)=U1(q)+∑kU1(q−k)[Π2(k)−Π1(k)]U2(k). (3.48)

In view of (3.44) this simplifies [up to terms that we neglect in what follows] to

 U2=U1+U1C12U2, (3.49)
 C12=∑k[Π2(k)−Π1(k)]. (3.50)

If and are defined by (3.47) with different cut-offs and , we have

 C12=m4πlnϵ(2)0ϵ(1)0(d=2). (3.51)

The right choice of the functions and implies that and differ only by sub-leading terms. Thus, we can expand the r.h.s. of (3.49) in powers of :

 U2=U1+U21C12+…. (3.52)

Formally, the expansion in terms of the pseudo-potential is perturbative only as long as we exclude the high-momentum contributions to the -body diagrams generating, upon complete summation, -body scattering amplitudes. In a dilute gas, the corresponding diagrams are small (contain extra powers of the gas parameter ) and are neglected in this manuscript.

The expression (3.46) has the meaning of mapping a dilute three-dimensional system with an arbitrary short-range interaction potential onto a system with the interaction potential (2.1), so that the pseudo-potentials of the two systems coincide. The same approach is possible (and popular) in 2D, the parameter being called the two-dimensional scattering length, since the mapping applies to scattering properties as well. For a given potential , the value of can be obtained either from the asymptotic behaviour of the pseudo-potential at appropriately small , or from the asymptotic behaviour of the scattering amplitude at  [28]. Both asymptotic behaviours are closely related to each other since the large- limit of the kernel in (3.41) coincides with the large- limit of the scattering kernel

 Πsc(k,k′)=mk′2−k2+i0 (3.53)

in the integral equation for ,

 f(k,k′)=U(k−k′)+∑qU(k−q)Πsc(q,k′)f(q,k′) (3.54) ≡U(k−k′)+∑qf(k,q)Πsc(q,k′)U(q−k′).

Moreover, a direct relationship between and is readily obtained by noticing that (3.54) and (3.41) imply [cf. (3.48)]

 f(k,k′)=U(k−k′)+∑qU(k−q)[Πsc(q,k′)−Π(q)]f(q,k′), (3.55)

which dramatically simplifies upon replacing with , in accordance with (3.44). In this case, is -independent, and for we get

 f(k′)=U+Uf(k′)∑q[Πsc(q,k′)−Π(q)]. (3.56)

Substituting (3.47) and (3.53) into (3.56) and performing the integral, we find

 f(k′)=2π/mln(√2mϵ0/k′)+2π/mU+iπ/2. (3.57)

Comparing this to the known hard-disk result (see, e.g., [29])

 f(k′)=2π/mln(2/ak′)−γ+iπ/2, (3.58)

where is Euler’s constant, we conclude that

 U=4π/mln(2/a2mϵ0)−2γ. (3.59)

In 3D, equation (3.56)—with given by (3.45)—yields

 f(k′)=U1−imk′U/4π, (3.60)

thus leading to (3.46) by comparison with known expression for the hard-sphere scattering amplitude.

## 4 Thermodynamic functions in the low-temperature region

### 4.1 Basic relations and notions

Explicitly calculating the lowest-order diagrams shown in figure 8 and utilizing (3.17) one finds

 Σ(P)=−2G(r=0,τ=−0)U+2n0U=2nU, (4.1)
 ~Σ(P)=n0U. (4.2)

We see that within the first approximation, both and turn out to be momentum- and frequency-independent. [It is easy to check that the next-order diagrams inevitably introduce momentum and frequency dependence to self-energies and drastically change the structure of the theory.] At this level of accuracy, the chemical potential equals to according to the Hugenholtz-Pines relation. As mentioned above, it is extremely convenient to use an effective chemical potential

 ~μ=μ−2nU, (4.3)

as a thermodynamic variable to characterize properties of the WIBG. [Note that is negative.] Within the same accuracy we can substitute with and simplify expressions for and to

 G(P)=−iξ+ϵ(k)+|~μ|ξ2+E2(k), (4.4)
 F(P)=|~μ|ξ2+E2(k), (4.5)
 E2(k)=ϵ(k)[ϵ(