Zero-temperature equation of state of a two-dimensional bosonic quantum fluid with finite-range interaction

TL;DR

This paper derives the zero-temperature equation of state for a two-dimensional bosonic quantum fluid with finite-range interactions, using a hydrodynamic approach and regularization techniques to handle divergences.

Contribution

It introduces a method to calculate the equation of state for 2D bosonic systems with finite-range interactions, including regularization of divergences.

Findings

Derived the 2D equation of state for finite-range interactions.

Calculated superfluid equations of motion at zero temperature.

Regularized ultraviolet divergences using an improved dimensional regularization.

Abstract

We derive the two-dimensional equation of state for a bosonic system of ultracold atoms interacting with a finite-range effective interaction. Within a functional integration approach, we employ an hydrodynamic parametrization of the bosonic field to calculate the superfluid equations of motion and the zero-temperature pressure. The ultraviolet divergences, naturally arising from the finite-range interaction, are regularized with an improved dimensional regularization technique.

Figures

Equations

$${\cal L}=\bar{\psi}(\vec{r},\tau)\bigg{(}\hbar\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu\bigg{)}\psi(\vec{r},\tau)\\ +\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},$$

$$S[\bar{\psi },\psi ]={\int }_0^{\beta \hslash }d\tau {\int }_{L^D}{d}^{D}r\mathcal{L}$$

$$\psi (\vec{r},\tau )=\sqrt{\rho (\vec{r},\tau )}{e}^{i\theta (\vec{r},\tau )},$$

$${\vec{v}}_s=\frac{\hslash }{m}\vec{\nabla }\theta .$$

$$\vec{\nabla }\times {\vec{v}}_s=0.$$

$$\mathcal{L}=-\mu \rho +i\hslash \rho {\partial }_{\tau }\theta +\frac{{\hslash }^2}{8m\rho }(\nabla \rho {)}^2+\frac{{\hslash }^2\rho }{2m}(\nabla \theta {)}^2+\frac{1}{2}\int d^{D}r{}^{\prime }\rho (\vec{r}{}^{\prime },\tau )V(|\vec{r}-\vec{r}{}^{\prime }|)\rho (\vec{r},\tau ),$$

$$\frac{\partial \rho }{\partial t}+\vec{\nabla }\cdot (\rho {\vec{v}}_s)=0$$

$$m\frac{\partial\vec{v}_{s}}{\partial t}+\vec{\nabla}\bigg{(}\frac{1}{2}mv_{s}^{2}-\mu+\int d^{D}r\,^{\prime}\ \rho(\vec{r}\,^{\prime},\tau)\ V(|\vec{r}-\vec{r}\,^{\prime}|)-\frac{\hbar^{2}}{2m\rho^{1/2}}\nabla^{2}(\rho^{1/2})\bigg{)}=0,$$

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

$$\mathcal{Z}=\int \mathcal{D}[\rho ,\theta ]e^{-\frac{S[\rho ,\theta ]}{\hslash }}.$$

$$\rho (\vec{r},\tau )={\rho }_0+\delta \rho (\vec{r},\tau ),$$

$$\begin{array}{rl}\mathcal{L}=& -\mu {\rho }_0-\mu \delta \rho +i\hslash {\rho }_0{\partial }_{\tau }\theta +\frac{{\hslash }^2}{8m{\rho }_0}(\nabla \delta \rho {)}^2+\frac{{\hslash }^2{\rho }_0}{2m}(\nabla \theta {)}^2+\\ & \frac{1}{2}\int d^D{r}^{\prime }V(|\vec{r}-\vec{r}{}^{\prime }|)({\rho }_0^2+\delta \rho (\vec{r},\tau )+\delta \rho (\vec{r}{}^{\prime },\tau )+\rho (\vec{r},\tau )\delta \rho (\vec{r}{}^{\prime },\tau )),\end{array}$$

$$\delta \rho (\vec{r},\tau )$$

$$\theta (\vec{r},\tau )$$

$$\delta \rho (\vec{k},{\omega }_n)$$

$$\theta (\vec{k},{\omega }_n)$$

$$S=S_0+S_g.$$

$$S_{0}=\beta\hbar L^{D}\big{(}-\mu\rho_{0}+\frac{1}{2}g_{0}\rho_{0}^{2}\big{)},$$

$$\Omega_{0}=\big{(}-\mu\rho_{0}+\frac{1}{2}g_{0}\rho_{0}^{2}\big{)}L^{D}.$$

$$S_{g}=\beta\hbar\sum_{\vec{k}\,\omega_{n}}\bigg{[}\frac{\hbar^{2}k^{2}\rho_{0}}{2m}\theta(k)\theta(-k)+\bigg{(}\frac{\hbar^{2}k^{2}}{8m\rho_{0}}+\frac{\tilde{V}(k)}{2}\bigg{)}\delta\rho(k)\delta\rho(-k)+\hbar\omega_{n}\theta(k)\delta\rho(-k)\bigg{]},$$

$$S_g=\frac{\hslash }{2}\sum _{\vec{k}{\omega }_n}\left(\begin{array}{cccc}\theta (k)& \theta (-k)& \delta \rho (k)& \delta \rho (-k)\end{array}\right)\text{M}(k)\left(\begin{array}{c}\theta (k)\\ \theta (-k)\\ \delta \rho (k)\\ \delta \rho (-k)\end{array}\right),$$

$$\textbf{M}(k)=\beta\left(\begin{array}[]{c c c c}0&\frac{\hbar^{2}k^{2}\rho_{0}}{m}&0&\hbar\omega_{n}\\ \frac{\hbar^{2}k^{2}\rho_{0}}{m}&0&-\hbar\omega_{n}&0\\ 0&-\hbar\omega_{n}&0&\frac{\hbar^{2}k^{2}}{4m\rho_{0}}+\tilde{V}(k)\\ \hbar\omega_{n}&0&\frac{\hbar^{2}k^{2}}{4m\rho_{0}}+\tilde{V}(k)&0\\ \end{array}\right).$$

$${\mathcal{Z}}_g=\prod _{\begin{array}{c}\vec{k}{\omega }_n\\ k_z>0\end{array}}[\det \text{M}(k){]}^{-1/2},$$

$${\Omega }_g=\frac{1}{2\beta }\sum _{\vec{k}{\omega }_n}\ln [{\beta }^2({\hslash }^2{\omega }_n^2+E_k^2)].$$

$$E_{k}=\sqrt[]{\frac{\hbar^{2}k^{2}}{2m}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}+2\rho_{0}\tilde{V}(k)\bigg{)}},$$

$${\rho }_0=\frac{\mu }{g_0}$$

$$E_{k}^{B}=\sqrt[]{\frac{\hbar^{2}k^{2}}{2m}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}+2\mu\frac{\tilde{V}(k)}{g_{0}}\bigg{)}}.$$

$$\Omega ={\Omega }_0+{\Omega }_g^{(0)}+{\Omega }_g^{(T)},$$

$${\Omega }_g^{(0)}=\frac{1}{2}\sum _{\vec{k}}{E}_k^B$$

$${\Omega }_g^{(T)}=\frac{1}{\beta }\sum _{\vec{k}}\ln (1-e^{-\beta E_{\vec{k}}^B})$$

$$P(\mu ,T=0)=\frac{{\mu }^2}{2g_0}-\frac{1}{2L^D}\sum _{\vec{k}}{E}_k^B.$$

$$P(\mu,T=0)=\frac{\mu^{2}}{2g_{0}}-\frac{1}{2}\int\frac{d^{D}k}{(2\pi)^{D}}\;\sqrt[]{\frac{\hbar^{2}k^{2}}{2m}\bigg{(}\frac{\hbar^{2}k^{2}}{2m}+2\mu\frac{\tilde{V}(k)}{g_{0}}\bigg{)}},$$

Full text

\articlenumber

x \doinum10.3390/—— \pubvolumexx

\TitleZero-temperature equation of state of a two-dimensional bosonic quantum fluid with finite-range interaction \AuthorAndrea Tononi 1 \AuthorNamesAndrea Tononi

\corresCorrespondence: andrea.tononi@phd.unipd.it

1 Introduction

A fluid perspective on the study of bosonic gases made of ultracold atoms may find an origin in the pioneering work of Lev Landau landau1941 , who correctly described the superfluid behavior of He-4 observed in the experiments kapitza1938 . From those years, several technical advances allowed a precise experimental control of the atomic gases, culminating in the achievement of Bose-Einstein condensation in 1995 anderson1995 ; davis1995 ; hulet1995 . Despite the rich phenomenology deriving from the variety of different trapping potentials dalfovo1999 , the theoretical and experimental study of uniform gases gives a fundamental insight into the intrinsic properties of the condensates. Particularly interesting is the case of two spatial dimensions, in which the quantum and thermal fluctuations play a fundamental role merminwagner ; hohenberg1967 , justifying the necessity of a beyond mean field theory. Historical results for a two-dimensional Bose gas were obtained by Schick, who calculated the thermodynamics of a gas of hard-spheres schick1971 , improved by Popov derivation of the equation of state for a weakly-interacting superfluid popov1972 . Recent works provide an extension of Popov approach pastukhov2018 , while the nonuniversal corrections to the equation of state in $D=2$, arising from a finite-range interaction between the atoms, have been studied salasnich2017 ; beane . Thank to the tunability of interparticle interactions, it is possible investigate the static and dynamical properties of homogeneous quantum fluids in $D$-spatial dimensions in regimes where finite-range corrections are relevant braaten ; salasnichcappellaro ; cappe . In this work, we provide an alternative derivation of the zero-temperature equation of state by adopting an explicit superfluid parametrization of the bosonic field. In particular, we develop an improved dimensional regularization technique to regularize the zero-temperature pressure of a bosonic quantum fluid, whose particles interact with a finite-range interaction.

2 Zero-temperature equation of state of a two-dimensional

bosonic quantum fluid

2.1 Superfluid parametrization of the bosonic field.

We introduce the Euclidean Lagrangian ${\cal L}$ of a uniform quantum fluid of bosonic particles with mass $m$, described by the complex field $\psi(\vec{r},t)$, namely nagaosa

[TABLE]

where $\hbar$ is Planck constant, $\mu$ is the chemical potential and we suppose that the particles interact with the isotropic two-body potential $V(|\vec{r}-\vec{r}\,^{\prime}|)$. The imaginary time $\tau$ is introduced for uniformity with a functional integration approach, but real time $t$ can be recovered in any moment performing the Wick rotation $\tau\rightarrow it$. According to the least action principle, the Euler-Lagrange equations of the system are obtained as the functional derivative of the action $S[\bar{\psi},\psi]$, which reads

[TABLE]

where $L^{D}$ is the volume in $D$ dimensions containing the particles and $\beta=1/(k_{B}T)$, with $T$ the absolute temperature and $k_{B}$ the Boltzmann constant. Until the end of this paper, for reasons connected to the dimensional regularization of the final results, we will not fix explicitly the spatial dimension to $D=2$. The minimization of the action of Eq. (2) gives the Gross-Pitaevski equation for the complex field $\psi(\vec{r},\tau)$, which constitutes the macroscopic wavefunction of the condensate leggett . In this work, however, we adopt a superfluid perspective through the following phase-amplitude parametrization of the bosonic field salasnichlaser

[TABLE]

where $\rho(\vec{r},\tau)=|\psi(\vec{r},\tau)|^{2}$ is a real field describing the system local density and $\theta(\vec{r},\tau)$ is the phase field, which must be included due to the complex nature of the order parameter $\psi(\vec{r},\tau)$. This field transformation allows us to introduce the superfluid velocity $\vec{v}_{s}$, which is proportional to the gradient of the phase, namely

[TABLE]

We emphasize that the phase field $\theta(\vec{r},\tau)$ is defined in the compact interval $[0,2\pi]$ and is therefore periodic of $2\pi$: this fact constitutes the origin of many topological phenomena in condensed-matter physics. Indeed, here we focus on two-dimensional systems, where the singularities of the phase field

• the vortices - are responsible for the Berezinski-Kosterlitz-Thouless (BKT) transition berezinski1971 ; kosterlitz1973 . However, in the following we will study only the zero-temperature properties, for which the vortex-antivortex phenomenology does not play a fundamental role and can be neglected. We will then assume that the domain of definition of the phase field $\theta(\vec{r},\tau)$ can be extended to $\mathbb{R}$ and that its spatial and time derivatives are well defined everywhere. In this case, the superfluid flow is irrotational, thus it has zero vorticity

[TABLE]

We now substitute the parametrization of Eq. (3) in the Lagrangian (1), obtaining

[TABLE]

where we omit for simplicity the dependence of the fields on their coordinates $(\vec{r},\tau)$. The minimization of the action (2), which now becomes a functional of $\rho$ and $\theta$: $S=S[\rho,\theta]$, leads to the Euler-Lagrange equations for these fields. Recovering the real time $t\rightarrow-i\tau$, we get the hydrodynamic equations

[TABLE]

and

[TABLE]

which are the continuity equation and the equation of motion of a superfluid with velocity $\vec{v_{s}}$ and number density $\rho$. Notice that, inside the parenthesis, Eq. (8) contains the quantum pressure term: if the density $\rho$ of the fluid is slowly varying and this contribution can be neglected, Eq. (8) reduces to the familiar Euler equation of motion for an irrotational fluid without viscosity. For a self-consistent derivation of these equations from the Gross-Pitaevski equation we refer the reader to Ref. coste .

2.2 Zero-temperature equation of state.

We now adopt the superfluid parametrization of Eq. (3) to derive the zero-temperature equation of state of the quantum fluid, namely a relation at $T=0$ between the pressure $P$ and the chemical potential $\mu$. For the finite-range interaction, an explicit implementation of this relation will be given in the next section.

In the grand canonical ensemble we calculate the pressure of the bosonic fluid as $P=-\Omega/L^{D}$, where $\Omega$ is the grand potential

[TABLE]

and ${\cal Z}$ is the grand canonical partition function, which, within a functional integration perspective, can be calculated as

[TABLE]

To perform the explicit functional integration of the Lagrangian (6), we rewrite the local density field $\rho(\vec{r},\tau)$ as

[TABLE]

where $\rho_{0}$ is the condensate density of the system in the broken-symmetry phase and $\delta\rho(\vec{r},\tau)$ is a real field describing the density fluctuations.

We substitute the field transformation (11) in the Lagrangian of Eq. (6), obtaining

[TABLE]

where we keep only terms up to second order in the fluctuation fields $\delta\rho(\vec{r},\tau)$ and $\theta(\vec{r},\tau)$, thus making a Gaussian (one-loop) approximation.

Considering the Lagrangian of Eq. (12) inside the action $S$, which now becomes the functional $S=S[\delta\rho,\theta]$, it is particularly convenient to express it in terms of the Fourier series of the fluctuation fields, namely

[TABLE]

where $\omega_{n}=2\pi n/(\beta\hbar)$ are the bosonic Matsubara frequencies. Notice that, since we are supposing that the phase field $\theta(\vec{r},\tau)$ is defined on $\mathbb{R}$, its Fourier components are non-numerable and can assume continuous values, thus they can be treated like ordinary functional integral variables. The action in the Fourier space is obtained by simply substituting these expressions in $S$ and using the definition of the $D+1$-dimensional delta function. Moreover, we also substitute the Fourier series $\tilde{V}(k)$ of the real space interaction potential and we define with $g_{0}$ the zero-range interaction strength $g_{0}=\tilde{V}(k=0)$. In this way, the action can be rewritten as the sum of two contributions

[TABLE]

The first is the action of the homogeneous system $S_{0}$, namely

[TABLE]

which does not depend on the functional integration variables: using Eqs. (9) and (10) one can employ $S_{0}$ to calculate $\Omega_{0}$, the mean field contribution to the grand potential

[TABLE]

The second contribution to the action $S$ is the Gaussian action $S_{g}$, which is given by

[TABLE]

where, for simplicity of notation, we define $\delta\rho(\pm k)=\delta\rho(\pm\vec{k},\pm\omega_{n})$ and $\theta(\pm k)=\theta(\pm\vec{k},\pm\omega_{n})$. Since $S_{g}$ is quadratic in the fluctuation fields $\delta\rho(k)$ and $\theta(k)$, one can rewrite it in the following matricial form

[TABLE]

where $\textbf{M}(k)$, the inverse of the propagator, is the $4\times 4$ matrix

[TABLE]

The functional integral of the real fluctuation fields $\theta(k)$ and $\delta\rho(k)$ can be performed explicitly altland2010 , obtaining the corresponding Gaussian grand canonical partition function $\mathcal{Z}_{g}$ as

[TABLE]

which, considering the definition of the grand potential of Eq. (9), leads to the Gaussian contribution to the grand potential

[TABLE]

Here, we find the gapless excitation spectrum $E_{k}$ of the quantum fluid in the form

[TABLE]

where, within a perturbative approach, $\rho_{0}$ is determined by the saddle point condition $\partial\Omega_{0}/\partial\rho_{0}=0$, which leads to

[TABLE]

and whose substitution in the excitation spectrum gives $E_{k}^{B}$, the renowned Bogoliubov spectrum bogoliubov1947

[TABLE]

The sum over the Matsubara frequencies $\omega_{n}$ in the Gaussian grand potential of Eq. (21) can be performed according to the prescriptions described in the Appendix, obtaining the grand potential as the sum of three contributions

[TABLE]

where $\Omega_{0}=-L^{D}\mu^{2}/(2g_{0})$ due to Eqs. (16) and (23), and

[TABLE]

is the zero-temperature Gaussian grand potential encoding quantum fluctuations, while

[TABLE]

is the finite-temperature Gaussian grand potential, encoding thermal fluctuations. Finally, we explicitly write the zero-temperature equation of state, namely we calculate the pressure as the opposite of the grand potential of Eq. (25) at $T=0$:

[TABLE]

In the thermodynamic limit of $L\rightarrow\infty$, the sum over $\vec{k}$ can be rewritten as a $D$-dimensional integral in momentum space $(2\pi/L)^{D}\sum_{\vec{k}}=\int d^{D}k$, and, substituting again the Bogoliubov spectrum (24), the equation of state becomes

[TABLE]

where the integral can be calculated after the explicit choice of $\tilde{V}(k)$.

2.3 Explicit implementation for finite-range interaction.

We now provide an explicit implementation of the zero-temperature equation of state (29) for a bosonic quantum fluid of particles interacting with the finite-range effective interaction

[TABLE]

where $g_{0}=\tilde{V}(k=0)$ is the usual zero-range interaction coupling, and

[TABLE]

is the first nonzero correction in the gradient expansion of an isotropic interaction potential $V(|\vec{r}|)$. At zero temperature, we expect that the finite-range corrections to the equation of state are detectable, but small with respect to the zero-range result of Ref. salasnich2016 . By using scattering theory in two spatial dimensions, these couplings can be linked with the s-wave scattering length $a_{s}$ and the characteristic range $R$ of the real interatomic two-body interaction tononi ; astrakharchik2009 ; salasnich2017

[TABLE]

where $n$ is the number density of the system in $D=2$.

The equation of state (29) becomes, with the finite-range interaction of Eq. (30)

[TABLE]

where we define the zero temperature Gaussian pressure $P_{g}^{(0)}$ as

[TABLE]

with

[TABLE]

Since the integrand function depends only on the modulus of the momentum $|\vec{k}|$, we rewrite the integral in $P_{g}^{(0)}$ using $D$-dimensional spherical coordinates, namely

[TABLE]

where $S_{D}=2\pi^{D/2}/\Gamma[D/2]$ is the solid angle in $D$-dimensions and $\Gamma[D/2]$ is the Euler Gamma function. In order to integrate this equation, we introduce the adimensional variable $t=\hbar^{2}k^{2}\lambda/(4m\mu)$, obtaining

[TABLE]

As a consequence of the substitution of the real interatomic potential with an effective interaction, the zero-temperature Gaussian pressure $P_{g}^{(0)}$ is ultraviolet divergent. In our framework, an efficient way to regularize $P_{g}^{(0)}$ is constituted by the technique of dimensional regularization thooft1972 . The basic idea of this approach is to rewrite a diverging integral in terms of the Euler beta and gamma functions, whose integral representation for $x,y,z>0$ is given by

[TABLE]

[TABLE]

Thank to the properties $B(x,y)=\Gamma(x)\Gamma(y)/\Gamma(x+y)$ and $\Gamma[z+1]=z\,\Gamma[z]$, one can extend the domain of definition of the gamma and beta functions by analytic continuation of their arguments $x,y,z$ also to negative values, which usually appear in many physical problems. However, despite this dimensional regularization procedure can be successfully used to regularize many ultraviolet diverging integrals, in our peculiar two-dimensional case the procedure described above would lead to a result containing the gamma function evaluated for negative integer values, which is again a diverging quantity. To avoid this residual divergence, we extend the dimension of the system to the complex value $\mathcal{D}=D-\varepsilon$, and we formally perform the integration of Eq. (37). We obtain

[TABLE]

in which the wavevector $\kappa$ is introduced for dimensional reasons. Notice how in $D=2$ and for $\varepsilon=0$ the Gaussian pressure is still divergent. To regularize it, we rely on the following small-$\varepsilon$ expansion of the gamma function kleinert2001

[TABLE]

where $\Psi(n+1)$ is the digamma function and $\Psi^{\prime}(n+1)$ is its derivative. Moreover, we express the exponentiation of a generic coefficient $x^{\varepsilon}$ for $\varepsilon\to 0$ as

[TABLE]

With this recipe, the Gaussian pressure $P_{g}^{(0)}$ in $D=2$ gives

[TABLE]

where $\gamma\approx 0.55722$ is the Euler-Mascheroni constant. Finally, we delete the $o(\varepsilon^{-1})$ divergence in the square bracket zeidler2009 and we rewrite the zero-temperature equation of state $P(\mu,T=0)$ of Eq. (29) as

[TABLE]

where we define the energy cutoff $\epsilon_{0}$ as

[TABLE]

The equation of state (44) improves the one derived for bosons with a zero-range interaction popov1972 by Popov, whose result can be reproduced by setting $\lambda=1$, i.e. $g_{2}=0$. We emphasize that, with a precise tuning of the interparticle interaction (see Ref. salasnich2017 for a detailed discussion), the finite-range corrections derived within our Gaussian approximation become larger than the zero-range beyond-Gaussian ones obtained by Mora and Castin moracastin . For weakly-interacting bosons with $na_{s}^{2}\ll 1$, where $a_{s}$ is the two-dimensional s-wave scattering length, we expect that the nonuniversal corrections of Eq. (44) arise for $R\geq a_{s}$, where $R$ is the characteristic range of the interaction. In this intermediate regime the neglection of higher order terms in the gradient expansion of Eq. (30) is justified but, at the same time, the finite-range contributions are of comparable size to the zero-range ones.

3 Conclusions

In this work we derive the two-dimensional zero-temperature equation of state for a bosonic quantum fluid with a generic isotropic interaction. The superfluid perspective is emphasized by performing the Gaussian functional integration within a phase-amplitude parametrization of the complex order parameter. For a system with zero-range interaction, we reproduce the classical result by Popov. Nonetheless, we apply a novel dimensional regularization recipe to reproduce the nonuniversal corrections for a finite-range interaction potential. We expect that, with a fine tuning of the experimental interaction parameters, the finite-range correction produce sizable corrections to the thermodynamics of the weakly-interacting superfluid. Our derivation of the zero-temperature equation of state is valid also for other interparticle interactions. In particular, the previous results can be extended for a quasi-two-dimensional system of dipolar bosons whose polarization direction is perpendicular to the plane of confinement. For a generic orientation, however, it is necessary to consider the dependence of the interaction on the in-plane angle between the particles and to include it consistently in the dimensional regularization procedure.

4 Appendix

We illustrate here the procedure to calculate the summation over the bosonic Matsubara frequencies $\omega_{n}$, which are defined as

[TABLE]

where $n\in\mathbb{Z}$ are integer numbers. The most common sum that one has to perform is in the form

[TABLE]

Using the properties of the logarithm and considering that the summation involves all $n\in\mathbb{Z}$ integers, both positive and negative, $I[\xi_{\vec{k}}]$ can also be rewritten in the useful form

[TABLE]

Taking the derivative of $I[\xi_{\vec{k}}]$ with respect to $\xi_{\vec{k}}$ in the Eq. (47) we get

[TABLE]

In the zero temperature limit, the difference

[TABLE]

becomes infinitesimal and we can substitute the sum over $n$ with an integral over $\omega$, obtaining

[TABLE]

which is the zero-temperature contribution to $I[\xi_{\vec{k}}]$. If the temperature is relatively low, but non-zero, we cannot substitute the sum in Eq. (49) with an integral, but we can rewrite it as

[TABLE]

and, using the identity

[TABLE]

we obtain

[TABLE]

We integrate this equation on $\xi_{\vec{k}}$ and, setting the arbitrary constant resulting from the indefinite integral to zero (it is not dependent on physical parameters), we finally obtain the result of the summation over the Matsubara frequencies

[TABLE]

which is used in this article to obtain Eq. (21).

Acknowledgements.

The author thanks Luca Salasnich and Alberto Cappellaro for useful discussions and suggestions. \conflictofinterests“The author declares no conflict of interest.”