Condensation and superfluidity of dilute Bose gases with finite-range interaction

We investigate an ultracold and dilute Bose gas by taking into account a finite-range two-body interaction. The coupling constants of the resulting Lagrangian density are related to measurable scattering parameters by following the effective-field-theory approach. A perturbative scheme is then developed up to the Gaussian level, where both quantum and thermal fluctuations are crucially affected by finite-range corrections. In particular, the relation between spontaneous symmetry breaking and the onset of superfluidity is emphasized by recovering the renowned Landau's equation for the superfluid density in terms of the condensate one.


Introduction
Since the seminal investigation on liquid Helium by H. Kamerlingh Onnes, the research on low-temperature physics has been focused on the understanding and engineering of superfluid states of matter [1]. Nowadays, superfluids and superconducting materials are at the core of a new technological revolution centered around the development of quantum devices [2,3,4].
From a theoretical point of view, an enormous effort has been devoted to finally provide a microscopic theory accounting for the transition between a normal state and a superfluid one, where dissipationless flow occurs. This research oscillated between the necessity of statistical mechanics to identify general mechanisms leading to supertransport and the interests of condensed matter theorists and material scientists on peculiar setups with unique properties. Concerning metallic superconductors, the theoretical investigation reached one of its peaks with the BCS theory and its consequent refinements [5,6].
Moving to the superfluid side, liquid Helium remains for decades the only efficient platform to probe the markers of superfluid behaviour, like the absence of viscous forces and the vorticity quantization [7,8]. The formulation of a microscopic theory for the superfluid phase of liquid Helium proves to be an exceptionally demanding task. Up to now, a reliable picture can be achieved via ab initio numerical simulations such as Path-Integral Monte Carlo algorithms [9]. Indeed, liquid Helium has to be classified as a strongly-correlated system, due to the experimental values at play for density and interaction strength. In this situation, even identifying a smallness parameter is non trivial, preventing the effective implementation of a perturbative expansion [10].
The two-fluid phenomenological approach by Landau [11,12] has been much more fruitful, enabling the derivation of a self-consistent hydrodynamic theory over few crucial assumptions. Remarkably, the original formulation of Landau did not rely upon an atomic point of view and, moreover, did not invoke any symmetry breaking related to the modern characterization of phase transitions [13].
Since the pioneering guess of London in 1938 [14,15], it has been a common approach to interpret the superfluid transition in terms of Bose-Einstein condensation, where a macroscopic fraction of Helium atoms can be described by a macroscopic wavefunction. Unfortunately, analytical result can be obtained only in the weakly interacting limit, which does not hold up to the experimental values for Helium. Applied to Helium, the resulting picture is only qualitative since, for instance, it does not even manage to capture the peculiar rotonic minimum of the excitation spectrum.
In 1995, the experimental realization of a Bose-Einstein condensates [13] changed the scenario in a crucial way: for the first time, the predictions of the Bogoliubov theory [16] have been checked in actual weakly-interacting quantum gases. At the same time, within a field-theory approach, it is possible to recover the Landau main results moving from a microscopic Lagrangian for ultracold atoms.
The several successful theoretical studies based on the Bogoliubov framework move from the crucial assumption that the true atom-atom interaction can be replaced by a contact (i.e. zero-range) pseudopotential whose strength is given by the s-wave scattering length a s [17,18]. The resulting thermodynamics is universal since there is no dependence on the potential shape, with only a s playing a relevant role. The same point can be made for transport quantities as the superfluid fraction. Despite the many achievements of this strategy, current experiments deal with higher density setups, reduced dimensionalities and more complex interactions [19,20]. Thus, it is pressing to extend the usual two-body zero-range framework in order to capture more realistic and interesting experimental regimes. Within a functional integration formalism, atoms are represented by a bosonic field whose dynamics is governed by a microscopic interacting Lagrangian density. The coupling constants of the finite-range theory can be determined in terms of the s-wave scattering parameters, namely a s and the corresponding effective range r e . In [21,22,23,24], the finite-range thermodynamics is derived up to the Gaussian level for a three dimensional uniform Bose gas, while the non-trivial case of two spatial dimensions is addressed in [25,26]. In Fig. 1 we report a visual summary of the major analytical approaches to modelling bosonic quantum gases. A similar analysis concerning the superfluid properties of a finite-range theory is still missing and it is the main subject of this work. By adopting a functional integration point of view as in [23,25], we are going to show that both condensate and superfluid depletion are modified by the finite-range character of the two-body interaction. Moreover, they are not independent from each other but the familiar Landau equation for the superfluid density can be formulated in terms of the condensate one.
The paper is organized as follows: in the next section we discuss the Landau twofluid model from a field-theory perspective moving from a microscopic Langrangian density. The key point consists in properly accounting the different response of normal and superfluid components to a Galilean boost [27]. Technically, this corresponds to performing a phase twist on the order parameter [28,29], which has to be considered throughout the whole perturbative expansion. By following this scheme, we will derive a Landau-like equation relating superfluid and condensate density. Our calculation are carried on in a generic dimension D, simplifying the application of dimensional regularization to heal the UV-divergent zero-point energy. We then specify to the cases D = 3 and D = 2.

Two-fluid model of a superfluid
The thermodynamic properties of a physical system can be described calculating the grand canonical partition function Z, which is related to the grand potential Ω by We consider a uniform D-dimensional Bose gas of identical cold atoms described by the complex scalar field ψ( r, τ ). We calculate the grand canonical partition function Z as the functional integral where is the Euclidean action and β = 1/(k B T ), with k B the Boltzmann constant and T the absolute temperature. We introduce the non-relativistic Lagrangian where µ is the chemical potential and V (| r − r |) is the generic interaction between bosons, assuming it is dependent only on the relative distance | r − r |.
In the context of the two-fluid model of a superfluid [11,12] we describe the fluid behavior of the system as composed by a mixture of a normal component and a superfluid component. In particular, we consider the normal part as a fluid current moving with velocity v with respect to the laboratory frame of reference. Working with imaginary time, we describe this motion by substituting the time derivative in the Lagrangian (4) with the Lagrangian fluid derivative Moreover, since the superfluid part does not exchange momentum with the normal part, we impose a superfluid current with a phase twist of the bosonic field [28] ψ( r, τ ) → e i m vs· r h ψ( r, τ ) Substituting these expressions in the Lagrangian (4) we obtain where we define the effective chemical potential µ e as [27] Considering that in the condensed phase the U(1) global symmetry is spontaneously broken we use the bosonic field parametrization ψ( r, τ ) = ψ 0 + η( r, τ ) where η( r, τ ) is the complex field describing the fluctuation around the uniform field configuration ψ 0 , which represents the order parameter of the condensate phase transition. Substituting the parametrization (9) in the action (3) and keeping only quadratic terms in the fluctuation field η( r, τ ) we obtain the homogeneous system action and we calculate the Gaussian action in the Fourier space as [30] where k = ( k, ω n ) are the D + 1 wavevectors, ω n = 2πn βh are the bosonic Matsubara frequencies and with Performing the functional integration of the Gaussian action (11) we obtain the partition function (2) and the grand potential (1). In particular, we find that the grand potential Ω is the sum of two contributions where is the order parameter grand potential and is the Gaussian contribution to the grand potential, where is the excitation spectrum of the Bose gas. The Hugenholtz-Pines theorem is guaranteed imposing the saddle point condition ∂Ω/∂ψ 0 = 0, for which the excitation spectrum becomes gapless and the chemical potential reads Moreover, we identify the condensate density as n 0 = ψ 2 0 . Notice that, within the Gaussian approximation of the action, the excitation spectrum does not contain the anomalous density, which is instead included by adopting the next-next-to-leading Hartree-Fock-Bogoliubov scheme [31,32,33]. The initial assumption that the real space interaction depends only on the distance between bosons implies that the interaction potential is left unchanged by a reflection of the momenta:Ṽ ( k) =Ṽ (− k). Due to this property we are able to perform the summation over Matsubara frequencies ω n in the Gaussian grand potential (16) [27,34], obtaining the grand potential Ω in the form where Ω 0 is given by Eq. (15), is the zero-temperature Gaussian contribution, written as the sum of noninteracting elementary excitations with spectrum E k (ψ 2 0 ), and is the finite-temperature Gaussian contribution. The beyond-mean-field Gaussian equation of state with finite-range interaction has been analyzed in previous papers [21,23,25]. For the sake of completeness we report the main results in the Appendix. We now calculate the superfluid density n s of the system with a self-consistent approach, employing the Gaussian grand potential (19). In particular, we will identify n s from the calculation of the total momentum density P of the fluid, which is obtained as follows where we take the derivative of the grand potential with respect to the velocity − v first and then we substitute the mean field value of the chemical potential µ e = g 0 ψ 2 0 . We find that the momentum density P is given by where, since the grand potential (19) is given by the sum of three sum contributions, the three terms of the momentum density P are defined accordingly. In particular, we find The Gaussian zero-temperature contribution where we define the zero-temperature number density contribution f (0) g (n 0 ) as The Gaussian thermal contribution P (T ) g to the momentum density (23) is more involved: Assuming that the difference between the velocity v of the normal fluid and the velocity v s of the superfluid is small, we can expand the exponential and, taking into account that some terms are zero for symmetry reasons, we obtain where f (T ) g (n 0 ) is the Gaussian thermal density contribution and we have defined the normal density of the fluid as Notice that, in the thermodynamic limit L D → ∞, the normal fluid density (30) is fully consistent with the familiar Landau result [12]. In conclusion, putting together the contributions (24), (25) and (28) we rewrite the momentum density P as In the square bracket we identify the number density n(n 0 , T ), expressed as a function of the condensate number density n 0 and the temperature T as follows [27] n(n 0 , We remark that we take the derivative of the grand potential with respect to the chemical potential first and then we substitute the mean field value of the chemical potential µ e = g 0 ψ 2 0 , with the identification for the condensate density n 0 = ψ 2 0 . This procedure can be justified considering that the same procedure is implemented to calculate the condensate fraction of a noninteracting Bose gas [35]. With this identification, we express the momentum density P as Finally, we identify the superfluid density n s as which allows us to express the momentum density P as the sum of the momentum density n s m v s of the superfluid part of the fluid and the momentum density n n (n 0 , T ) m v of the normal part of the fluid, namely We emphasize that Eq. (34) constitutes the main result of this paper, since it highlights the non-trivial relationship between the superfluid density n s , the condensate density n 0 and the temperature T . This result may be regarded as the explicit formulation, at a Gaussian level, of the Josephson relation [36]. Notice that our number density (32) and the superfluid fraction (34) are equivalent to the result obtained with Beliaev diagrammatic technique in reference [37] if we approximate n 0 ≈ n in Eqs. (26), (29) and (30) and we consider the zero-range interactionṼ ( k) = g 0 . In the following section we implement the superfluid density calculation for bosons with finite-range interaction in three-and two-dimensional systems.

Superfluid density of bosons with finite-range interaction
In order to obtain explicit formulas for the superfluid density n s , in this section we shall implement Eq. (34) in the three-and in the two-dimensional Bose gas, considering the explicit form of the interaction V ( r). The usual approximation to study a weaklyinteracting Bose gas of ultracold atoms is constituted by the zero-range interaction V ( r) = g 0 δ( r), which in the Fourier space gives Here we improve this approximation, considering the finite-range effective interactioñ which is obtained adding to the zero-range interaction strength g 0 the first nonzero correction in the gradient expansion of the real interaction potential V (| r − r |), namely g 2 k 2 , where we define In the three dimensional case, the values of the couplings g 0 and g 2 are determined with the scattering theory in terms of the s-wave scattering length a s and the effective range r e as follows [23] g 0 = 4πh 2 a s m , In the two-dimensional case we choose the zero-range interaction coupling g 0 according to reference [38] and we derive the coupling g 2 from the definition of the characteristic range R = 2 |g 2 /g 0 | discussed in [25], obtaining We now explicitly implement the superfluid density n s calculation for the finiterange effective interaction (37). According to Eq. (34), n s is obtained by subtracting the normal density n n (n 0 , T ) from the number density n(n 0 , T ). We first calculate the number density n(n 0 , T ), which is given as the sum of the three contributions of Eq. (32). The first contribution is the condensate density n 0 . The second contribution is the zero-temperature Gaussian contribution to the normal density, namely where S D = 2π D/2 /Γ[D/2] is the D-dimensional solid angle, in which Γ[x] is the Euler gamma function and where we define the excitation spectrum with λ(g 2 , n 0 ) = 1 + 4m h 2 g 2 n 0λ (g 2 , n 0 ) = 1 + 2m h 2 g 2 n 0 Notice that the excitation spectrum E k (n 0 ) reproduces the familiar Bogoliubov spectrum if the zero-range interaction is restored by putting g 2 = 0. Since f (0) g (n 0 ) is ultraviolet divergent, we will regularize it using dimensional regularization [39,40]. In particular, we obtain an adimensional integral using the integration variable t = h 2 k 2 λ(g 2 , n 0 )/(4mg 0 n 0 ), then we extend the spatial dimension D to the complex value D = D − ε. We remark that this additional step is needed because the dimensional regularization procedure is not always able to heal the ultraviolet divergence of the integrals [35]. After the integration, we obtain f (0) g (n 0 ) in the form where B(x, y) is the Euler Beta function and κ is an ultraviolet cutoff wavevector introduced for dimensional reasons.
The third term in the number density of Eq. (32) is the finite-temperature Gaussian contribution f (T ) g (n 0 ), which, unlike the zero-temperature one, is convergent. However, this integral can be calculated analitically only in the low-temperature regime, where it is useful to introduce the integration variable x = βE k (n 0 ). Doing so in Eq. (29) in which the finite-range interaction (37) is substituted, we get Before substituting the spatial dimension D to calculate explicitly the number density contributions obtained in the previous equations, let us also calculate the normal density n n (n 0 , T ), which is given by Eq. (30) where the finite-range interaction (37) is substituted. In analogy with the finite-temperature density contribution f (T ) g (n 0 ), the normal density n n (n 0 , T ) can be calculated analytically only in the low-temperature regime: as before we introduce the integration variable x = βE k (n 0 ), obtaining where k(x) is given again by Eq. (46). We now calculate explicitly the number density n(n 0 , T ) and the superfluid density n s (n 0 , T ) in D = 3 and in D = 2

D = 3
As a preliminar result for the superfluid density n s , we employ Eq. (32) to calculate the density n (n 0 , T ) of bosons with finite-range interaction. In D = 3, the zerotemperature Gaussian contribution f (0) g (n 0 ) is regularized simply by the means of dimensional regularization, therefore we put ε = 0 into the equation (44). In the limit of small g 2 we get This result allows us to calculate the zero-temperature number density n(n 0 , T = 0) as the sum of the condensate density n 0 and the zero-temperature density contribution f (0) g (n 0 ). In particular, we substitute explicitly the expressions of the couplings g 0 and g 2 as given by Eq. (39), obtaining n(n 0 , T = 0) as a function of the three-dimensional s-wave scattering length a s and the effective range r e : n(n 0 , T = 0) = n 0 1 + which differs from the result of reference [22] by a factor 2 in the finite-range correction.
In Fig. 2 we compare our finite-range condensate fraction n 0 /n (black dot-dashed line) obtained from the numerical solution of Eq. (49) with the zero-range result (blue solid line) obtained putting r e = 0 and the classical result of Bogoliubov (red dashed line) [16]. Let us also calculate the finite-temperature contribution f (T ) g (n 0 ) in D = 3, given by the integration of Eq. (45) in the low-temperature regime The finite-temperature number density n(n 0 , T ), according to Eq. (32), is given by the sum of the condensate density n 0 and the Gaussian density contributions of Eqs. (48) and (50) n (n 0 , T ) = n 0 + 1 3π 2 We also rewrite the general form of n(n 0 , T ) in terms of the three-dimensional gas parameter na 3 s and the effective range r e , employing the explicit form of the couplings g 0 and g 2 given by Eq. (39), namely This equation can be used to express the condensate fraction n 0 /n explicitly in the very weakly-interacting regime in which a s , r e → 0, where one can approximate in the second and third subleading terms the condensate density n 0 with the density n, since the phenomenon of quantum depletion is absent in the noninteracting zero-temperature limit and these terms are finer corrections with respect to the first. We obtain where we have rescaled the temperature in terms of T BEC = (2πh 2 n 2/3 )/(mk B ζ(3/2) 2/3 ) for noninteracting bosons, and ζ(x) is the Riemann zeta function. We emphasize that, at T = 0 and for a zero-range interaction for which r e = 0, the Bogoliubov result for the condensate quantum depletion is reproduced [16]. Finally, we calculate the normal density n n of Eq. (47) substituting D = 3 and considering the low-temperature regime, in which we get The superfluid density n s for a system of bosons interacting with the finite-range interaction is obtained substituting Eqs. (49), (50) and (54) in Eq. (34), namely As for the condensate fraction, we obtain an explicit expression of the superfluid fraction n s /n in the very weakly-interacting limit in which one can approximate the condensate density n 0 with the number density n in the subleading terms, rewriting the previous equation as where we substitute also the expressions for g 0 and g 2 of Eq. (39). Notice that, while the superfluid fraction n s /n = 1 at zero temperature, the condensate fraction is n 0 /n < 1 due to the quantum depletion. Obviously, Eq. (56) is reliable only in the deep T /T BEC 1 regime but, more generally, one must consider the full solution of Eq. (47).
In Fig. 3 we report the superfluid fraction n s /n as a function of the scaled temperature T /T BEC for three values of the ratio r e /a s at fixed gas parameter na 3 s by numerically solving Eq. (47). Fig. 3 shows that a positive r e /a s slightly enhances the superfluid fraction while the opposite occurs for negative values.

D = 2
Here we formally follow the same path we have introduced in the three-dimensional case, obtaining the number density n at zero-temperature and employing it to calculate the superfluid density n s . In the two-dimensional case we calculate the regularized zerotemperature density contribution f (0) g (n 0 ) from Eq. (44) where D = 2 is substituted. Notice that f (0) g (n 0 ) is obtained as a sum of the two terms inside the square bracket of (44): while the first term is finite, the divergence of the second term is healed by performing a Taylor expansion around = 0 and deleting the o( −1 ) divergence [41,42]. The regularized zero-temperature density contribution reads where we identify the ultraviolet energy scale 0 as If the finite-range interaction strength g 2 , as we suppose, constitutes a small correction of the zero-range term g 0 , one can expand the previous equation for small values of the adimensional parameter 2m h 2 g 2 n 0 , thus In this limit, the zero-temperature density n (n 0 , T = 0) reads n (n 0 , T = 0) = n 0 + 1 4π In analogy with the three-dimensional case, we derive an explicit -but approximated -formula for the zero-temperature condensate fraction n 0 /n in D = 2 considering the very weakly-interacting regime in which g 0 , g 2 → 0. As before, in this regime we can approximate n 2 0 ≈ n 0 n in the third term of Eq. (60) and n 0 ≈ n inside the logarithm, since it is a subleading term with respect to the first one. Let us come back to ε 0 , by chosing the cutoff as κ = 2π/a s . Then, by taking g 0 and g 2 as in Eq. (40), we obtain the approximated condensate fraction n 0 /n as a function of the two-dimensional gas parameter na 2 s and the characteristic range R, namely Notice that the zero-range interaction result by Schick for the condensate fraction n 0 /n [43] is easily reproduced by setting R = 0. Working outside the very weakly-interacting limit, one can also obtain the zero temperature condensate fraction from the numerical solution of Eq. (60). We report it as the black dot-dashed line in Fig. 4, in comparison with our zero-range result (blue solid line) and the result by Schick (red dashed line), which is reproduced in the weakly-interacting regime in which na 3 s 1. Finally, following the three-dimensional case, we may want to calculate also the finitetemperature density contribution f (T ) g (n 0 ). However, substituting D = 2 in Eq. (45) we find that is infrared divergent, therefore cannot be regularized with dimensional regularization. This result is indeed correct and reflects the absence of Bose-Einstein condensation at finite-temperature in two-dimensional systems [44], as already pointed out in [45]. The two-dimensional normal density n n (n 0 , T = 0) of bosons with finite-range interaction is obtained from the integration of Eq. (47), in which we expand the integrand in the low-temperature limit, namely Since the thermal contribution to the density f (T ) g (n 0 ) is divergent in D = 2, we cannot express the superfluid density n s as a function of the condensate density n 0 at a finite temperature T . Therefore we employ the zero-temperature number density of Eq. (60) to obtain the condensate density in the implicit form n 0 = n 0 (n, T = 0), namely as a function of the density n. In this way we calculate the superfluid density substituting Eq. (63) into Eq. (34) which we expect to provide a reliable approximation in the low-temperature regime, remembering that -outside it -our approach would in any case be incorrect due to the Berezinski-Kosterlitz-Thouless transition [46]. Finally, we calculate an approximated expression of the superfluid fraction n s /n in the very weakly-interacting regime where g 0 , g 2 → 0, in which we can approximate the condensate density inside Eq. (64) as n 0 (n, T = 0) ≈ n. In the context of this approximation, we also substitute in Eq. (64) the explicit form of the parameter λ(g 2 , n 0 ) of Eq. (43). Moreover, we remember that the two-dimensional interaction strengths g 0 and g 2 are given by Eq. (40), obtained in terms of the s-wave scattering length a s and the characteristic range R of the interatomic potential with the scattering theory. The approximated superfluid fraction n s /n reads where we rescale the result in terms of the temperature T * =h 2 n/(mk B ) of quantum degeneracy.

Conclusions
We have used finite-temperature one-loop functional integration to reproduce the density momentum equation of the two-fluid model. An analytical relationship between the density n and the condensate density n 0 has been obtained at zero-temperature and in the low-temperature limit. This result has been used to express the low-temperature superfluid density n s as a function of n 0 and T for bosons with finite-range interaction, which can be regarded as an explicit implementation of the Josephson relation. We analyze thoroughly the D = 3 and D = 2 case, but our approach could be applied also in D = 1, where we expect to reproduce the Lieb-Liniger theory except in the strong coupling regime [47].
We expect that our theory is meaningful under the conditions of diluteness of the bosonic gas, for which na D s 1 and nr D e 1. Moreover, finite-temperature results must be considered in the limit k B T /(g 0 n) 1, for which the mean thermal energy is much lower than the gas healing length. Our finite-range corrections to the condensate fraction n 0 /n and to the superfluid fraction n s /n can be detected in D = 3 in the regime a s /r e ≤ 1 and in D = 2 for a s /R ≤ 1 but not where they are much lower than 1. In that case, the higher order terms which we are neglecting in the gradient expansion of the interaction potential (37) become relevant. Notice that, in D = 3 and in D = 2, these may represent different regimes. In fact, while the effective range value r e can be tuned by the means of a Feshbach resonance, the characteristice range R is fixed, being essentially a geometric property of the real two-body interaction potential between the atoms. Indeed, we expect that R is proportional to the Van der Waals radius of the atoms and it can be numerically computed, following its definition, using a model two-body potential V ( r).
An extension of this work consists in the numerical calculation of the thermal integrals f (T ) g (n 0 ) and n n (n 0 , T ) outside the zero-temperature limit in D = 2, which we have considered for obtaining analytical results. In any case, we expect that our predictions fail to describe the superfluid fraction of the two-dimensional Bose gas at sufficiently high temperature, due to the occurrence of the BKT transition. In this case it is needed a profound rethinking of our approach, including explicitly in the bosonic field parametrization the contribution of vortex configurations of the phase field, which cause the BKT topological phase transion.
In the three-dimensional case, it is sufficient to substitute the dimension D = 3 and put ε = 0 in Eq. (67), to get the regularized zero-temperature Gaussian grand potential as This equation of state corrects the one obtained in Ref. [25]: here λ(g 2 , µ e /g 0 ) appears also inside the logarithm. Moreover, it is important to stress that, in the case of zerorange interaction, Mora and Castin [48] were able to extend Eq. (70) by including a beyond-Gaussian term. This next-next-to-leading extension in the finite-range case is highly non trivial and it surely deserves a separate detailed investigation.