Beliaev technique for a weakly interacting Bose gas

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.


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 na 3 ∼ 10 −6 , where n is the number density of the gas and a is the s-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 quasicondensate 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-T 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 quasiparticle 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  . Sketch of the density matrix dependence on distance in the WIBG at temperatures below the fluctuation region. After an initial drop the density matrix levels-off at the healing length and undergoes slow decay up to a distance where the superfluid hydrodynamic approach becomes applicable. The U(1) symmetry breaking terms make the density matrix to converge to a finite value (in any dimension) at the hydrodynamic length scales.
results include sub-leading corrections ‡. 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 s-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 subleading 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 T c . 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 R 0 is much smaller than the distance between the particles: In principle, the theory works also at R 0 n −1/d , 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 where λ T 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 λ d T n 1. An extension of the theory to the Boltzmann hightemperature regime λ T R 0 , 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 temperaturedependent 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.

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 m 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 T and the effective chemical potentialμ. The latter is always negative, so thatμ ≡ −|μ|.

Coupling constant
The interaction is characterized by an (effective) coupling constant U . In one dimension U is the zero-momentum Fourier component of the bare potential U(r). In two and three dimensions, U is naturally expressed in terms of the s-wave scattering length a; the latter being defined as the radius of hard disk/sphere potential with low-energy scattering properties identical to those of the potential U(r). A representational freedom allows one to use slightly different expressions for U , without sacrificing the order of accuracy. Our particular choice is in three dimensions, and respectively, where ǫ(k) is the particle dispersion relation. In one dimension Π ≡ 0. Within the representational freedom, different choices of the (Π, U )-pair are straightforwardly connected with each other by (3.49)-(3.50).
In the 2D case, the coupling constant U , and, correspondingly, function Π(k), slowly evolve with changing effective chemical potential and temperature: 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 ǫ 0 is used in both (2.3) and (2.5). The only reason for introducing the order-unity factor e/2 is to cast the zerotemperature 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 a is not straightforward, one can resort directly to the integral equation (3.41) defining U in terms of U and Π. Numeric solution of this equation can be obtained within bold diagrammatic Monte Carlo technique [25], which gives the scattering length as a byproduct of the calculation.

Thermodynamics in the superfluid region
Thermodynamic quantities are obtained in parametric form, as functions of the pair of parametersμ and T . The main relation is for the number density: where is the Bogoliubov-type quasiparticle dispersion relation and 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 The pressure and entropy density are obtained as from which the energy density follows from the general relation ε = µn + T s − p.
At zero temperature the integrals can be evaluated analytically, as discussed in subsection 4.4.
In continuous space, where ǫ(k) = k 2 /2m, the superfluid density is given compactly by Landau's formula (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.

Single-particle density matrix. Quasi-condensate
In the superfluid region, the single-particle density matrix can be parameterized as By definition, ρ(∞) is the condensate density. The functionρ(∞) is finite in any dimensions. In 1D and at finite temperature in 2D, the function Λ(r) diverges at r → ∞, consistently with the general fact of absence of condensate in those systems. Nevertheless, at distances at which the functionρ(r) saturates to its asymptotic valuẽ ρ(∞), the function Λ(r) 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 ∼ 1, 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.

Thermodynamics in the normal region
The density n ≡ n(μ, T ) in the normal region is given by The pressure is obtained by For the entropy and energy densities we have the relations

Accuracy control and fluctuation region
A necessary condition for the above-outlined relations to apply is the smallness of the parameter (2.23) At T nU , the parameter γ 0 directly controls the systematic error of the theory, and the condition γ 0 ≪ 1 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 nonlinear and are not captured by our theory. The closeness to the fluctuation region is described by the dimensionless parameter (2.24) In 2D and 3D systems, µ c (T ) 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), µ (1) c ≡ 0. The theory applies as long as |x| ≫ 1, getting progressively less accurate with decreasing |x|. At |x| 1, 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.

The Hamiltonian and approach
The diagrammatic technique for bosons can be derived in terms of functional integrals over the classical complex-valued fields ψ(r, τ ), propagating in the imaginary time from τ = 0 to τ = β ≡ 1/T (we set = 1), and subject to the β-periodicity constraint with respect to the variable τ [5]. The classical-field grand-canonical Hamiltonian has a form where is the ideal-gas term, with m the particle mass, and is the pairwise interaction term. The H 1 term contains linear and bilinear terms associated with the chemical potential, µ, external potential, V (r), and the field η(r) which explicitly breaks the U(1) symmetry, For a homogeneous system V (r) = 0. 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 |η| → 0. We will therefore continue referring to µ as the chemical potential. Standard finite-temperature treatments of the WIBG typically suffer from infrared 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 U (k) ≈ U (0) 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 Figure 2. Graphical objects representing the single particle propagator G (0) (r − r ′ ), the external field V (r) − µ , the symmetry breaking fields η * (r) and η(r), the condensate lines ψ in (r) and ψout(r), and the interaction U (r − r ′ ) from left to right. Here, and in all other figures in the paper, we will only write down the coordinates or the momenta of the diagrammatic elements, omitting the frequency or imaginary time dependence.
-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, H 0 . At the end of the day, we will use µ − 2nU (0) and T as two variables describing the thermodynamic ensemble.

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 where ǫ(k) is the single particle dispersion relation. We typically assume a parabolic dispersion relation ǫ(k) = k 2 /2m, but our final answers (excluding formulae for superfluid properties which explicitly invoke Galilean invariance) are valid for arbitrary ǫ(k) 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: ξ ≡ ξ s = 2sπT (s = 0, ±1, ±2, . . .) with T 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 H int and H 1 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 nonperturbed 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 ψ in (r) and ψ out (r), 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 ψ in (r) then reads: Here Θ in 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 [V (r) − µ] vertex. The subscript 'in' reminds that Θ in has an extra incoming particle line. Similarly, The fact that G (0) (r) is a real even function of its argument implies that The diagrammatic expansion for the tail is identical to that for the condensate wave function defined as an anomalous average (3.10) Correspondingly, the condensate density is defined as n 0 (r) = |ψ 0 (r)| 2 . From now on we will write ψ 0 for ψ in and ψ * 0 for ψ out , and use the notions of condensate lines and tails on equal footing. Of course, in the limit of η → 0 speaking of the condensate wave-function is meaningful only when the long-range fluctuations of the phase are small.

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 Θ in [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, G(ξ, r 1 , r 2 ), defined as the sum-with a global minus sign, which is a mere convention-of all diagrams which have an incoming G (0) -line with frequency ξ to point r 1 and an outgoing G (0) -line (with the same frequency) from point r 2 . The anomalous Green's function, F in (ξ, r 1 , r 2 ), by definition, has an incoming G (0) -line with frequency ξ to point r 1 and another incoming G (0) -line with frequency −ξ, by conservation of frequency, to point r 2 . The anomalous Green's function F out (ξ, r 1 , r 2 ) is a counterpart of the function F in (ξ, r 1 , r 2 ): instead of two incoming it has two outgoing G (0) -lines; one with frequency ξ from point r 1 and another with frequency −ξ from point r 2 . The symmetry under exchanging the end points of the anomalous Green's functions immediately implies the following relations Since complex conjugation is equivalent to changing the sign of the Matsubara frequency and direction of propagation we also have 14) The physical meaning of G, F in , and F out follows from the structure of the twopoint correlation functions in the imaginary-time-coordinate representation: These relations can be readily checked by expanding the averages into diagrammatic series. The special case of (3.15), corresponding to τ 1 → τ 2 − 0 and r 1 = r 2 , relates the local density to the normal Green's function and the condensate density: The Beliaev-Dyson equations then read with the standard definition of self-energies Σ's as sums of diagrams which can not be cut through a single G or F 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.

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 This equation reduces to the Gross-Pitaevskii equation at low enough temperature when the leading term in Θ in is ∝ |ψ 0 (r)| 2 ψ 0 (r). In the homogeneous case at η → 0, when ψ 0 and Θ in are coordinate-independent, equation (3.20) simplifies to generalizing relations (21.3)-(21.5), and µ = 0|∂(H int /V )/∂n 0 |0 , in [28] to finite temperatures. There is also an exact relation between Σ 11 , Σ 02 , Θ in , and ψ 0 (see also [4]). Here (and only here!) we assume that all diagrams for Σ's are in terms of G (0) and condensate lines. Let D (l) in be the sum of diagrams contributing to Θ in with l incoming and l − 1 outgoing condensate lines. Then (for ξ = 0) because each diagram, with l − 1 incoming condensate lines, contributing to Σ 11 (r, r ′ ) produces-upon integration over r ′ with the weight ψ 0 (r ′ )-a diagram contributing to D (l) in , and there are l such diagrams contributing to Σ 11 . An identical argument leads to By subtracting (3.23) from (3.22) we obtain In the homogeneous case ψ 0 = ψ * 0 . Then, in the η → 0 limit, equations(3.24) and (3.21) can be combined to yield the Hugenholtz-Pines relation [4] (and its finitetemperature version [3]) (3.25)

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 which readily translates into phase of the condensate wave function The only difference with the previous discussion is that now we have to associate a finite momentum ±k 0 carried by the condensate lines and modify the momentum conservation laws accordingly. Then (3.21) and (3.25) become It is convenient (for transparency of expressions which follow) to combine frequency and momentum into a single "4-momentum" variable P = (ξ, k) and to introduce an auxiliary momentum P ′ = (−ξ, 2k 0 − k). 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) [In a more comprehensive notation scheme one has to mention momenta of both incoming lines in F in/out (P, P ′ ).] 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 (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 k 0 does not change the analytical expression for propagators. Hence, by inverting the direction of all the lines, including the condensate ones, it follows that (For the normal Green's function, G(P ), and the self-energy, Σ 11 (P ), 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 Σ ≡ Σ 11 we obtain The solution in terms of self-energies reads 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 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 Γ(0, 0, 0). For the higher-order corrections, equation (3.40) involves three different length scales: the size of the potential, R 0 , the healing length of the Bose condensed system, and the de Broglie wavelength, while only the non-perturbative physics at the scale R 0 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 Γ(0, 0, 0)), and then systematically investigate the difference between this object and Γ. To achieve this goal let us define the pseudo-pontential U (q) by the equation (see figure 7) where Π(k) is such that when k is much larger than the inverse healing length or thermal momentum the value of Π(k) approaches that of the product G(P 1 +K)G(P 2 − K) summed over the frequency of the 4-vector K. That is, In this case, Γ(P 1 , P 2 , Q) ≈ U (q), 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 U , 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 P 1 and P 2 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 where ξ (K) is the frequency of the 4-momentum K.
As long as we are not interested in the non-universal ultraviolet corrections, we can replace U (q) with U (0), the systematic error introduced by the replacement being controlled by the following dimensional estimate One may wonder how this estimate is reconciled with the momentum dependence of Γ, which is of the order qR 0 in 3D, and of the order 1/ ln(1/qR 0 ) in 2D. The solution Figure 7. The diagrammatic expression for the pseudo-potential (left hand side) involves the bare potential (first term on the right hand side) and a modified two-particle propagator Π(k).
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 U . In 3D, a natural choice for Π(k) is in which case we have the usual [27,28] (see also below) with a the s-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 where ǫ 0 is an arbitrary low-energy cutoff. The particular value of ǫ 0 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 ǫ 0 such that nU ǫ 0 n/m is equally reasonable in terms of accuracy, provided the temperature is not much larger than n/m; moreover, replacing ǫ 0 with ck is also acceptable. (An optimal choice of ǫ 0 in the regime T ≫ n/m will be discussed in subsection 7.2.) There is also no need to introduce Π in d = 1. We will assume that formally Π(k) ≡ 0 in d = 1 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, U 1 and U 2 -corresponding to Π 1 and Π 2 , respectively-to each other. We notice that (3.41) implies In view of (3.44) this simplifies [up to O(q 2 R 2 0 ) terms that we neglect in what follows] to The right choice of the functions Π 1 and Π 2 implies that U 1 and U 2 differ only by sub-leading terms. Thus, we can expand the r.h.s. of (3.49) in powers of |U 1 C 12 | ≪ 1: Formally, the expansion in terms of the pseudo-potential is perturbative only as long as we exclude the high-momentum contributions to the (M > 3)-body diagrams generating, upon complete summation, M -body scattering amplitudes. In a dilute gas, the corresponding diagrams are small (contain extra powers of the gas parameter na 3 ) 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 a being called the two-dimensional scattering length, since the mapping applies to scattering properties as well. For a given potential U, the value of a can be obtained either from the asymptotic behaviour of the pseudo-potential U at appropriately small ǫ 0 , or from the asymptotic behaviour of the scattering amplitude f (k, k ′ ) at k, k ′ → 0 [28]. Both asymptotic behaviours are closely related to each other since the large-k limit of the kernel Π(k) in (3.41) coincides with the large-k limit of the scattering kernel Moreover, a direct relationship between U and f (k, k ′ ) is readily obtained by noticing that (3.54) and (3.41) imply [cf. (3.48)] which dramatically simplifies upon replacing U (q) with U (0), in accordance with (3.44). In this case, f (k, k ′ ) is k-independent, and for f ≡ f (k ′ ) we get Substituting (3.47) and (3.53) into (3.56) and performing the integral, we find Comparing this to the known hard-disk result (see, e.g., [29]) where γ = 0.5772 . . . is Euler's constant, we conclude that In 3D, equation (3.56)-with Π(k) given by (3.45)-yields We see that within the first approximation, both Σ andΣ turn out to be momentumand 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 µ = 2nU − n 0 U according to the Hugenholtz-Pines relation. As mentioned above, it is extremely convenient to use an effective chemical potential 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 G and F to , (4.5) A note is in order here. Below we will calculate higher-order corrections to the chemical potential which are necessary for the construction of the accurate equation of state, see (4.9) and (4.11). However, it is still possible to useμ in the definitions of the propagators and the spectrum of elementary excitations while keeping the accuracy of the entire scheme intact, see, subsection 5 where we study systematic errors involved in approximations. Here we mention briefly that the anomalous average contribution to thermodynamic properties is always small, and thus further corrections to F are negligible. The same is true for G at low temperature. At temperatures T ≫ n 0 U , on the other hand, one can ignore tiny modifications in the spectrum because of an additional small parameter n 0 U/T .
For the total density we have where is the Bose-Einstein occupation number for the mode with the energy ε. Formula (4.7) includes the leading and sub-leading terms. To get an expression for the chemical potential with the same degree of accuracy, we need to take into account higher-order diagrams and go beyond the leading-order expressions for Σ andΣ. This is because in the absence of interaction the chemical potential is identically zero and thus the leading term is proportional to U . There are three second-order diagrams contributing to Θ. The first one is the anomalous Green's function convoluted with the bare interaction potential (see figure 9. The second one is the 'sunrise' diagram with the proper correction (3.43) for the ladder structure and the third one is similar to the sunrise diagram but with propagators connecting different vertices (see the two diagrams in the second row in figure 9). The latter two can be safely neglected away from the fluctuation regime because they involve an additional small parameter γ 0 or γ T , see Section 5 for the analysis and definitions of the fluctuation region and parameters γ in different spatial dimensions. It is worth noting that consistently taking into account contributions of the two neglected diagrams in the condensed regime would require simultaneously going to higher orders in the expansions for self-energies, figure 8. The latter, however, is impossible without sacrificing the attractive paradigm of independent quasiparticles with Bogoliubov dispersion.
Keeping the leading and the largest sub-leading diagrams for Θ, we get We have no choice but to use the bare potential here because all ladder diagrams leading to the pseudopotential vertex are already absorbed in the anomalous Green's function, by construction of the latter. This feels unsatisfactory only at first glance since simple formal manipulations allow one to express (4.9) in terms of U alone. Let us introduce an auxiliary function ∆(q) defined by the integral equation (shown graphically in figure 10) . Diagrams contributing to the chemical potential µ of a condensed gas, up to the second order. The last two diagrams [that actually have to be corrected according to (3.43)] are smaller than the previous ones and will be neglected; see the text for more detail. It can be used, in combination with the definition of the pseudopotential Eq.(3.41), to transform the last two terms in (4.9) where ∆ 0 = q ∆(q) . (4.12) The last approximate equality in (4.11) comes from U (q) ≈ U (0) and the observation that ∆(q) vanishes at momenta much smaller than 1/R 0 . The smallness of the parameter |U Πk d | at momenta k ≪ 1/R 0 allows one to expand ∆ [see (4.10)] in the series: For the effective chemical potential we thus get where within our order of accuracy we can take We find it convenient to introduce a quantity with the dimension of density, With this quantity the form of certain thermodynamic relations simplifies. For example, using (4.16) we can replace n 0 with n * in (4.7) and arrive at the result n = n * + n ′ , (4.17) This completes the self-consistent theory because we obtain a closed set of relations which define n = n(µ, T ) in the parametric form-given someμ, or n * , one calculates n from (4.17) and (4.18) and then determines µ from (4.3). The integral in (4.18) is convergent not only in 3D, but also in 2D and 1D. Hence, at this point the quantity η can be set equal to zero. We want to emphasize the fact that all specific expressions derived in this section feature a two-parametric representational freedom, within the same order of accuracy. First, it is possible to useμ with or without higher-order corrections in the spectrum of elementary excitations and propagators. Second, there is a freedom of choosing the function Π(k). For example, one could be tempted to absorb the first two terms in the integral in (4.18) into the definition of U , such that n ′ (T = 0) gets identically equal to zero, and n * (T = 0) equals to the total number of particles. We, however, do not see any merit in this protocol, because U becomes dependent onμ and cannot be considered as a fixed external parameter. Finally, this and analogous "improvements" do not make the theory more/less accurate since the difference is of the same order as omitted diagrams.
A remark is in order here concerning the two-dimensional case, where the value of U cannot be defined irrespectively of the system density. Even in this case, we can proceed with formally independent parameter ǫ 0 in (3.47) till we arrive at the final answers along with the estimates of neglected higher-order terms (the latter being essentially ǫ 0 -dependent). After that, the value of ǫ 0 is selected in such a way that the order-of-magnitude values of neglected terms are minimal.

Pressure
To derive relations for other thermodynamic properties, we start with the pressure as a function of T andμ. Using the general thermodynamic formula n = ∂p(µ, T ) ∂µ , (4.19) and adopting-throughout this subsection-the convention that T is treated as a fixed constant, so that partial derivatives with respect to either µ orμ can be replaced with ordinary ones, we write (after doing straightforward integrations by parts). The integral in (4.23) is readily done by noticing that and thus where ǫ(k) makes the first term convergent. The result of the integration is and we finally get

Entropy and energy
Whenever n, µ, and p are specified as functions of (T, x), where x is a quantity of arbitrary nature, the expressions for entropy per unit volume, s, and energy per unit volume, ε, are readily found from the following two generic thermodynamic relations: With x ≡μ we thus find

Explicit integrations. T = 0 case
It is useful to note that the following two integrals (we useμ ≡ −|μ|), Here we take into account that within our level of accuracy we can replace |μ| → µ in the arguments of I and I (d) 2 , since the integrals are responsible for sub-leading corrections, while the sub-sub-leading terms are beyond our control.
Performing the integrations, one finds: In 2D, we have a freedom of fine-tuning ǫ 0 to simplify the form of the answer. A natural choice, especially convenient for the T = 0 limit, is to set At T = 0 this translates into very compact grand-canonical expressions [29] |μ| = µ , n = µ U (d = 2, T = 0) , (4.46) It is seen that the simplicity of the form of (4.46) comes at the expense of more sophisticated (fine-tuned) form of the effective interaction (4.48). If, instead, one would opt to simplify the form for the effective interaction, getting rid of sublogarithmic terms: U = (4π/m)/ ln(1/a 2 mµ), then (4.46) would acquire the generic form (4.35)-(4.36). Needless to note that this does not change the sum of leading and sub-leading terms in the equations of state since, by construction, the result cannot depend on the specific choice of ǫ 0 .

Superfluid density
The standard way of calculating the superfluid density within the quasi-particle picture is based on Landau's formula for the normal component density, n n , and the relation n s = n − n n . Here we employ a more general approach based on the current induced by the gradient of the condensate wave-function. By definition, the superfluid density is the linear-response coefficient relating the persistent-current density to the gradient of the phase, ϕ, of the (complex) order parameter field: In the last equality we assumed that ϕ = k 0 · r. On the other hand, the average current can be calculated from the microscopic operator expression in terms of the condensate density and the Green's function where N (k, k 0 ) = G(τ = −0, k)| k0 . (4.51) In the limit of k 0 → 0 we obtain the required relation At first glance equation (4.52) does not resemble Landau's formula at all, since the first term is given by the condensate density, not n. However, the structure of the normal average contribution is such that it has a part which completes n 0 to the total density and a part which equals (with minus sign) to the normal density component. To derive this result, we start with the expression for the Green's function and solve for the roots of the denominator in (4.53) to rewrite the Green's function identically as Now we express frequency sums through the Bose functions to arrive at The rest of the calculation is straightforward. In view of the limit to be taken in (4.52) it is sufficient to keep only terms linear in [ǫ(k) − ǫ(k− 2k 0 )] ≈ 2k · k 0 /m in R 1,2 , α, and N (k)| k0 . We omit the algebra of doing the expansion which leads to the result . (4.58) Finally, using (4.7) to exclude n 0 and integrating by parts to simplify the expression we arrive at Landau's formula (4.59)

Expansion parameters. Estimates for higher-order terms
Let us analyze the structure of small parameters which control the applicability of the diagrammatic technique presented above. As we will see a posteriori, it is sufficient to consider two characteristic limits: (i) the T = 0 case and (ii) the finite-T contributions of diagrams with zero frequencies (i.e., classical-field contributions). The analysis is based on dimensional estimates of the diagrams, the fact that an extra interaction vertex increases the total number of propagators by two, and that the largest contributions come from small momenta (and frequencies), where the normal and anomalous propagators behave as and we do not need to distinguish between them. We will further assume that there is an infra-red momentum cutoff k 1 ≪ k 0 where is directly related to the inverse healing length. The symmetry-breaking field η also changes the behaviour of propagators at k ≪ k 0 , but for our purposes we do not need the explicit form of propagators at finite η because we will keep k 1 large enough to neglect effects originating from finite η.
At T < T c there are two kinds of generic vertices in higher-order diagrams, namely (i) full vertices, where four propagators meet, and (ii) vertices with one condensate line (vertices with two and three condensate lines are all accounted for by the lowest-order diagrams for the self-energies and can not be part of the diagrammatic expansion). A naive definition of the diagram order as the total number of vertices proves inconvenient, since contributions of condensate vertices turn out to be larger than contributions of full vertices. Below we will see that the difference is exactly compensated by the fact that condensate vertices (generically) come in pairs. This suggests to define the diagram order as the sum of the total number of full vertices and half the total number of condensate vertices (plus 1/2 for anomalous diagrams).

T = 0 case
At T = 0, the structure for the dimensionless factor associated with adding an extra full vertex to a diagram is with δ(∆k) and δ(∆ξ) representing the δ-functions taking care of momentum and frequency conservation laws. The dimensional estimate, following from (5.3) and (5.1), reads The structure for the dimensionless factor associated with adding an extra pair of condensate vertices is  [18,19,20] (as usual, a zero power of the momentum cutoff should be understood as a logarithm), and emphasize the usefulness of the cutoff-enforcing field η. For the diagrammatic expansion to be consistent, we need k 1 to be large enough in order to guarantee γ 0 ≪ 1. On the other hand, we do not want the field η to significantly affect the physics of the system, and thus need k 1 ≪ k 0 . Hence, the smallest possible expansion parameter one can afford without distorting the physics corresponds to Here we took into account that |μ(T = 0)| ∼ nU , and expressed U in terms of a in 3D. It is clear that γ 0 , equation (5.7), is the actual expansion parameter for all local thermodynamic quantities. These quantities, by referring to large length scales at T = 0, should mainly have contributions from wavevectors k 0 . Technically, it means that infrared divergencies of individual diagrams have to cancel each other in the final answers and the resulting integrals become convergent at the wavevectors k ∼ k 0 , leading to a well defined expansion in powers of γ 0 (5.7). [Since dimensional estimates do not distinguish between pure powers and powers with logarithmic prefactors, we extend the meaning of the term 'in powers' to the powers with logarithmic pre-factors.] Note that all our relations between the thermodynamic quantities that do not vanish in the T = 0 limit (specifically, n, µ, p, ε, and n 0 ), are accurate up to the firstorder corrections in γ 0 , so that their systematic errors are of the order of γ 2 0 (with possible logarithmic pre-factors). The relation for the entropy is accurate only up to the leading term, which actually has the same quasi-particle origin as the sub-leading terms in all non-vanishing at T → 0 quantities. The same is true for the heat capacity that behaves similarly to entropy at T → 0.

Finite-T zero-frequency contributions
The dimensionless factor associated with adding an extra full vertex to a diagram consisting of zero-frequency propagators scales as For the pair of condensate vertices we have Similarly to the T = 0 case (assuming that |μ| ∼ n 0 U , we get In complete analogy with the T = 0 case, the largest possible k 1 (enforced by η) cannot exceed √ mn 0 U. Assuming that k 1 ≤ √ mn 0 U results in This expression shows that at T ≪ nU the parameter γ 0 dominates over γ T , and thus this temperature regime is equivalent to T = 0. In the crossover regime, T ∼ nU , both γ 0 and γ T are of the same order. In the regime T ≫ nU , the classical-field parameter γ T dominates, and becomes of order unity when n 0 gets sufficiently small. The latter situation corresponds to the fluctuation region, where the perturbative theory becomes inadequate, since it fails to properly describe the non-linear long-wave fluctuations of the classical component of the quantum field. Before hitting the region γ T ∼ 1, one has to switch to the description in terms of scaling functions [15,17] which are universal for all U (1) weakly-interacting theories, see section 10. At T ≫ nU , all our relations for the thermodynamic quantities are accurate up to the first-order corrections in γ T , so that their systematic errors are of the order of γ 2 T (with possible logarithmic pre-factors). Note that at T ≫ nU the leading terms for all the thermodynamic quantities, except for the chemical potential (and with a reservation for the density in 1D where the corresponding condition is T ≫ nU/γ 0 ), are the same as for the ideal gas.
To summarize our analysis, we present below an order-of-magnitude estimate of systematic uncertainties due to omitted higher-order terms for each thermodynamic quantity. It is convenient to start with µ and recall that the omitted term in (4.9) comes from the 'sunrise' diagram for Θ. A dimensional analysis of this diagram gives In the case of entropy, it is convenient to start with the regime T nU when-using (4.29) to relate the uncertainty in s to ∆p, ∆n, and ∆µ-we find We see that at T ∼ nU , in contrast to the behaviour of other thermodynamic functions, the uncertainty in s scales only as the first power of γ 0 (simply because s itself gets of order ∼ γ 0 ). This scaling persists down to T → 0: 18) because, at T ≪ nU , the thermodynamics of the system corresponds to the generic low-temperature behaviour of superfluids, where the leading temperaturedependent contributions are due to phonons. Hence, the first-order correction to the sound velocity, which is of the order of γ 0 as is seen from the expression for the energy, translates into the order-γ 0 correction to the entropy. At T ≪ nU , when thermodynamics is exhausted by dilute non-interacting phonons, the correction to the sound velocity can be found directly from the zero-temperature compressibility, and the accuracy of the expression for s (and other thermodynamic quantities) can be improved. At T ∼ nU , however, the order-γ 0 correction to the quasiparticle dispersion law goes beyond the Bogoliubov ansatz [30]. Note also that improving the value of the sound velocity in the T ≪ nU , while rendering the answers for thermodynamic quantities more accurate, is inconsistent with retaining the Bogoliubov form of the spectrum for all the momenta. Finally, equation (4.32) in combination with the above results yields the estimate for the uncertainty of energy:

Fluctuational contributions
In all spatial dimensions, there is an interval (in terms of temperature, if density is kept fixed, or in terms of chemical potential at fixed temperature, etc.) where γ T becomes of order unity, and the systematic perturbative description breaks down. In 3D and 2D this happens in the vicinity of the superfluid phase transition point. In 1D there is no superfluidity in the strict sense of the word, and no phase transition occurs, but the picture remains very similar to that of 2D and 3D case for a weakly interacting gas, as explained in the next subsection. A detailed discussion of the fluctuation region, including accurate expressions for the critical points, will be presented later, in section 9. Meanwhile, here we want to utilize the fact that the results (5.13)- (5.19) allow one to make order-of-magnitude estimates of the non-perturbative fluctuation contributions to the thermodynamic functions. These estimates are quite important as they show the degree to which the perturbative results of the previous sections are inaccurate in the fluctuation region. As we will see, for some quantities the fluctuation corrections are smaller than the leading ideal-gas contributions. By continuity, the order-of-magnitude estimates for fluctuation contributions can be obtained from (5.13)-(5.19) by simply setting γ T ∼ 1, while for the quantities themselves we can use their mean-field critical expressions. This way we arrive at the following results For the superfluid and condensate densities we get an obvious answer that in both cases the fluctuation contributions are on the order of 100%.

Specifics of 1D system
As long as a 1D system is weakly interacting, i.e. γ T ≪ 1, the notion of the order parameter field with well defined amplitude and the two-component (normal + superfluid) description remain physically meaningful. Since superfluidity is a topological phenomenon, it can be destroyed only by topological defects-phase slips. At γ T ≪ 1 the phase slips are rare events which do not contribute significantly to the local thermodynamic quantities; in correlation functions, phase slips show up only at length scales much larger than 1/k 0 , where their effect can be described at the hydrodynamic level.
When the temperature reaches the characteristic scale of the parameter γ T becomes of order unity. Apart from the lack of a genuine phase transition, the physics in the temperature range T ∼ T fluct is close to that of the fluctuation region in 2D and 3D. Only the long-wave classical-field subsystem of the original quantum-field experiences strong non-linear fluctuations. This leads to nonperturbative contributions to the system thermodynamics which are universal for all weakly interacting one-dimensional U(1) systems, and can be described by universal scaling functions in direct analogy with the fluctuation regions in 2D and 3D [15,17]. At T ≫ T (1D) fluct , the low-momenta part of the classical-field component gets depleted, so that the non-linearity of interactions becomes weak and accurately accounted for within the normal-gas mean-field picture.

Off-diagonal correlations
Beliaev's diagrammatic technique allows one to calculate correlation functions up to distances much larger than the correlation radius r 0 ∼ 1/k 0 . However, addressing the asymptotic long-range behaviour of off-diagonal correlation functions requires special techniques properly accounting for long-wave fluctuations of the Goldstone mode, i.e. the phase of the superfluid order parameter. Since the diagrammatic and hydrodynamic description overlap, one can straightforwardly extend the diagrammatic calculation of off-diagonal correlators obtained up to large enough distances, r ≫ r 0 to arbitrarily large r's. We proceed in the spirit of Popov's hydrodynamic approach [19], but with a significant simplification that the hydrodynamic treatment takes place on top of Beliaev's techniques, and thus does not require any additional modifications. The simplification is possible due to the fact that at large distances only phase fluctuations are important and their effect can be factored out as where, under coarse-graining up to momentum k 1 , is the long-wave part of the phase field. As before, the momentum cutoff k 1 ≪ k 0 should be small enough to guarantee (i) the statistical independence of the fluctuations of Φ from the fluctuations ofψ, and (ii) the Gaussian character of phase fluctuations in the hydrodynamic regime. (A particular value of k 1 is not important and does not appear in final expressions.) In view of (i) and (ii), we have for an m-particle correlator, where the subscripts label the space-time variables. Using (6.2) it can be written as whereK is obtained from K by substituting ψ j →ψ j , for all j's. At distances |r j −r s | ≫ 1/k 0 the correlatorK remains constant, and the remaining dependence on distance is exclusively due to the fluctuations of Φ. [The correlatorK does depend on the coarse-graining momentum k 1 and this dependence is crucial for compensating the k 1 -dependence of the function Λ; in fact its origin is directly related to Λ when coarse-graining is stopped at larger momenta and then taken further to k 1 .] The independence ofK on coordinates and times in the asymptotic limit immediately leads to the relation between correlators at different sets of variables, X and X ′ , provided both are in the asymptotic region. Now if X ′ is within reach of diagrammatic expansions, then (6.6) allows one to "extrapolate" K( X ′ ) to K( X) at arbitrary X in the asymptotic domain. By the structure of (6.6), the difference Λ( X ′ )−Λ( X) is independent of k 1 , as opposed to individual Λ's. Here we assumed for simplicity that all variables in X and X ′ are large enough. The same analysis is readily adapted for cases when only some of the coordinates in the m-particle correlation function are in the asymptotic regime; only these coordinates have to be mentioned in (6.6). The Gaussian character of the field Φ implies, see (6.5), that where r sj = r s − r j , τ sj = τ s − τ j (the same for primed variables), and We see that for the extrapolation of correlation functions to arbitrarily large distances one needs to know only Ξ(r, τ ; r ′ , τ ′ ) which is defined through the average |Φ(k, ξ)| 2 . The latter is readily found from Popov's hydrodynamic action leading to |Φ(k, ξ)| 2 = (n s /m)k 2 + κξ 2 −1 . (6.10) In fact, non-zero frequencies (accounting for the quantization of the fluctuations of the phase) are relevant only at T ≪ nU , meaning that the parameter κ in (6.10) can be taken in its T = 0 limit, At T nU , only the ξ = 0 term should be left in (6.8).
Let us consider now the single-particle density matrix ρ(r) = ψ * (r)ψ(0) . (6.12) Taking into account that ρ(r) ≡ n + G(r = 0, τ = −0) − G(r, τ = −0), from the diagrammatic expressions for G we have (at small and intermediate r's) We can use this expression in any dimension at distances significantly exceeding 1/k 0 . In 2D at finite temperature and in 1D at any temperature, the expression (6.13) ultimately becomes inaccurate, and we have to rely on the above-described procedure of extrapolation. In the asymptotic regime we have where, at T ≪ nU , with c = n s /mκ the sound velocity at T = 0 (here n s = n), while at T nU : Within the relevant orders, equation (6.13) corresponds to the expansion of the exponential in (6.14) in powers of Ξ, up to the term ∝ Ξ included. This immediately suggests that within the same accuracy one can exponentiate Ξ-terms in (6.13) to extend its domain of applicability to much larger distances. The other advantage of proceeding this way is having a physical definition of the quasi-condensate density as the amplitude of the order parameter field in the long-wave limit [13]. The equations that follow have the same accuracy as (6.13) though they do contain artificial higherorder terms which arise from factorizing the correlation function. More specifically, we single out terms in (6.13) which on large distances reproduce Ξ, and then exponentiate them: whereρ(r) = n − d(r), and We have added and subtracted [1 + 2N E ] [μǫ/2E 2 ] to (6.13) to ensure that (i) both ρ(r) and Λ(r) are free from ultra-violet and infra-red divergencies in all dimensions, and (ii) that only small momenta k < k 0 contribute to the phase correlator Λ(r). By comparing (6.19) and (6.15) we see that they coincide at large distances up to leading terms. It means that in 2D and in 1D at zero temperature, equation (6.17) can be trusted up to exponentially large scales, while in 1D at finite temperature T ≥ |μ| it works at least up to distances ∼ n s /mT .
We are now in position to define the quasi-condensate density as the limiting value ofρ(r → ∞): The physical meaning of this relation is the amplitude squared of the order parameter field at large distances. There is a certain degree of freedom in attributing terms which do not result in the power-law or exponential decay of ρ(r) to d(r) or to Λ(r). The idea behind our choice is three-fold: (i) equations(6.17), (6.18), and (6.19) are valid as written in all spatial dimensions; (ii) at large distances Λ(r) has the structure of hydrodynamic phase correlations; (iii) the final expressions have the simplest form possible within the same accuracy. It is also worth mentioning that n qc is a quantity which controls all physical processes happening in the WIGB at short distances at low temperature, e.g. m-body recombination rates. In all spatial dimensions it plays the same role as the condensate density in 3D system as long as one is interested in length scales not much larger than the healing length [6,13]. Finally, the condensate density is defined from It is not accidental that the ultimate long-wave length property of the superfluid system is determined last. It was always an unpleasant feature of numerous meanfield treatments that the crucial parameter determining physics at short scales was linked to n 0 .

Normal region
As was mentioned in the Introduction, we assume that the temperature is not very high, so that the condition (1.2) is preserved and the quantity U remains the only parameter characterizing the interaction between the particles.

Thermodynamic functions
In the normal region, where there are no anomalous correlators, we have to deal only with the Dyson equation for the Green's function G in terms of self-energy Σ 11 . Within the leading-order approximation, Σ 11 = 2nU (in 2D this formula implies an adequate choice of the parameter ǫ 0 ≡ ǫ 0 (T ) for the effective interaction U , discussed in subsection 7.2), and the expression for G yields the self-consistent relation for the number density, For the pressure we find , and (7.5) specify n, µ, and p as functions of (T,μ). Utilizing (4.29), with x ≡μ, and then (4.30), we find for the entropy and energy Figure 11. Second-order contributions to the self-energy in the normal regime, including the sunrise diagram.

Effective interaction in the normal regime in 2D
In the quasi-condensate region, the parameter ǫ 0 defining the effective interaction U can vary over a wide range of values thanks to the retained second-order correction that produces counter-terms compensating for the arbitrary choice. Since in the normal region we confine ourselves to the first-order expression for the self-energy, Σ = 2nU , we need to choose the value of ǫ 0 which minimizes the omitted second-order contributions originating from the sunrise diagram, see figure 11. Due to momentum independence of the pseudo-potential line, the third (fourth) diagram is identical to the first (second) one; we thus consider only the first two diagrams and multiply the result by a factor of 2. The parameter ǫ 0 can be much smaller or much larger than T . In either case the leading term comes from the logarithmic ultraviolet contribution from the product of two propagators, yielding the result where k T is the thermal momentum. This leads to the expression Comparing the first two terms in the brackets with (3.51) and (3.52), we see that if we take the value of ǫ 0 substantially away from T , then the leading second-order correction to the expression Σ = 2nU amounts to renormalizing the value of the pseudo-potential U in such a way that it corresponds to ǫ 0 ∼ T . This brings us to the conclusion that ǫ 0 ∼ T is the optimal choice for ǫ 0 , in which case omitting the diagrams shown in figure 11 is justified by the parameter mU ≪ 1.

Expansion parameter
At temperatures above the fluctuation region where the renormalized chemical potential remains small |μ| ≪ T the expansion parameter is defined by the infrared behaviour of the zero-frequency lines. The dimensionless factor, γ T , associated with adding an extra full interaction vertex to a diagram is given by (5.8), but now with G ∼ T /ǫ(k), resulting in the estimate whereμ fluc is the size of the fluctuation region in terms of the chemical potential. At temperatures much higher than the degeneracy temperature T c ∼ n 2/d /m wherẽ µ ≈ −T ln(nλ d T ), all corrections are small in the parameter (7.11)

Pseudo-Hamiltonian
In Section 4 we have obtained the thermodynamics of the system in the quasicondensate region, see (4.17)-(4.18), (4.28), (4.31), and (4.32) which specify n, µ, p, s, and ε as functions of two independent variables,μ and T , in the form of integrals involving theμ-dependent quasiparticle spectrum, E(k), given by (4.6). So far, we did not discuss the physical meaning of the function E(k), since this was not necessary for our derivations. On the other hand, equation (4.31) for the entropy clearly suggests that, within our approximation, the system is equivalent to a system of non-interacting bosonic quasiparticles with dispersion E(k), described by the Hamiltoniañ  It is important, however, to remember that E 0 is a function of bothμ and T , satisfying In view of this temperature dependence, we refer to (8.1) as a pseudo-Hamiltonian, rather than an effective Hamiltonian. The pseudo-Hamiltonian description can be applied to the normal region as well. In this case, were n ≡ n(μ, T ) is specified by (7.1). Theμ-dependent spectrum of quasi-particles is now given by (7.2).

Fluctuation region
In the fluctuation region the long-wave part of the classical-field component of the quantum field experiences strong (non-perturbative) fluctuations. It is exclusively due to these fluctuations that the diagrammatic technique for the quantum field loses an expansion parameter. The special role of the classical-field component is evident from the fact that on approach to the fluctuation region the leading contribution to the expansion parameter γ T is associated with zero-frequency propagators, see subsections 5.2 and 7.3. The diagrammatic technique built on zero-frequency propagators (with perturbative parts of self-energies included) and the interaction represented by pseudopotential U corresponds to the diagrammatic expansion of the Gibbs distribution for the in momentum space truncated classical field with the Hamiltonian functional Here k ′ ≪ k T = √ mT is a truncation momentum (to be discussed later), and µ ′ ≡ µ ′ (k ′ ) is the reduced chemical potential obtained from the bare one by subtracting relevant perturbative contributions from all the harmonics with k > k ′ , including the quantum ones. The bottom line is that describing the quantum gas in the fluctuation region reduces to solving the classical-field problem (9.1)-(9.3), in its fluctuation region.
In d = 2, 3, numerical solutions to the problem (9.1)-(9.3)-in terms of scaling functions for thermodynamic quantities-are available in the literature and have already been used to accurately describe weakly interacting 2D and 3D quantum gases in the fluctuation region [15,16,17]. The same is possible in 1D, but to the best of our knowledge it has not been done so far.
In the fluctuation region, an important quantity is the momentumk ≪ k T separating weakly and strongly coupled classical modes. For a quantitatively accurate description of the fluctuating classical-field sub-system, the momentum k ′ separating classical modes of interest from the rest of the modes has to be much larger thank. However, as long as we are interested in the order-of-magnitude estimates, it is safe and convenient to set k ′ ∼k, so that all the modes we are dealing with in (9.1)-(9.3) are strongly coupled. To estimatek we note that with k ′ ∼k all three terms in the Hamiltonian (9.3) have to be of the same order of magnitude, that is is the long-wavelength contribution to the total density, andμ ≡ µ ′ (k ′ →k). Forñ we haveñ ∼k d nk, and sincek is separating strongly coupled long-wave harmonics from slightly perturbed short-wave ones, the order-of-magnitude estimate for nk follows by continuity from the ideal system formula: Substituting this back into (9.4)-(9.5) yields The estimates (9.8)-(9.9) imply the following parameterization of the grand-canonical equation of state in the fluctuation region and its vicinity: (9.11) Here λ (d) (x) is a dimensionless scaling function of a dimensionless scaling variable x.
In 2D and 3D, the quantities n Similarly, superfluid and condensate densities-in the dimensions in which they are meaningful-are parameterized as and The functions λ(x), f s (x), and f 0 (x) are universal for all weakly interacting U(1)symmetric systems in the given dimension (no matter quantum or classical, continuous space or lattice); the numerical data for them is available for d = 2, 3 in [15,17]. By numerically solving the problem (9.1)-(9.3), it was established that, up to higher-order corrections in the parameter γ 0 , [15,16,17] [15,16,17] Note that the leading terms in the above relations, as well as order-of-magnitude estimates for the values of sub-leading terms, are readily obtained from the condition γ T ∼ 1; the accurate numerical treatment of the problem (9.1)-(9.3) is necessary to fix the values of the dimensionless constants.
By re-writing (9.15), (9.14) in terms of the critical temperature as a function of density, we get  where T (0) c is the critical temperature of the ideal 3D gas. Now it is easy to estimatẽ k. The relevant temperature for having strong fluctuations is given by The estimate (9.7) at T = T fluct then yields We see thatk ≪ k T fluct in all cases, which justifies the statement that the physics of the fluctuation region is that of a (strongly interacting) classical field.
Putting aside n s and n 0 , which are most sensitive to non-perturbative fluctuations of the classical field ψ, the next two quantities which are sensitive to fluctuations, especially in 1D and 2D, are the density and the chemical potential-see estimates (5.20). For other quantities-as is seen from (5.21)-neglecting fluctuation corrections will not result in a significant error. It is also important that up to a few dimensionless constants characterizing sub-leading contributions to critical values of p, ε, and s, the fluctuation corrections to these quantities are expressed in terms of the same function λ (d) (x) and its integral. Indeed, equation (4.19) implies the following relation for p: and then with (4.29)-(4.30) one obtains similar relations for s and ε.

Comparison with numerical results
To see how accurate our description is for the density matrix we performed Monte carlo simulations of weakly interacting two-and one-dimensional systems where phase fluctuations are the strongest. We deliberately aim at systems with substantial decay of the density matrix to see the difference between the perturbative treatment and exact solution. In figure 12 we present data for the two dimensional square-lattice system described by the Bose-Hubbard Hamiltonian where t is the hopping amplitude between the nearest neighbor sites. Model (10.1) was simulated for U/t = 0.25, T /t = 0.5 and filling factor n = 0.5 on a lattice with L × L = 200 × 200 points by using the Worm Algorithm approach [31]. Formally, all expressions in this manuscript relating energy spectrum, particle density, and density matrix remain valid for an arbitrary dispersion relation with ǫ(k) → k 2 /2m at small momenta. Moreover, the theory is supposed to work as written for systems with periodic boundary conditions provided the sums over momenta are understood as L −d k =0 , where k = 2πn/L and n is a d-dimensional integer. To suppress statistical noise, simulation data for ρ(r) were collected to spherically symmetric bins. We are using exactly the same procedure to present theoretical data. Thus, the only difference between the essentially exact (up to statistical errors) simulation data and the theory is due to the finite value of U . Even for relatively large values of U and temperature T ∼ 0.6T c we observe a remarkable accuracy of (6.13). One can get an idea of the systematic theoretical error by comparing curves derived from (6.13) and (6.17). [The latter is obtained by exponentiating the phase correlator at large distances].
In figure 13 we apply our theory to the Lieb-Liniger model of the one-dimensional Bose gas with contact interactions. Since exact correlation functions at finite temperature are not known, we performed Monte carlo simulations using recently developed techniques for continuous space systems [32,33]. It is convenient to introduce the on-dimensional characteristic length l 0 = 2/mU as the unit of length. The limit of weak interactions is obtained then by considering large densities na ≫ 1, and we have chosen nl 0 = 8 corresponding to the power law-decay of off-diagonal correlations at zero temperature with the exponent α = 1/π √ 2nl 0 = 1/4π. The upper curve in figure 13 shows results for low temperature mT /(2πn 2 ) = 0.001, or T /nU ≈ 0.025 when thermal phase fluctuations are barely influencing the data for the simulated system size L/l 0 = 25. When temperature is increased to mT /(2πn 2 ) = 0.004 (lower curve in figure 13) we clearly see the bimodal decay of the density matrix and the crossover from the power-law to exponential decay. Contrary to the 2D case, equation (6.17) captures the actual behaviour better than (6.13) which is not so surprising given the large suppression of ρ(r) at large distances.
Finally, figures 17 and 18 show the agreement of Monte Carlo results for the energy (squares) and the pressure (triangles) at temperatures T > T c with theoretical curves from (7.7) and (7.5) (solid lines).
In the following we compare thermodynamic functions obtained in Sections 4 and 7 with exact quantum Monte Carlo results of a weakly interacting three-dimensional system. Let us start with the case of a homogeneous system with a small parameter na 3 = 10 −6 . In figures 14 and 15 we compare Monte Carlo results for the energy (squares) per particle and pressure (triangles) at temperatures T < T c , with the theoretical curves from (4.32) and (4.28) (solid lines). The analytical expressions are in very good agreement with numerical results up to temperatures T ∼ T c and there is no need to switch to the description in terms of universal functions [15]. This is consistent with the estimate (5.21) showing that fluctuation contributions to energy and pressure are very small. The same does not apply to the superfluid and condensate densities for which the perturabtive expansion in U is not valid as fluctuations of the order parameter cannot be neglected on approach to the critical temperature. In  Quantum Monte Carlo results for superfluid fraction (squares) and condensate fraction (triangles) for the three-dimensional Bose gas with na 3 = 10 −6 . Solid and dashed lines are a combination of theoretical expressions from (4.59) for the superfluid density and (4.7) for condensate density with preexisting classical Monte Carlo results [15] for the fluctuation region (the arrow is pointing to the gap in the curves where the two results are matched).  figure (16) squares and triangles are the results of quantum Monte Carlo simulations for the superfluid and condensate fraction, respectively. Solid and dashed lines are the theoretical expressions from (4.59) for the superfluid density and (4.7) for the condensate density, respectively, combined with preexisting classical Monte Carlo results [15] for the fluctuation region. At a temperature T 0.6T c we switch to the universal description, as γ (full) T [see (5.12)] becomes of the order of unity. For T ∼ T c the remaining discrepancy between classical and quantum Monte Carlo results is due to finite interaction strength. One should not be misled by the small value of the gas parameter in the fluctuation region since the proper parameter controlling the size of the fluctuation region in temperature, ∆T /T c , and the accuracy of the universal description involves a large numerical prefactor and is rather na 3 /10 −5 , see [15].

Conclusions
We have shown that Beliaev's diagrammatic technique regularized by adding a small term explicitly breaking the U(1) symmetry to the Hamiltonian yields a simple and controllable description of the weakly interacting Bose gas in any dimension. This approach is especially convenient for obtaining thermodynamic functions that are not sensitive (up to higher-order corrections) to the long-range phase fluctuations. The symmetry breaking term introduces a small gap in the otherwise gapless Goldstone mode and suppresses long-range fluctuations of the phase of the order parameter, thereby introducing a genuine condensate density and thus substantially simplifying the description of lower-dimensional systems. In this sense, the symmetrybreaking trick is an alternative to Berezinskii's finite-size trick and Popov's special (hydrodynamic) treatment of long-wave parts of the fields. At the end of the calculation, the symmetry breaking term is set to zero.
Another important feature of our approach is that it completely avoids such notions as seed-or quasi-condensate for its construction. The quasi-condensate density defined as the modulus squared of the order parameter field at distances larger then the healing length is calculated as the long-wave property of the theory, rather than serving as an input parameter. It is only natural to have the theory which handles the short-range physics first and uses it for dealing with the low-energy physics next. As far as long-range off-diagonal correlations are concerned, the approach produces accurate answers for the off-diagonal correlation functions up to distances where further evolution of the correlators is controlled by generic hydrodynamic relations, and thus can be accurately extrapolated to arbitrarily large distances.
We confined our analysis to the most typical and non-trivial (in 2D and 3D) case of dilute gas, when the size of the potential R 0 is much smaller than the distance between the particles. In the (Boltzmann) high-temperature regime, we have restricted ourselves to low enough temperatures at which de Broglie wavelength remains much larger than R 0 and the answers do not depend on the details of the pair potential. Under the two above-mentioned conditions, the effective inter-particle interaction is described by a single parameter U . Technically, generalization of the theory to the case of R 0 n −1/d is straightforward. In this case, the weakness of interaction literally means Born interaction potential, and complete summation of ladder diagrams in 2D and 3D becomes irrelevant. The description of a normal gas at λ T n −1/d is a well-studied topic that goes beyond the scope of the present paper.
The pair of variablesμ = µ − 2nU and T form the most convenient set of parameters determining the rest of thermodynamic quantities, including the condensate density, if any. The structure of the answers corresponds to a picture of non-interacting bosonic quasiparticles with the Bogoliubov spectrum (4.6), on top of the 'vacuum' with temperature-dependent energy. This allows one to introduce the bilinear bosonic pseudo-Hamiltonian (8.1) with Bogoliubov spectrum. We argue that within the bilinear pseudo-Hamiltonian ansatz our results cannot be further improved at least at T nU . In the T ≪ nU limit, when thermodynamics is exhausted by dilute non-interacting phonons, using the improved value of the sound velocity obtained directly from the zero-temperature compressibility, yields more accurate answers for thermodynamic quantities, but is inconsistent with retaining the Bogoliubov form of the spectrum for all the momenta. We also show that any attempt to improve the self-consistent mean-field description of the superfluid region, aimed at getting rid of the spurious first-order phase transition, is senseless as long as the fluctuation-induced non-perturbative shift of the critical temperature is not accounted for. In fact, there is no need to extrapolate the pseudo-Hamiltonian description to the fluctuation region around the critical temperature in view of the availability of the numerical answers for universal scaling functions describing all the weakly interacting U(1) systems at |T − T c |/T c ≪ 1. With these data, the two analytical descriptions-in the normal and superfluid phases-are readily connected across the fluctuation region. By comparing our results to first-principles Mote Carlo data, we observe that the most vulnerable region of parameters is |T − T c |/T c ≪ 1, where the description is accurate only when the gas parameter an 1/3 is as small as ∼ 10 −2 .