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

A. Tononi  Dipartimento di Fisica e Astronomia “Galileo Galilei” and CNISM,
Università di Padova, Via Marzolo 8, I-35131 Padova, Italy  A. Cappellaro
Dipartimento di Fisica e Astronomia “Galileo Galilei” and CNISM, Università di
Padova, Via Marzolo 8, I-35131 Padova, Italy  L. Salasnich  Dipartimento di
Fisica e Astronomia “Galileo Galilei” and CNISM, Università di Padova, Via
Marzolo 8, I-35131 Padova, Italy  Istituto Nazionale di Ottica (INO) del
Consiglio Nazionale delle Ricerche (CNR), Via Nello Carrara 1, I-50019 Sesto
Fiorentino, Italy

###### Abstract

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.

Keywords: Bose-Einstein condensation, superfluidity, finite-range interaction,
quantum field theory, ultracold atoms.

##  1 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.

Figure 1: Flowchart diagram showing some the key theoretical developments
related to our analysis of three- and two-dimensional Bose gases with a
finite-range interaction potential.

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
two-fluid 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$.

##  2 Two-fluid model of a superfluid

The thermodynamic properties of a physical system can be described calculating
the grand canonical partition function ${\cal Z}$, which is related to the
grand potential $\Omega$ by

$\Omega=-\frac{1}{\beta}\ln({\cal Z}).$ (1)

We consider a uniform $D$-dimensional Bose gas of identical cold atoms
described by the complex scalar field $\psi(\vec{r},\tau)$. We calculate the
grand canonical partition function ${\cal Z}$ as the functional integral

${\cal Z}=\int{{\cal
D}[\bar{\psi},\psi]\;e^{-\frac{S[\bar{\psi},\psi]}{\hbar}}},$ (2)

where

$S[\bar{\psi},\psi]=\int_{0}^{\beta\hbar}d\tau\int_{L^{D}}d^{D}r\;{\cal
L}(\bar{\psi},\psi)$ (3)

is the Euclidean action and $\beta=1/(k_{B}T)$, with $k_{B}$ the Boltzmann
constant and $T$ the absolute temperature. We introduce the non-relativistic
Lagrangian

$\displaystyle{\cal L}=$
$\displaystyle\bar{\psi}(\vec{r},\tau)\bigg{(}\hbar\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu\bigg{)}\psi(\vec{r},\tau)$
(4) $\displaystyle+\frac{1}{2}\int d^{D}r\,^{\prime}\
|\psi(\vec{r}\,^{\prime},\tau)|^{2}\ V(|\vec{r}-\vec{r}\,^{\prime}|)\
|\psi(\vec{r},\tau)|^{2},$

where $\mu$ is the chemical potential and $V(|\vec{r}-\vec{r}\,^{\prime}|)$ is
the generic interaction between bosons, assuming it is dependent only on the
relative distance $|\vec{r}-\vec{r}\,^{\prime}|$.

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 $\vec{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

${\partial_{\tau}}\rightarrow{\partial_{\tau}}-i\vec{v}\cdot\vec{\nabla}$ (5)

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]

$\psi(\vec{r},\tau)\rightarrow
e^{i\frac{m\vec{v}_{s}\cdot\vec{r}}{\hbar}}\psi(\vec{r},\tau)$ (6)

Substituting these expressions in the Lagrangian (4) we obtain

$\displaystyle{\cal
L}=\bar{\psi}(\vec{r},\tau)\bigg{(}\hbar\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu_{e}+(\vec{v}-\vec{v}_{s})\cdot(-i\hbar\vec{\nabla})\bigg{)}\psi(\vec{r},\tau)$
$\displaystyle+\frac{1}{2}\int d^{D}r\,^{\prime}\
|\psi(\vec{r}\,^{\prime},\tau)|^{2}\ V(|\vec{r}-\vec{r}\,^{\prime}|)\
|\psi(\vec{r},\tau)|^{2},$ (7)

where we define the effective chemical potential $\mu_{e}$ as [27]

$\mu_{e}=\mu-\frac{1}{2}m\vec{v}_{s}\cdot(\vec{v}_{s}-2\vec{v})$ (8)

Considering that in the condensed phase the $\mbox{U}(1)$ global symmetry is
spontaneously broken we use the bosonic field parametrization

$\psi(\vec{r},\tau)=\psi_{0}+\eta(\vec{r},\tau)$ (9)

where $\eta(\vec{r},\tau)$ is the complex field describing the fluctuation
around the uniform field configuration $\psi_{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 $\eta(\vec{r},\tau)$ we obtain the homogeneous system action

$S_{0}=\beta\hbar L^{D}(-\mu_{e}\psi_{0}^{2}+\frac{1}{2}g_{0}\psi_{0}^{4})$
(10)

and we calculate the Gaussian action in the Fourier space as [30]

$S_{g}[\bar{\eta},\eta]=\frac{\hbar}{2}\sum_{k}\left(\bar{\eta}(k)\
\eta(-k)\right)\ {\cal M}\ \left(\begin{array}[]{c}\eta(k)\\\
\bar{\eta}(-k)\end{array}\right)$ (11)

where $k=(\vec{k},\omega_{n})$ are the $D+1$ wavevectors,
$\omega_{n}=\frac{2\pi n}{\beta\hbar}$ are the bosonic Matsubara frequencies
and

${\cal M}=\left(\begin{array}[]{c c }{\cal M}_{11}&{\cal M}_{12}\\\ {\cal
M}_{21}&{\cal M}_{22}\\\ \end{array}\right)$ (12)

with

$\displaystyle{\cal M}_{11}$
$\displaystyle=-i\hbar\omega_{n}+\frac{\hbar^{2}k^{2}}{2m}-\mu_{e}+g_{0}\psi_{0}^{2}+\psi_{0}^{2}\,{\tilde{V}}(\vec{k})+\hbar(\vec{v}-\vec{v}_{s})\cdot\vec{k}$
$\displaystyle{\cal M}_{22}$
$\displaystyle=+i\hbar\omega_{n}+\frac{\hbar^{2}k^{2}}{2m}-\mu_{e}+g_{0}\psi_{0}^{2}+\psi_{0}^{2}\,{\tilde{V}}(\vec{k})-\hbar(\vec{v}-\vec{v}_{s})\cdot\vec{k}$
$\displaystyle{\cal M}_{12}$ $\displaystyle={\cal
M}_{21}=\psi_{0}^{2}{\tilde{V}}(\vec{k})$ (13)

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 $\Omega$ is the sum of two contributions

$\Omega=\Omega_{0}+\Omega_{g}$ (14)

where

$\Omega_{0}=L^{D}(-\mu_{e}\psi_{0}^{2}+\frac{1}{2}g_{0}\psi_{0}^{4})$ (15)

is the order parameter grand potential and

$\Omega_{g}=\frac{1}{2\beta}\sum_{\vec{k},\omega_{n}}\ln[\beta^{2}(E_{\vec{k}}(\psi_{0}^{2})^{2}-(i\hbar\omega_{n}-\hbar(\vec{v}-\vec{v}_{s})\cdot\vec{k}\
)^{2})]$ (16)

is the Gaussian contribution to the grand potential, where

$E_{\vec{k}}(\psi_{0}^{2})=\sqrt[]{\left(\frac{\hbar^{2}k^{2}}{2m}-\mu_{e}+g_{0}\psi_{0}^{2}+\psi_{0}^{2}\,{\tilde{V}}(\vec{k})\right)^{2}-\big{(}\psi_{0}^{2}\,{\tilde{V}}(\vec{k})\big{)}^{2}}$
(17)

is the excitation spectrum of the Bose gas. The Hugenholtz-Pines theorem is
guaranteed imposing the saddle point condition
$\partial\Omega/\partial\psi_{0}=0$, for which the excitation spectrum becomes
gapless and the chemical potential reads

$\mu_{e}=g_{0}\psi_{0}^{2}$ (18)

Moreover, we identify the condensate density as $n_{0}=\psi_{0}^{2}$. 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:
${\tilde{V}}(\vec{k})={\tilde{V}}(-\vec{k})$. Due to this property we are able
to perform the summation over Matsubara frequencies $\omega_{n}$ in the
Gaussian grand potential (16) [27, 34], obtaining the grand potential $\Omega$
in the form

$\Omega=\Omega_{0}+\Omega_{g}^{(0)}+\Omega_{g}^{(T)}$ (19)

where $\Omega_{0}$ is given by Eq. (15),

$\Omega_{g}^{(0)}=\frac{1}{2}\sum_{\vec{k}}E_{\vec{k}}(\psi_{0}^{2})$ (20)

is the zero-temperature Gaussian contribution, written as the sum of
noninteracting elementary excitations with spectrum
$E_{\vec{k}}(\psi_{0}^{2})$, and

$\Omega_{g}^{(T)}=\frac{1}{\beta}\sum_{\vec{k}}\ln(1-e^{-\beta(E_{\vec{k}}(\psi_{0}^{2})+\hbar(\vec{v}-\vec{v}_{s})\cdot\vec{k})})$
(21)

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 $\vec{\cal P}$ of the fluid, which is obtained as follows

$\vec{\cal
P}=\frac{1}{L^{D}}\frac{\partial\Omega(\mu_{e},\psi_{0})}{\partial(-\vec{v})}\bigg{|}_{{\mu_{e}=g_{0}\psi_{0}^{2}=g_{0}n_{0}}}$
(22)

where we take the derivative of the grand potential with respect to the
velocity $-\vec{v}$ first and then we substitute the mean field value of the
chemical potential $\mu_{e}=g_{0}\psi_{0}^{2}$. We find that the momentum
density $\vec{\cal P}$ is given by

$\vec{\cal P}=\vec{\cal P}_{0}+\vec{\cal P}_{g}^{(0)}+\vec{\cal P}_{g}^{(T)}$
(23)

where, since the grand potential (19) is given by the sum of three sum
contributions, the three terms of the momentum density $\vec{\cal P}$ are
defined accordingly. In particular, we find

$\vec{\cal P}_{0}=n_{0}m\vec{v}_{s}$ (24)

The Gaussian zero-temperature contribution $\vec{\cal P}_{g}^{(0)}$ is

$\vec{\cal P}_{g}^{(0)}=f_{g}^{(0)}(n_{0})\;m\vec{v}_{s}$ (25)

where we define the zero-temperature number density contribution
$f_{g}^{(0)}(n_{0})$ as

$f_{g}^{(0)}(n_{0})=\frac{1}{2L^{D}}\sum_{\vec{k}}\frac{1}{E_{\vec{k}}(n_{0})}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}+n_{0}{\tilde{V}}(\vec{k})\bigg{)}$
(26)

The Gaussian thermal contribution $\vec{\cal P}_{g}^{(T)}$ to the momentum
density (23) is more involved:

$\vec{\cal
P}_{g}^{(T)}=\frac{1}{L^{D}}\sum_{\vec{k}}\frac{1}{e^{\beta(E_{\vec{k}}(\psi_{0}^{2})+\hbar(\vec{v}-\vec{v}_{s})\cdot\vec{k})}-1}\bigg{[}\frac{\partial
E_{\vec{k}}(\psi_{0}^{2})}{\partial\mu_{e}}\
\frac{\partial\mu_{e}}{\partial(-\vec{v})}-\hbar\vec{k}\bigg{]}\bigg{|}_{{\mu_{e}=g_{0}\psi_{0}^{2}=g_{0}n_{0}}}$
(27)

Assuming that the difference between the velocity $\vec{v}$ of the normal
fluid and the velocity $\vec{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

$\vec{\cal
P}_{g}^{(T)}=f_{g}^{(T)}(n_{0})\;m\vec{v}_{s}+n_{n}(n_{0},T)\;m(\vec{v}-\vec{v}_{s})$
(28)

where $f_{g}^{(T)}(n_{0})$ is the Gaussian thermal density contribution

$f_{g}^{(T)}(n_{0})=\frac{1}{L^{D}}\sum_{\vec{k}}\frac{1}{E_{\vec{k}}(n_{0})}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}+n_{0}{\tilde{V}}(\vec{k})\bigg{)}\frac{1}{e^{\beta
E_{\vec{k}}(n_{0})}-1}$ (29)

and we have defined the normal density of the fluid as

$n_{n}(n_{0},T)=\frac{\beta\hbar^{2}}{mDL^{D}}\sum_{\vec{k}}k^{2}\
\frac{e^{\beta E_{\vec{k}}(n_{0})}}{(e^{\beta E_{\vec{k}}(n_{0})}-1)^{2}}$
(30)

Notice that, in the thermodynamic limit $L^{D}\rightarrow\infty$, 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 $\vec{\cal P}$ as

$\vec{\cal
P}=[n_{0}+f_{g}^{(0)}(n_{0})+f_{g}^{(T)}(n_{0})]m\vec{v}_{s}+n_{n}(n_{0},T)\;m(\vec{v}-\vec{v}_{s})$
(31)

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]

$\displaystyle
n(n_{0},T)=-\frac{1}{L^{D}}\frac{\partial\Omega(\mu_{e},\psi_{0},T)}{\partial\mu_{e}}\bigg{|}_{\mu_{e}=g_{0}n_{0}}=n_{0}+f_{g}^{(0)}(n_{0})+f_{g}^{(T)}(n_{0})$
(32)

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 $\mu_{e}=g_{0}\psi_{0}^{2}$, with the identification
for the condensate density $n_{0}=\psi_{0}^{2}$. 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 $\vec{\cal P}$ as

$\vec{\cal P}=n(n_{0},T)m\vec{v}_{s}+n_{n}(n_{0},T)\;m(\vec{v}-\vec{v}_{s})$
(33)

Finally, we identify the superfluid density $n_{s}$ as

$n_{s}=n(n_{0},T)-n_{n}(n_{0},T)$ (34)

which allows us to express the momentum density $\vec{\cal P}$ as the sum of
the momentum density $n_{s}m\vec{v}_{s}$ of the superfluid part of the fluid
and the momentum density $n_{n}\left(n_{0},T\right)m\vec{v}$ of the normal
part of the fluid, namely

$\vec{\cal P}=n_{s}m\vec{v}_{s}+n_{n}\left(n_{0},T\right)m\vec{v}$ (35)

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}\approx n$ in Eqs. (26), (29) and (30)
and we consider the zero-range interaction ${\tilde{V}}(\vec{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.

##  3 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(\vec{r})$. The usual approximation to study a weakly-interacting Bose gas
of ultracold atoms is constituted by the zero-range interaction
$V(\vec{r})=g_{0}\,\delta(\vec{r})$, which in the Fourier space gives

${\tilde{V}}(\vec{k})=g_{0}$ (36)

Here we improve this approximation, considering the finite-range effective
interaction

${\tilde{V}}(\vec{k})=g_{0}+g_{2}k^{2}$ (37)

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(|\vec{r}-\vec{r}\,^{\prime}|)$, namely $g_{2}k^{2}$, where we
define

$g_{0}=\int d^{D}r\ V(r),\qquad g_{2}=-\frac{1}{2}\int d^{D}r\ r^{2}\ V(r)$
(38)

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}=\frac{4\pi\hbar^{2}a_{s}}{m},\qquad
g_{2}=\frac{2\pi\hbar^{2}a_{s}^{2}r_{e}}{m}$ (39)

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\,\sqrt{|g_{2}/g_{0}|}$
discussed in [25], obtaining

$g_{0}=\frac{4\pi\hbar^{2}}{m}\frac{1}{|\ln(na_{s}^{2})|}\qquad
g_{2}=\frac{\pi\hbar^{2}}{m}\frac{R^{2}}{|\ln(na_{s}^{2})|}$ (40)

We now explicitly implement the superfluid density $n_{s}$ calculation for the
finite-range 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

$f_{g}^{(0)}(n_{0})=\frac{S_{D}}{2(2\pi)^{D}}\int_{0}^{+\infty}dk\,k^{D-1}\frac{1}{E_{k}(n_{0})}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}\tilde{\lambda}(g_{2},n_{0})+g_{0}n_{0}\bigg{)}$
(41)

where $S_{D}={2\pi^{D/2}}/{\Gamma[D/2]}$ is the $D$-dimensional solid angle,
in which $\Gamma[x]$ is the Euler gamma function and where we define the
excitation spectrum

$E_{k}(n_{0})=\sqrt[]{\frac{\hbar^{2}k^{2}}{2m}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}\lambda(g_{2},n_{0})+2n_{0}g_{0}\bigg{)}}$
(42)

with

$\displaystyle\lambda(g_{2},n_{0})=1+\frac{4m}{\hbar^{2}}g_{2}n_{0}\qquad\tilde{\lambda}(g_{2},n_{0})=1+\frac{2m}{\hbar^{2}}g_{2}n_{0}$
(43)

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_{g}^{(0)}(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=\hbar^{2}k^{2}\lambda(g_{2},n_{0})/(4mg_{0}n_{0})$, then we extend the
spatial dimension $D$ to the complex value ${\cal D}=D-\varepsilon$. 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_{g}^{(0)}(n_{0})$ in the
form

$\displaystyle
f_{g}^{(0)}(n_{0})=\frac{\kappa^{\varepsilon}\lambda(g_{2},n_{0})^{1/2}}{4\Gamma({\cal
D}/2)}\
\bigg{(}\frac{mg_{0}n_{0}}{\pi\hbar^{2}\lambda(g_{2},n_{0})}\bigg{)}^{{\cal
D}/2}\
\bigg{[}2\frac{\tilde{\lambda}(g_{2},n_{0})}{\lambda(g_{2},n_{0})}\cdot$
$\displaystyle B\bigg{(}\frac{{\cal D}+1}{2},\frac{-{\cal
D}}{2}\bigg{)}+B\bigg{(}\frac{{\cal D}-1}{2},\frac{2-{\cal
D}}{2}\bigg{)}\bigg{]}$ (44)

where $B(x,y)$ is the Euler Beta function and $\kappa$ 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_{g}^{(T)}(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=\beta E_{k}(n_{0})$. Doing so in Eq. (29) in which the finite-
range interaction (37) is substituted, we get

$f_{g}^{(T)}(n_{0})=\frac{S_{D}}{(2\pi)^{D}k_{B}T}\int_{0}^{+\infty}dx\
\frac{dk(x)}{dx}\frac{k(x)^{D-1}}{x(e^{x}-1)}\bigg{(}\frac{\hbar^{2}k(x)^{2}}{2m}\lambda(g_{2},n_{0})+n_{0}g_{0}\bigg{)}$
(45)

with

$k(x)=\sqrt[]{\frac{2mn_{0}g_{0}}{\hbar^{2}\lambda(g_{2},n_{0})}}\
\sqrt[]{-1+\sqrt[]{1+\frac{(k_{B}T)^{2}\lambda(g_{2},n_{0})x^{2}}{n_{0}^{2}g_{0}^{2}}}}$
(46)

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_{g}^{(T)}(n_{0})$, the normal
density $n_{n}\left(n_{0},T\right)$ can be calculated analytically only in the
low-temperature regime: as before we introduce the integration variable
$x=\beta E_{k}(n_{0})$, obtaining

$n_{n}\left(n_{0},T\right)=\frac{\beta\hbar^{2}S_{D}}{mD(2\pi)^{D}}\int_{0}^{+\infty}dx\
\frac{dk}{dx}(x)\ k(x)^{D+1}\ \frac{e^{x}}{(e^{x}-1)^{2}}$ (47)

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$

###  3.1 $D=3$

As a preliminar result for the superfluid density $n_{s}$, we employ Eq. (32)
to calculate the density $n\left(n_{0},T\right)$ of bosons with finite-range
interaction. In $D=3$, the zero-temperature Gaussian contribution
$f_{g}^{(0)}(n_{0})$ is regularized simply by the means of dimensional
regularization, therefore we put $\varepsilon=0$ into the equation (3). In the
limit of small $g_{2}$ we get

$f_{g}^{(0)}(n_{0})=\frac{1}{3\pi^{2}}\
\bigg{(}\frac{mg_{0}n_{0}}{\hbar^{2}}\bigg{)}^{3/2}\bigg{[}1-\frac{12m}{\hbar^{2}}g_{2}n_{0}\bigg{]}$
(48)

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_{g}^{(0)}(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}\ \bigg{[}1+\frac{8}{3\
\sqrt[]{\pi}}(n_{0}a_{s}^{3})^{1/2}-64\
\sqrt[]{\pi}\frac{r_{e}}{a_{s}}(n_{0}a_{s}^{3})^{3/2}\bigg{]}$ (49)

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].

Figure 2: Zero-temperature condensate fraction $n_{0}/n$ as a function of the
gas parameter $na_{s}^{3}$ in $D=3$. Here we represent the finite-range
condensate fraction (black dot-dashed line) for the effective range value
$r_{e}/a_{s}=-10$ in comparison with the result for zero-range interaction
(blue solid line) and the Bogoliubov’s result (red dashed line). Notice that,
with a weak dependence on the choice of $r_{e}$, the finite-range correction
becomes relevant for values of the gas parameter $na_{s}^{3}$ greater than
$10^{-3}$.

Let us also calculate the finite-temperature contribution $f_{g}^{(T)}(n_{0})$
in $D=3$, given by the integration of Eq. (45) in the low-temperature regime

$f_{g}^{(T)}(n_{0})=\frac{(k_{B}T)^{2}}{12(n_{0}g_{0})^{1/2}}\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{3/2}\bigg{(}1-\frac{\pi^{2}(k_{B}T)^{2}\lambda(g_{2},n_{0})}{20(n_{0}g_{0})^{2}}\bigg{)}+o(k_{B}T)^{5}$
(50)

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\left(n_{0},T\right)=n_{0}+\frac{1}{3\pi^{2}}\
\bigg{(}\frac{mg_{0}n_{0}}{\hbar^{2}}\bigg{)}^{3/2}\bigg{[}1-\frac{12m}{\hbar^{2}}g_{2}n_{0}\bigg{]}+\frac{(k_{B}T)^{2}}{12(n_{0}g_{0})^{1/2}}\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{3/2}$
(51)

We also rewrite the general form of $n(n_{0},T)$ in terms of the three-
dimensional gas parameter $na_{s}^{3}$ and the effective range $r_{e}$,
employing the explicit form of the couplings $g_{0}$ and $g_{2}$ given by Eq.
(39), namely

$n\left(n_{0},T\right)=n_{0}+\frac{8(n_{0}a_{s})^{3/2}}{3\sqrt{\pi}}\frac{1-8\pi^{2}n_{0}a_{s}^{2}r_{e}}{(1+8\pi^{2}n_{0}a_{s}^{2}r_{e})^{2}}+\frac{(k_{B}T)^{2}}{24\sqrt{\pi}(n_{0}a_{s})^{1/2}}\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{2}$
(52)

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}\rightarrow 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

$\displaystyle\frac{n_{0}}{n}=$ $\displaystyle 1-\frac{8}{3\
\sqrt[]{\pi}}(na_{s}^{3})^{1/2}\bigg{[}1-24\pi(na_{s}^{3})\,\frac{r_{e}}{a_{s}}+\frac{\pi^{2}}{16\zeta(3/2)^{4/3}(na_{s}^{3})^{2/3}}\cdot$
(53)
$\displaystyle\bigg{(}\frac{T}{T_{BEC}}\bigg{)}^{2}\bigg{(}1-\frac{\pi^{2}}{80\zeta(3/2)^{4/3}(na_{s}^{3})^{2/3}}\bigg{(}\frac{T}{T_{BEC}}\bigg{)}^{2}\big{(}1+8\pi\frac{r_{e}}{a_{s}}(na_{s}^{3})\big{)}\bigg{)}\bigg{]}$

where we have rescaled the temperature in terms of
$T_{BEC}=(2\pi\hbar^{2}n^{2/3})/(mk_{B}\zeta(3/2)^{2/3})$ for noninteracting
bosons, and $\zeta(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

$n_{n}=\frac{2\pi^{2}}{45}\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{3/2}\frac{(k_{B}T)^{4}}{(n_{0}g_{0})^{5/2}}\bigg{(}1-\frac{5\pi^{2}\lambda(g_{2},n_{0})(k_{B}T)^{2}}{2(n_{0}g_{0})^{2}}\bigg{)}+o(k_{B}T)^{7}$
(54)

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

$\displaystyle n_{s}=$ $\displaystyle n_{0}+\frac{1}{3\pi^{2}}\
\bigg{(}\frac{mg_{0}n_{0}}{\hbar^{2}}\bigg{)}^{3/2}+\frac{(k_{B}T)^{2}}{12(n_{0}g_{0})^{1/2}}\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{3/2}-$
(55)
$\displaystyle\frac{(8\lambda(n_{0})+3)\pi^{2}}{360}\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{3/2}\frac{(k_{B}T)^{4}}{(n_{0}g_{0})^{5/2}}$

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

$\displaystyle\frac{n_{s}}{n}=$ $\displaystyle
1-\frac{\pi^{7/4}}{45(na_{s}^{3})^{5/6}}\bigg{(}\frac{T}{T_{BEC}}\bigg{)}^{4}\bigg{[}1-\frac{5\pi^{2}}{8\,\zeta(3/2)^{4/3}(na_{s}^{3})^{2/3}}\bigg{(}\frac{T}{T_{BEC}}\bigg{)}^{2}\cdot$
(56)
$\displaystyle\bigg{(}1+8\pi\frac{r_{e}}{a_{s}}(na_{s}^{3})\bigg{)}\bigg{]}$

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}\ll 1$ regime but,
more generally, one must consider the full solution of Eq. (47).

Figure 3: Superfluid fraction $n_{s}/n$ in $D=3$ for the gas parameter value
$na_{s}^{3}=10^{-4}$, obtained as the numerical integral of Eq. (47). The
zero-range interaction result is reported as the red dashed line, while the
finite-range corrections are reported as the solid blue line (for
$r_{e}/a_{s}=10^{2}$) and the black dot-dashed line (for
$r_{e}/a_{s}=10^{-2}$).

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_{s}^{3}$ 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.

###  3.2 $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 zero-temperature density
contribution $f_{g}^{(0)}(n_{0})$ from Eq. (3) where $D=2$ is substituted.
Notice that $f_{g}^{(0)}(n_{0})$ is obtained as a sum of the two terms inside
the square bracket of (3): while the first term is finite, the divergence of
the second term is healed by performing a Taylor expansion around $\epsilon=0$
and deleting the $o(\epsilon^{-1})$ divergence [41, 42]. The regularized zero-
temperature density contribution reads

$f_{g}^{(0)}(n_{0})=\frac{mg_{0}n_{0}}{4\pi\hbar^{2}}\frac{1}{\lambda(g_{2},n_{0})^{3/2}}\bigg{[}\tilde{\lambda}(g_{2},n_{0})+\frac{2m}{\hbar^{2}}g_{2}n_{0}\ln\bigg{(}\frac{2\epsilon_{0}}{g_{0}n_{0}}\lambda(g_{2},n_{0})\bigg{)}\bigg{]}$
(57)

where we identify the ultraviolet energy scale $\epsilon_{0}$ as

$\epsilon_{0}=\frac{4\pi\hbar^{2}\kappa^{2}}{me^{\gamma}}$ (58)

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
$\frac{2m}{\hbar^{2}}g_{2}n_{0}$, thus

$f_{g}^{(0)}(n_{0})=\frac{1}{4\pi}\
\frac{mg_{0}n_{0}}{\hbar^{2}}+\frac{1}{2\pi}\
\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{2}g_{0}g_{2}n_{0}^{2}\
\bigg{[}\ln\bigg{(}\frac{2\epsilon_{0}}{g_{0}n_{0}}\bigg{)}-2\bigg{]}$ (59)

In this limit, the zero-temperature density $n\left(n_{0},T=0\right)$ reads

$n\left(n_{0},T=0\right)=n_{0}+\frac{1}{4\pi}\
\frac{mg_{0}n_{0}}{\hbar^{2}}+\frac{1}{2\pi}\
\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{2}g_{0}g_{2}n_{0}^{2}\
\bigg{[}\ln\bigg{(}\frac{2\epsilon_{0}}{g_{0}n_{0}}\bigg{)}-2\bigg{]}$ (60)

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}\rightarrow 0$. As before, in this regime we can approximate
$n_{0}^{2}\approx n_{0}n$ in the third term of Eq. (60) and $n_{0}\approx n$
inside the logarithm, since it is a subleading term with respect to the first
one. Let us come back to $\varepsilon_{0}$, by chosing the cutoff as
$\kappa=2\pi/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_{s}^{2}$ and the characteristic range $R$,
namely

$\frac{n_{0}}{n}=1-\frac{1}{|\ln(na_{s}^{2})|}-\frac{2\pi(R/a_{s})^{2}}{|\ln(na_{s}^{2})|^{2}}\
\ln\bigg{(}\frac{8\pi^{2}|\ln(na_{s}^{2})|}{na_{s}^{2}e^{\gamma+2}}\bigg{)}$
(61)

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_{s}^{3}\ll 1$.

Figure 4: Condensate fraction $n_{0}/n$ for bosons in $D=2$ at $T=0$, reported
as a function of the gas parameter $na_{s}^{2}$. The blue solid line is the
condensate fraction for bosons with zero-range interaction (Z-R), obtained
from the numerical solution of Eq. (60) with $R=0$. The black dot-dashed line
is the condensate fraction for bosons with finite-range interaction (F-R)
given by Eq. (60) with a characteristic range $R$ value given by $R=2a_{s}$
and a ultraviolet cutoff $\kappa={2\pi}/{a_{s}}$. The red dashed line is the
analytical formula of Schick [43], i.e. Eq. (61) with $R=0$. The inset
highlights the small differences between our theoretical scheme and the Schick
one in the weakly-interacting regime.

Finally, following the three-dimensional case, we may want to calculate also
the finite-temperature density contribution $f_{g}^{(T)}(n_{0})$. However,
substituting $D=2$ in Eq. (45) we find that

$f_{g}^{(T)}(n_{0})=\frac{mk_{B}T}{2\pi\hbar^{2}}\int_{0}^{+\infty}dx\
\frac{1}{e^{x}-1}$ (62)

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}\left(n_{0},T=0\right)$ 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

$\displaystyle n_{n}\left(n_{0},T=0\right)=$
$\displaystyle\frac{3m}{2\pi\hbar^{2}}\frac{(k_{B}T)^{3}}{(n_{0}g_{0})^{2}}\bigg{(}\zeta(3)-$
(63)
$\displaystyle\frac{15\zeta(5)\lambda(g_{2},n_{0})}{(n_{0}g_{0})^{2}}(k_{B}T)^{2}\bigg{)}+o(k_{B}T)^{6}$

Since the thermal contribution to the density $f_{g}^{(T)}(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}\left(n,T=0\right)$,
namely as a function of the density $n$. In this way we calculate the
superfluid density substituting Eq. (63) into Eq. (34)

$\displaystyle n_{s}=$ $\displaystyle
n-\frac{3m}{2\pi\hbar^{2}}\frac{(k_{B}T)^{3}}{(g_{0}\,n_{0}\left(n,T=0\right))^{2}}\bigg{(}\zeta(3)-$
(64)
$\displaystyle\frac{15\zeta(5)\lambda(g_{2},n_{0}\left(n,T=0\right))}{(g_{0}\,n_{0}\left(n,T=0\right))^{2}}(k_{B}T)^{2}\bigg{)}$

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}\rightarrow
0$, in which we can approximate the condensate density inside Eq. (64) as
$n_{0}\left(n,T=0\right)\approx n$. In the context of this approximation, we
also substitute in Eq. (64) the explicit form of the parameter
$\lambda(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

$\displaystyle\frac{n_{s}}{n}=1-\frac{3}{32\pi^{3}}|\ln(na_{s}^{2})|^{2}\bigg{(}\frac{T}{T^{*}}\bigg{)}^{3}\bigg{[}\zeta(3)\,-\frac{15\zeta(5)}{16\pi^{2}}\bigg{(}\frac{T}{T^{*}}\bigg{)}^{2}|\ln(na_{s}^{2})|^{2}\cdot$
$\displaystyle\bigg{(}1+4\pi\frac{R^{2}}{a_{s}^{2}}\frac{na_{s}^{2}}{|\ln(na_{s}^{2})|}\bigg{)}\bigg{]}$
(65)

where we rescale the result in terms of the temperature
$T^{*}=\hbar^{2}n/(mk_{B})$ of quantum degeneracy.

##  4 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_{s}^{D}\ll 1$ and $nr_{e}^{D}\ll 1$. Moreover,
finite-temperature results must be considered in the limit $k_{B}T/(g_{0}n)\ll
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}\leq 1$ and in $D=2$ for $a_{s}/R\leq 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(\vec{r})$.

An extension of this work consists in the numerical calculation of the thermal
integrals $f_{g}^{(T)}(n_{0})$ and $n_{n}\left(n_{0},T\right)$ 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.

## Acknowledgements

The authors thank Francesco Ancilotto and Flavio Toigo for useful discussions.
LS acknowledges for partial support the FFABR grant of Italian Ministry of
Education, University and Research.

## Appendix

In this appendix, we calculate the zero-temperature Gaussian grand potential
per unit of volume $\Omega^{(0)}/L^{D}$, which is given by Eqs. (15) and (20)
in which the mean field condition Eq. (18) is substituted, namely

$\frac{\Omega^{(0)}}{L^{D}}=-\frac{\mu_{e}^{2}}{2g}+\frac{S_{D}}{2(2\pi)^{D}}\int_{0}^{+\infty}dk\,k^{D-1}\
\sqrt[]{\frac{\hbar^{2}k^{2}}{2m}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}\lambda(g_{2},\mu_{e}/g_{0})+2\mu_{e}\bigg{)}}$
(66)

where $\lambda$ is defined in Eq. (43). Performing dimensional regularization
of the ultraviolet divergence accordingly to the superfluid density
calculation, we obtain the regularized zero-temperature Gaussian grand
potential $\Omega_{g}^{(0)}/L^{D}$

$\displaystyle\frac{\Omega^{(0)}}{L^{D}}=-\frac{\mu_{e}^{2}}{2g}-\frac{\mu_{e}^{\frac{D+2}{2}}}{2\
\pi^{\frac{D+1}{2}}}\
\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{D/2}\frac{1}{\lambda(g_{2},\mu_{e}/g_{0})^{\frac{D+1}{2}}}\bigg{[}1+$
$\displaystyle\frac{\varepsilon}{2}\ln\bigg{(}\frac{\pi\hbar^{2}\kappa^{2}\lambda(g_{2},\mu_{e}/g_{0})}{m\mu_{e}}\bigg{)}+o(\varepsilon^{2})\bigg{]}\frac{\Gamma[(D-\varepsilon+1)/2]\
\Gamma[(\varepsilon-D-2)/2]}{\Gamma[(D-\varepsilon)/2]}$ (67)

In the three-dimensional case, it is sufficient to substitute the dimension
$D=3$ and put $\varepsilon=0$ in Eq. (Appendix), to get the regularized zero-
temperature Gaussian grand potential as

$\frac{\Omega^{(0)}}{L^{3}}=-\frac{\mu_{e}^{2}}{2g}-\frac{\mu_{e}^{2}}{2g}\frac{8}{15\pi^{2}}\
\bigg{(}\frac{m}{\hbar^{2}}\bigg{)}^{3/2}\
\frac{\mu_{e}^{5/2}}{{\lambda(g_{2},\mu_{e}/g_{0})}^{2}}$ (68)

which reproduces a previously known result [23]. In $D=2$ we expand Eq.
(Appendix) and retain only $o(\varepsilon^{0})$ terms, then we identify the
energy cutoff $\epsilon_{0}$ as

$\epsilon_{0}=\frac{4\pi\hbar^{2}\kappa^{2}}{m\exp(\gamma+1-\frac{4\pi\hbar^{2}{\lambda(g_{2},\mu_{e}/g_{0})}^{3/2}}{mg_{0}})}$
(69)

to get the zero-temperature Gaussian grand potential

$\frac{\Omega^{(0)}}{L^{2}}=-\frac{m\mu_{e}^{2}}{8\pi\hbar^{2}{\lambda(g_{2},\mu_{e}/g_{0})}^{3/2}}\
\bigg{[}\ln\bigg{(}\frac{\epsilon_{0}}{\mu_{e}}\lambda(g_{2},\mu_{e}/g_{0})\bigg{)}+\frac{1}{2}\bigg{]}$
(70)

This equation of state corrects the one obtained in Ref. [25]: here
$\lambda(g_{2},\mu_{e}/g_{0})$ appears also inside the logarithm. Moreover, it
is important to stress that, in the case of zero-range 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.

## References

## References

  * [1] Tilley D R and Tilley J 1990 Superfluidity and Superconductivity (CRC Press)
  * [2] Pezzé L, Smerzi A, Oberthaler M K, Schmied M and Treutlein P 2018 Rev. Mod. Phys. 90 035005
  * [3] Streltsov A, Adesso G and Plenio M 2017 Rev. Mod. Phys. 89 041003
  * [4] Chang D E, Douglas J S, Gonzalez-Tudela A, Hung C L and Kimble H J 2018 Rev. Mod. Phys. 90 031002
  * [5] Bardeen J, Cooper L N and Schrieffer J R 1957 Phys. Rev. 106 1
  * [6] Larkin A and Varlamov A 2005 Theory of fluctuations in superconductors (Oxford: Oxford University Press)
  * [7] Hohenberg P C and Martin P C 1965 Annals of Physics 281 636-705
  * [8] Anderson P W 1966 Rev. Mod. Phys. 38 2
  * [9] Ceperley D M 1995 Rev. Mod. Phys. 67 279
  * [10] Schmitt A 2014 Introduction to Superfluidity: Field-theoretical Approach and Applications (Switzerland: Springer)
  * [11] Landau L D 1941 Phys. Rev. 60 356
  * [12] Landau L D and Lifshitz E M 1981 Statistical Physics - Part 2 (Oxford: Pergamon Press)
  * [13] Leggett A J 2006 Quantum Liquids: Bose Condensation and Cooper Pairing in Condensed-matter Systems (Oxford: Oxford University Press)
  * [14] London F 1938 Nature 141 643
  * [15] London F 1938 Phys. Rev. 54 937
  * [16] Bogoliubov N N 1947 J. Phys. U.S.S.R. 11 23
  * [17] Huang K 1987 Statistical Mechanics (Wiley & sons)
  * [18] Stringari S and Pitaevskii L P 2016 Bose-Einstein Condensation and Superfluidity (Oxford University Press)
  * [19] Hadzibabic Z and Dalibard J 2011 Rivista del Nuovo Cimento 34 389
  * [20] Dalibard J, Gerbier F, Juzeliunas G and Öhberg P 2011 Rev. Mod Phys. 83 1523
  * [21] Braaten E, Hammer H W and Hermans S 2001 Phys. Rev. A 63 063609
  * [22] Andersen J O 2004 Rev. Mod. Phys. 76 599
  * [23] Cappellaro A and Salasnich L 2017 Phys. Rev. A  95 033627
  * [24] Yukalov V I and Yukalova E P 2016 Laser Phys. 26 045501
  * [25] Salasnich L 2017 Phys. Rev. Lett.  118 130402
  * [26] Beane S R 2018 Eur. Phys. J. D 72 55
  * [27] Schakel A M J 2008 Boulevard of Broken Symmetries (Singapore: World Scientific)
  * [28] Fisher M E, Barber M and Jasnow D 1973 Phys. Rev. A  8 1111
  * [29] Taylor E, Griffin A, Fukushima N and Ohashi Y, (2006), Phys. Rev. A 74 063626
  * [30] Nagaosa N 1999 Quantum Field Theory in Condensed Matter Physics (Springer: Verlag)
  * [31] Boudjemâa A and Benarous M 2011 Phys. Rev. A 84 043633
  * [32] Boudjemâa A 2012 Phys. Rev. A 86 043608
  * [33] Boudjemaa A and Guebli N 2017 J. Phys. A: Math Theor. 50 425004
  * [34] Altland A and Simons B D 2010 Condensed Matter Field Theory  (Cambridge: Cambridge University Press)
  * [35] Salasnich L and Toigo F 2016 Phys. Rep. 640 1-29
  * [36] Josephson B D 1966 Phys. Lett.  21 608
  * [37] Svistunov B, Babaev E, Prokof’ev N 2015 Superfluid States of Matter (CRC Press: Boca Raton)
  * [38] Astrakharchik G E, Boronat J, Casulleras J, Kurbakov I L, and Lozovik Yu E 2009 Phys. Rev. A 79 051602
  * [39] ’t Hooft G and Veltman M 1972 Nucl. Phys. B 44 189-213
  * [40] Leibbrandt G 1975 Rev. Mod. Phys.  47 849
  * [41] Kleinert H, Schulte-Frohlinde V 2001 Critical Properties of $\phi^{4}$ Theories (Singapore: World Scientific)
  * [42] Zeidler E 2009 Quantum Field Theory II: Quantum Electrodynamics (Berlin: Springer)
  * [43] Schick M 1971 Phys. Rev. A  3 1067
  * [44] Mermin N D and Wagner H 1966 Phys. Rev. Lett. 17 1133
  * [45] Stoof H T C, Dickerscheid D B M and Gubbels K 2009 Ultracold Quantum Fields (Dordrecht: Springer)
  * [46] Kosterlitz L M and Thouless D J 1973 J. Phys. C 6 1181
  * [47] Cappellaro A and Salasnich L 2017 Phys. Rev. A 96 063610
  * [48] Mora C and Castin Y 2009 Phys. Rev. Lett. 102 180404