Excitations and stability of weakly interacting Bose gases with multi-body interactions

TL;DR

This paper analyzes the excitations and stability of weakly interacting Bose gases with multi-body interactions, exploring contact and non-local cases, and applying findings to realistic systems like dipolar and Rydberg-dressed atoms.

Contribution

It provides a comprehensive analysis of multi-body interactions in Bose gases using Bogoliubov approximation, including contact, non-local, and realistic system applications.

Findings

Quasi-particle spectra and stability criteria derived for multi-body contact interactions.

Effects of non-local 3-body interactions on Bose gas excitations analyzed.

Impacts of 3-body terms in dipolar and Rydberg-dressed systems studied.

Abstract

We consider weakly interacting bosonic gases with local and non-local multi-body interactions. By using the Bogoliubov approximation, we first investigate contact interactions, studying the case in which the interparticle potential can be written as a sum of N-body {\delta}-interactions, and then considering general contact potentials. Results for the quasi-particle spectrum and the stability are presented. We then examine non-local interactions, focusing on two different cases of 3-body non-local interactions. Our results are used for systems with 2- and 3-body {\delta}-interactions and applied for realistic values of the trap parameters. Finally, the effect of conservative 3-body terms in dipolar systems and soft-core potentials (that can be simulated with Rydberg dressed atoms) is also studied.

Equations

$$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2m}{\nabla }^2\psi +V_{ext}\psi +f(\rho )\psi \text{mbox},$$

$$\displaystyle\hat{H}=\int\!\mathrm{d}\mathbf{r}\,\hat{\Psi}^{\dagger}(\mathbf{r})\,\biggl{(}-\frac{\hbar^{2}\nabla^{2}}{2m}\biggr{)}\,\hat{\Psi}(\mathbf{r})+\frac{1}{N!}\int\!\mathrm{d}\mathbf{r}_{1}\cdots\mathrm{d}\mathbf{r}_{N}\,\hat{\Psi}^{\dagger}(\mathbf{r}_{1})\cdots\hat{\Psi}^{\dagger}(\mathbf{r}_{N})\,U(\mathbf{r}_{1},\ldots,\mathbf{r}_{N})\,\hat{\Psi}(\mathbf{r}_{N})\cdots\hat{\Psi}(\mathbf{r}_{1})\mbox{,}$$

$$U({\mathbf{r}}_1,\dots ,{\mathbf{r}}_N)\equiv U_N\prod _{\begin{array}{c}i<j\end{array}}\delta ({\mathbf{r}}_i-{\mathbf{r}}_j)\text{mbox},$$

$$\hat{H}=\int\!\mathrm{d}\mathbf{r}\,\hat{\Psi}^{\dagger}(\mathbf{r})\,\biggl{(}-\frac{\hbar^{2}\nabla^{2}}{2m}\biggr{)}\,\hat{\Psi}(\mathbf{r})+\frac{U_{N}}{N!}\int\!\mathrm{d}\mathbf{r}\,\biggl{(}\hat{\Psi}^{\dagger}(\mathbf{r})\biggr{)}^{\!N}\!\!\!\cdot\!\biggl{(}\hat{\Psi}(\mathbf{r})\biggr{)}^{\!N}\mbox{.}$$

$$\hat{\Psi }(\mathbf{r})=\sum _{\mathbf{p}}{\Psi }_{\mathbf{p}}(\mathbf{r}){\hat{a}}_{\mathbf{p}}\text{mbox},$$

$$\hat{H}=\widehat{H_0}+{\hat{H}}_I\text{mbox},$$

$${\hat{H}}_0=\sum _{\mathbf{p}}{ϵ}_p^0{\hat{a}}_{\mathbf{p}}^{†}{\hat{a}}_{\mathbf{p}}\text{mbox},$$

$$\hat{H}_{I}=\frac{U_{N}}{\Omega^{N-1}N!}\!\sum_{\bigl{(}\sum_{i}\mathbf{p}_{i}=\sum_{i}\mathbf{p}_{i}^{\prime}\bigr{)}}\!\!\!\!\!\hat{a}_{\mathbf{p}_{1}^{\prime}}^{\dagger}\cdots\hat{a}_{\mathbf{p}_{N}^{\prime}}^{\dagger}\hat{a}_{\mathbf{p}_{N}}\cdots\hat{a}_{\mathbf{p}_{1}}\mbox{,}$$

$$\sum_{\bigl{(}\sum_{i}\mathbf{p}_{i}=\sum_{i}\mathbf{p}_{i}^{\prime}\bigr{)}}\!\equiv\sum_{\begin{subarray}{c}\mathbf{p}_{1},\ldots,\mathbf{p}_{N}\\ \mathbf{p}_{1}^{\prime},\ldots,\mathbf{p}_{N}^{\prime}\end{subarray}}\!\!\delta_{\mathbf{p}_{1}+\ldots+\mathbf{p}_{N},\,\mathbf{p}_{1}^{\prime}+\ldots+\mathbf{p}_{N}^{\prime}}\mbox{.}$$

$${\hat{a}}_0\sim \sqrt{N_0}\text{mbox}.$$

$$\hat{H}=\sum_{\mathbf{p}\neq 0}\epsilon_{p}^{0}\,\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+\frac{U_{N}}{\Omega^{N-1}N!}\bigl{(}\sqrt{N_{0}}\bigr{)}^{2N}\!+\frac{U_{N}}{\Omega^{N-1}N!}\bigl{(}\sqrt{N_{0}}\bigr{)}^{2N-2}\sum_{\mathbf{p}\neq 0}\biggl{\{}\binom{N}{2}\bigl{(}\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}+\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}\bigr{)}+N^{2}\,\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}\biggr{\}}\mbox{.}$$

$$\hat{H}=\frac{U_{N}N_{0}^{N}}{\Omega^{N-1}N!}+\sum_{\mathbf{p}\neq 0}\biggl{[}\epsilon_{p}^{0}+U_{N}n_{0}^{N-1}\frac{N^{2}}{N!}\biggr{]}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+U_{N}n_{0}^{N-1}\frac{N(N-1)}{2N!}\sum_{\mathbf{p}\neq 0}\biggl{[}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}+\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}\biggr{]}\mbox{.}$$

$$\hat{N}=\sum _{\mathbf{p}}{\hat{a}}_{\mathbf{p}}^{†}{\hat{a}}_{\mathbf{p}}\text{mbox},$$

$$\hat{H}=\frac{U_{N}n^{N-1}}{N!}N_{T}+\sum_{\mathbf{p}\neq 0}\Biggl{(}\epsilon_{p}^{0}+\frac{N(N\!-\!1)}{N!}U_{N}n_{0}^{N-1}\Biggr{)}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+\sum_{\mathbf{p}\neq 0}\frac{N(N\!-\!1)}{2N!}U_{N}n_{0}^{N-1}\biggl{(}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}+\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}\biggr{)}\mbox{.}$$

$$\hat{H}=\frac{U_N{N}_T^N}{{\Omega }^{N-1}N!}+\sum _{\begin{array}{c}\mathbf{p}\ne 0\\ (\mathbf{p}>0)\end{array}}{\hat{H}}_{\mathbf{p}}\text{mbox},$$

$$\hat{H}_{\mathbf{p}}=\biggl{(}\!\epsilon_{p}^{0}\!+\!\frac{N(N\!-\!1)}{N!}U_{N}n_{0}^{N\!-\!1}\biggr{)}\!\bigl{(}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+\hat{a}_{-\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}\bigr{)}+\frac{N(N\!-\!1)}{N!}U_{N}n_{0}^{N-1}\bigl{(}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}+\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}\bigr{)}\mbox{.}$$

$${\hat{a}}_{\mathbf{p}}$$

$${\hat{a}}_{-\mathbf{p}}$$

$$\hat{H}=\frac{U_{N}}{N!}n^{N-1}N_{T}+\sum_{\mathbf{p}\neq 0}\epsilon_{p}\,\hat{\alpha}_{\mathbf{p}}^{\dagger}\hat{\alpha}_{\mathbf{p}}-\frac{1}{2}\sum_{\mathbf{p}\neq 0}\biggl{(}\epsilon_{p}^{0}+X^{(N)}-\epsilon_{p}\biggr{)}\mbox{,}$$

$$X^{(N)}\equiv \frac{N(N-1)}{N!}{U}_N{n}_0^{N-1}\text{mbox},$$

$${ϵ}_p=\sqrt{({ϵ}_p^0{)}^2+2X^{(N)}{ϵ}_p^0}\text{mbox}.$$

$$\hat{H}=\int\!\mathrm{d}\mathbf{r}\,\hat{\Psi}^{\dagger}(\mathbf{r})\,\biggl{(}-\frac{\hbar^{2}\nabla^{2}}{2m}\biggr{)}\,\hat{\Psi}(\mathbf{r})+\sum_{\ell=2}^{N}\frac{U_{\ell}}{\ell!}\int\!\mathrm{d}\mathbf{r}\,\biggl{(}\hat{\Psi}^{\dagger}(\mathbf{r})\biggr{)}^{\!\ell}\!\!\!\cdot\!\biggl{(}\hat{\Psi}(\mathbf{r})\biggr{)}^{\!\ell}\mbox{,}$$

$$\hat{H}_{I}=\sum_{\ell=2}^{N}\frac{U_{\ell}}{\ell!}\frac{N_{0}^{\ell}}{\Omega^{\ell-1}}+\sum_{\ell=2}^{N}\sum_{\mathbf{p}\neq 0}\frac{U_{\ell}}{\ell!}\frac{N_{0}^{\ell-1}}{\Omega^{\ell-1}}\biggl{\{}\binom{\ell}{2}\bigl{(}\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}+\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}\bigr{)}+\ell^{2}\,\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}\biggr{\}}\mbox{.}$$

$$\hat{H}=\sum_{\ell=2}^{N}\frac{U_{\ell}}{\ell!}\frac{N_{T}^{\ell}}{\Omega^{\ell-1}}+\sum_{\mathbf{p}\neq 0}\epsilon_{p}^{0}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+\sum_{\ell=2}^{N}\sum_{\mathbf{p}\neq 0}\Biggr{\{}\frac{\ell(\ell-1)}{\ell!}U_{\ell}\,n_{0}^{\ell-1}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+\frac{\ell(\ell-1)}{2{\ell!}}U_{\ell}\,n_{0}^{\ell-1}\Biggl{(}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}+\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}\Biggr{)}\Biggr{\}}\mbox{.}$$

$$\hat{H}=\sum_{\ell=2}^{N}\frac{U_{\ell}}{\ell!}\frac{N_{T}^{\ell}}{\Omega^{\ell-1}}+\sum_{\begin{subarray}{c}\mathbf{p}\neq 0\\ (\mathbf{p}>0)\end{subarray}}\biggl{\{}\!\biggl{(}\epsilon_{p}^{0}+X\biggr{)}\biggl{(}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}+\hat{a}_{-\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}\biggr{)}+X\biggl{(}\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}+\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}\biggr{)}\!\biggr{\}}\mbox{,}$$

$$X=\sum _{\ell =2}^N\frac{\ell (\ell -1)}{\ell !}{U}_{\ell }{n}_0^{\ell -1}=\sum _{\ell =2}^N{X}^{(\ell )}\text{mbox}.$$

$$\tanh 2t=\frac{X}{{ϵ}_p^0+X}\text{mbox},$$

$$u_{p}^{2}=\frac{1}{2}\biggl{(}\frac{\xi_{p}}{\epsilon_{p}}+1\biggr{)}\mbox{,}\qquad v_{p}^{2}=\frac{1}{2}\biggl{(}\frac{\xi_{p}}{\epsilon_{p}}-1\biggr{)}\mbox{,}$$

$${ϵ}_p=\sqrt{({ϵ}_p^0{)}^2+2X{ϵ}_p^0}\text{mbox}.$$

$$\hat{H}=\sum_{\ell=2}^{N}\frac{U_{\ell}}{\ell!}\frac{N_{T}^{\ell}}{\Omega^{\ell-1}}+\sum_{\mathbf{p}\neq 0}\epsilon_{p}\,\hat{\alpha}_{\mathbf{p}}^{\dagger}\hat{\alpha}_{\mathbf{p}}-\frac{1}{2}\sum_{\mathbf{p}\neq 0}\Bigl{(}\epsilon_{p}^{0}+X-\epsilon_{p}\Bigr{)}\mbox{.}$$

$${\mathfrak{s}}^2=\frac{X}{m}\text{mbox}.$$

Full text

I Introduction

The Bogoliubov theory of weakly interacting Bose gases bogoliubov1947 provides an essential tool to investigate the effect of interactions in bosonic systems zagrebnov2001 and it plays a key role in the study of properties of Bose-Einstein condensates pethick2002 ; pitaevskii2016 . Its results for the ground-state energy and condensate depletion are in agreement in the weakly interacting limit with the findings obtained by other methods subsequently developed, including rigorous treatments lieb2005 . An important point is that also when it does not (quantitatively) work, it is useful to have results which are of guide for a qualitative understanding, as in the case of Helium$4$ annett2004 ; pitaevskii2016 , or to have estimates of the ground-state energy, as in the case of the $1D$ Lieb-Liniger model for small couplings lieb1961 . Moreover, being a self-consistent approach, where the number of condensate particles has to be self-consistently determined, it gives information on the issue whether there is condensation or not, as in low-dimensional systems pethick2002 ; pitaevskii2016 . Finally, the Bogoliubov transformation used to diagonalize the quadratic Hamiltonian obtained by the Bogoliubov approximation is used in a variety of other systems, including spin-wave theory of antiferromagnets manousakis1991 and superconductors schrieffer1964 ; bogoliubov1959 .

In this paper we study the generalization of the Bogoliubov theory to local and non-local/finite-range $3$-body and general multi-body interactions. Our reasons for such an investigation are the following:

• i)

We are firstly motivated by the the interest in studying the effects that the presence of $3$-body terms in experiments with ultracold atoms may induce on their equilibrium and dynamical properties, including the quasi-particle spectrum, with the goal to quantify how large are such effects.

• ii)

More generally, when (local or non-local/finite-range) $2$-body terms are present jointly with higher-body contributions (as $3$-body ones), they may compete to make the system stable or unstable and it is of interest to determine stability conditions and the spectrum of the quasi-particles. A typical example is given by an attractively interacting Bose gas, having $a<0$, which can be made stable by a repulsive $3$-body term.

• iii)

Another motivation is provided by the dipolar gas in presence of a $3$-body term xi16 ; bisset15 . The effect of $3$-body interaction terms in dipolar systems can be very interesting. As an example, it was recently shown that for a harmonically trapped dilute dipolar condensate with a $3$-body short-range interaction the system exhibits a condensate state and a droplet state blakie2016 , discussing how the droplet crystal may be an excited state arising from heating as the system crosses the phase transition. In the following we derive within the Bogoliubov theory the stability condition in presence of general multi-body local interactions (considering the case of $3$-body non-local interactions as well), then we discuss in detail some specific examples.

• iv)

A tool to study ultracold strongly interacting systems such as $1$D Bose gas lieb1961 or unitary Fermi gases zwerger2012 is to introduce effective local interactions. Recent examples are provided by the recent study of the monopole excitations for the $1$D Bose gas choi2015 using an effective Gross-Pitaevskii equation of the form

[TABLE]

where $V_{ext}$ is the external potential, $\rho=\mid\psi\mid^{2}$ is the density and the non-linear term $f(\rho)$ is extracted from the solution of the Bethe ansatz integral equations for the (homogeneous) $1D$ Bose gas lieb1961 ; yang1969 . Another example is the study of small-amplitude Josephson oscillations of a 6Li unitary Fermi gas in a double well potential valtolina2015 , where experimental data were compared with an equation of the form (1) with $\rho$ the pair density, $V_{ext}$ the double well one-body potential and $f(\rho)$ extracted from Monte Carlo numerical results valtolina2015 (an example of a possible parametrization of $f(\rho)$ across the BEC-BCS crossover is in manini2005 ). It is clear that in the weakly interacting limit it is $f(\rho)\propto\rho$: this corresponds for the $1D$ Bose gas to the limit $\gamma\to 0$, where $\gamma$ is the Lieb-Liniger coupling constant lieb1961 , and for fermions in the BEC-BCS crossover to the BEC limit $a\to 0^{+}$, $a$ being the scattering length. When deviations from the weakly interacting limit are incorporated through a function $f$ which is no longer proportional to $\rho$, if the function $f$ admits a series expansion of the form $f=\sum_{n}c_{n}\rho^{n}$, then there are effective multi-body local interactions (corresponding to integer values $n\geq 2$). Therefore, to treat such multi-body (albeit effective) interaction terms one needs to study in the Bogoliubov theory such higher-body terms, as we do systematically below.

The paper is organized as follows. In Section II we consider general multi-body contact interactions for a homogeneous weakly interacting gas treated in the Bogoliubov approximation. We consider first the case of an $N$-body $\delta$-interaction and successively we consider a local interaction which can be expanded in series of general $N$-body terms. For this class of interactions we compute the spectrum, the stability condition and the ground-state energy (also for the $1D$ case). We finally present results for contact interactions that cannot be expanded in series. In Section III we discuss the case of non-local interactions: after briefly reviewing the well-known case of a $2$-body non-local interaction, we study two different cases of $3$-body non-local interactions, and a comparison between these two cases is performed with a Gaussian pair-wise interaction. The results of Sections II and III are used in Section IV in the case of a model with $2$- and $3$-body contact interactions and the obtained findings are applied to possible realistic values of the trap parameters. In Section V we discuss some further realistic interaction potentials of interest for current experimental setups with ultracold atoms. We present results for a $2$-body non-local potential plus $2$- and $3$-body $\delta$-interactions, with applications to dipolar systems, e.g. magnetic atoms and polar molecules, and soft-core potentials, that can be simulated with Rydberg dressed atoms, to discuss the effect of a $3$-body interaction term. Finally we present our conclusions in Section VI, while more technical material is presented in the Appendices.

II Contact interactions

In this Section we consider general local $\delta$-interparticle potentials with multi-body interactions.

II.1 $N$-body interaction

We start by considering a model for a gas of $N_{T}$ bosonic particles interacting only via local repulsive $N$-body $\delta$-interactions.

The general Hamiltonian for $N$-body interactions reads

[TABLE]

where the bosonic field operators satisfy the canonical commutation relations $\left[\hat{\Psi}(\mathbf{r}),\hat{\Psi}(\mathbf{r^{\prime}})\right]=0=\left[\hat{\Psi}^{\dagger}(\mathbf{r}),\hat{\Psi}^{\dagger}(\mathbf{r^{\prime}})\right]$ and $\left[\hat{\Psi}(\mathbf{r}),\hat{\Psi}^{\dagger}(\mathbf{r^{\prime}})\right]=\delta(\mathbf{r}-\mathbf{r^{\prime}})$. We assume $N$-body local contact interactions having the form:

[TABLE]

where $U_{N}$ is a coefficient of dimension $[U_{N}]=[E]\cdot[L]^{3N-3}$. In the usual case of 2-body $\delta$-interaction one has $U_{2}=\frac{4\pi\hbar^{2}a}{m}$, where $m$ is the mass of the bosons pethick2002 ; pitaevskii2016 . With potential (3) the Hamiltonian reads

[TABLE]

Using for the field operator the expansion

[TABLE]

where $\Omega=L^{3}$ is the volume of the system (chosen to be a cube of side $L$ with periodic boundary conditions) and ${\Psi}_{\mathbf{p}}(\mathbf{r})=\frac{1}{\sqrt{\Omega}}e^{i\mathbf{k}\cdot\mathbf{r}}$, with $\mathbf{p}=\hbar\mathbf{k}$, the Hamiltonian assumes the form

[TABLE]

with $\hat{H_{0}}$ given as usual by

[TABLE]

where $\epsilon_{p}^{0}=p^{2}/2m$, and the interaction part reading as

[TABLE]

with the definition

[TABLE]

Now, proceeding as usual, we make the Bogoliubov prescription bogoliubov1947 ; zagrebnov2001

[TABLE]

As usual this implies that $N_{0}\sim N_{T}$, where $N_{0}$ is the condensate number, or the largest eigenvalue of the one-body density matrix pitaevskii2016 . We then proceed by neglecting in $\hat{H}_{I}$ products of 3 or more $\hat{a}_{\mathbf{p}}^{\dagger}$ with $\mathbf{p}\neq 0$.

To start with, we consider $3$-body interactions, that is, $N=3$. From (9) it follows the total momentum conservation. Therefore, one has to arrange $2$ nonzero momenta between $3$ initial and $3$ final possible momenta. To enumerate all the possibilities, we may start by considering all the momenta of the creation operators equal to zero:

[TABLE]

i.e. the contribution to $\hat{H}_{I}$ is $N_{0}^{2}\,\hat{a}_{\mathbf{p}}\hat{a}_{-\mathbf{p}}$ multiplied by $3$, that is the combinatorial multiplicity, $3=\binom{3}{2}$, as illustrated in the table above. For a generic $N$, it is therefore straightforward to conclude that the coefficient $\binom{3}{2}$ has to be replaced by $\binom{N}{2}$.

Going ahead, the next possible arrangements are

[TABLE]

from which one can infer that the corresponding contribution to $\hat{H}_{I}$ is $N_{0}^{2}\,\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{\mathbf{p}}$ with multiplicity $9$. For a general $N$ the multiplicity is $N^{2}$.

Finally, the remaining possibilities are when all the momenta of the annihilation operators are vanishing:

[TABLE]

whose contribution to the operatorial part is $N_{0}^{2}\,\hat{a}_{\mathbf{p}}^{\dagger}\hat{a}_{-\mathbf{p}}^{\dagger}$ with multiplicity 3, as in the first case considered above; hence, for a general $N$, the multiplicity is $\binom{N}{2}$.

Thus, the Hamiltonian (6) in the Bogoliubov approximation for a general $N$ reads

[TABLE]

Defining the density $n\equiv\frac{N_{T}}{\Omega}$ and the condensate fraction $n_{0}\equiv\frac{N_{0}}{\Omega}$, one gets

[TABLE]

Introducing the total particle number operator

[TABLE]

and enforcing the total number conservation (or, in other terms, subtracting the chemical potential pethick2002 ), one finally obtains

[TABLE]

We can rewrite $\hat{H}$ as

[TABLE]

where the sum on $\mathbf{p}>0$ indicates that it has to be taken over one half of momentum space, and

[TABLE]

The next step is to perform the Bogoliubov transformation

[TABLE]

where $u_{p}^{2}-v_{p}^{2}=1$ (general formulas for the coefficients $u_{p}$ and $v_{p}$ are given in the next Subsection). We obtain

[TABLE]

where the following quantity has been introduced:

[TABLE]

so that the quasi-particle spectrum is given by

[TABLE]

Of course, for interactions involving only $N$-body $\delta$-interactions the stability depends just on the sign of $X^{(N)}$: if $X^{(N)}$ is positive (negative), the argument in the square root of (19) is positive for all $\mathbf{p}$ (negative for small $\mathbf{p}$), and the system is stable (unstable). When more interactions are present, then one has to impose for stability a suitable combination of the parameters $U_{N}$ to be positive, as discussed in the next Subsection.

II.2 Sum of multi-body contact interactions

We consider in this Subsection a model where there is a sum of $2$-body, $3$-body,…,$N$-body $\delta$-interactions (where $N$ is arbitrary). In order to generalize the formulas presented in Section II.1, we consider $N_{T}$ bosons interacting via contact repulsive interactions described by the Hamiltonian

[TABLE]

where the $l$-th parameter $U_{l}$ has physical dimension $[E]\!\cdot\![L]^{3(\ell-1)}$, whose strength refers to $l$-body interaction. Again using Eq. (5), the Hamiltonian (20) reads $\hat{H}=\hat{H_{0}}+\hat{H}_{I}$ with $\hat{H_{0}}$ given by Eq. (7) and

[TABLE]

Proceeding as in Sec. (II.1), the Hamiltonian takes the form

[TABLE]

It follows

[TABLE]

where

[TABLE]

The quasi-particles are introduced according to Eqs. (16), and it has to be $u_{p}^{2}-v_{p}^{2}=1$ in order to guarantee the commutation relations \bigl{[}\hat{a}_{\pm\mathbf{p}},\hat{a}_{\pm\mathbf{p^{\prime}}}^{\dagger}\bigr{]}=\delta_{\mathbf{p}\mathbf{p^{\prime}}}. With the parametrization $u_{p}=\cosh{t}$, $v_{p}=\sinh{t}$, one gets

[TABLE]

from which follows

[TABLE]

where $\xi_{p}=\epsilon_{p}^{0}+X$ and

[TABLE]

The diagonalization yields

[TABLE]

The excitation spectrum (27) for small $p$ gives $\epsilon_{p}=\mathfrak{s}p$, with the sound velocity $\mathfrak{s}$ given by

[TABLE]

The stability condition can be deduced from the sign of $X$, stability requiring $X>0$.

II.3 Depletion at $T=0$ and ground-state energy

The density of particles in the excited states is

[TABLE]

so that the depletion fraction can be written as

[TABLE]

The previous expression is the usual one from the $2$-body contact interaction with the substitution $U_{2}n_{0}\to X$: notice however that if one wants to use it to determine self-consistently $n_{0}/n$ via the relation $n_{0}/n=1-n_{ex}/n$, one has to take into account the dependence of the coefficient $X$ (entering $\mathfrak{s}$) on the condensate density $n_{0}$ according to Eq. (24).

To compute the ground-state energy $E_{0}$ in $3D$ one has to regularize the large-$p$ divergence pethick2002 ; pitaevskii2016 . The correct way to write the ground-state energy is

[TABLE]

where we used $n_{0}\approx n$ and the sum over $\mathbf{p}$ is up to the cut-off scale pethick2002 . Finally we get:

[TABLE]

II.4 $1D$ case

The computation presented in the previous Subsections applies as well to the one-dimensional Hamiltonian

[TABLE]

(of course no finite condensate fraction is obtained in $1D$). Denoting in the $1D$ case the length of the system by $L$ and the particle density by $\rho\equiv\frac{N_{T}}{L}$, one gets the ground-state energy

[TABLE]

where $\epsilon_{p}^{0}=p^{2}/2m$, \epsilon_{p}=\sqrt{\bigl{(}\epsilon_{p}^{0}\bigr{)}^{2}+2X\epsilon_{p}^{0}}\, and

[TABLE]

In this way the ground-state energy becomes

[TABLE]

Defining

[TABLE]

after calculating the integral in Eq. (37), which converges to a finite nonzero value, the ground-state energy per particle in the Bogoliubov approximation is found to be

[TABLE]

This result is the generalization up to $N$-body contact interactions, of the Bogoliubov result obtained in lieb1961 for only $2$-body repulsive $\delta$-interactions in $1D$, which for small values of $\gamma_{2}$ is in agreement with the exact result lieb1961 .

II.5 General multi-body contact interactions

In this Subsection we briefly discuss two further generalizations of the Bogoliubov theory to two local interaction potentials.

In Section II.1 we considered a Hamiltonian of the form

[TABLE]

with $\ell=N$ integer and larger or equal than $2$. If $\ell$ is a real number (with $\ell>1$) and using in Eq. (40) $\Gamma(l+1)$ instead of $\ell!$ ($\Gamma$ being the Gamma function), then one can show that in the Bogoliubov approximation the following form for the Hamiltonian still holds:

[TABLE]

where

[TABLE]

The results of Section II.2 still hold with $X^{(\ell)}$ instead of $X$, i.e. the quasi-particle spectrum is given by $\epsilon_{p}=\sqrt{(\epsilon_{p}^{0})^{2}+2X^{(\ell)}\epsilon_{p}^{0}}$.

Finally, we may consider a general contact Hamiltonian of the form

[TABLE]

with $\hat{\rho}=\hat{\rho}(\mathbf{r})=\hat{\Psi}^{\dagger}(\mathbf{r})\hat{\Psi}(\mathbf{r})$ and ${\cal F}$ a function of the density operator, with the normal ordering to be taken in the second term of the right-hand side of (43). Details are given in Appendix A. One gets

[TABLE]

where ${\cal E}/\Omega={\cal F}\!\left({n}\right)$ and the parameter $X$ entering the quasi-particle spectrum $\epsilon_{p}=\sqrt{(\epsilon_{p}^{0})^{2}+2X\epsilon_{p}^{0}}$ given by

[TABLE]

III Non-local interactions

In the previous Section we considered general multi-body contact interactions, showing that the well-known results for the $2$-body $\delta$-interactions in the homogeneous case are generalized in the Bogoliubov approximation by the substitution $U_{2}n_{0}\to X$, where $X$ is given by Eq. (24) [or, according to the considered case, by (42) or (45)]. In principle, one should determine self-consistently $n_{0}/n$ from Eq. (30), but since the Bogoliubov approximation works when $n_{0}\approx n$, then one can make the substitution $n_{0}\to n$ in $X$, resulting for a sum of $N$-body contact interactions in the substitution $U_{2}n\to\sum_{\ell}U_{\ell}\ell(\ell-1)n^{\ell-1}/\ell!$. The case of higher-body non-local interactions is instead different and the final result (e.g., the quasi-particle spectrum) depends on the specific form of the interactions. In the following we explicitly show this for two different cases of $3$-body non-local interactions, cases that we treat after briefly recalling the corresponding well-known results for $2$-body non-local interactions.

III.1 $2$-body non-local potential

We start considering a Hamiltonian for a gas of $N_{T}$ bosons in a region of volume $\Omega$ and interacting via a $2$-body local repulsive, non-local potential $V_{2}(\mathbf{r})$:

[TABLE]

with $V_{2}(\mathbf{r})=V_{2}(-\mathbf{r})$. In momentum space, the Hamiltonian (46) reads $\hat{H}=\hat{H_{0}}+\hat{H}_{I}$ with $\hat{H_{0}}$ given by Eq. (7) and

[TABLE]

With the change of variables $\mathbf{r}=\mathbf{r}_{1}-\mathbf{r}_{2}$, $\mathbf{R}=\frac{\mathbf{r}_{1}+\mathbf{r}_{2}}{2}$ and using the Bogoliubov approximation, the Hamiltonian (46) becomes

[TABLE]

where

[TABLE]

having introduced the Fourier transform $V_{\mathbf{p}}$:

[TABLE]

The Hamiltonian (49) is readily diagonalized obtaining

[TABLE]

where the excitation spectrum is now

[TABLE]

To conclude this Section we mention that a detailed discussion of the non-local interactions in Bogoliubov approximation is reported in the recent paper yukalov16 , while a study of the two-body problem with arbitrary finite-range interactions on a lattice is in valiente2010 .

III.2 $3$-body non-local potentials

According to Eq. (2), for $3$-body interactions the Hamiltonian reads in general $\hat{H}=\hat{H}_{0}+\hat{H}_{I}$, where

[TABLE]

and

[TABLE]

In the following we consider two different kinds of non-local potentials and derive their excitation spectrum.

III.2.1 Potential as a sum of terms with $2$ factors

We consider a potential of the form

[TABLE]

where the dimension of the $V$ entering Eq. (55) is $[V]=[E]^{1/2}$ [with $V(\mathbf{r})=V(-\mathbf{r})$]. The Hamiltonian in the Bogoliubov approximation reads

[TABLE]

where we used the convention (50) for the Fourier transform (and we denote $V_{\mathbf{p}=0}$ by $V_{0}$). After some further manipulations, the final result is

[TABLE]

where

[TABLE]

with

[TABLE]

After diagonalizing (57), the quasi-particle energy spectrum $\epsilon_{\mathbf{p}}$ is seen to be

[TABLE]

III.2.2 Potential as a product of $3$ factors

The potential (55) is a sum of three terms, each of them given by the possible pairs which can be formed between the particles. One can also consider a potential which is the product of the three pair interactions:

[TABLE]

where the single factor has now dimension $[V]=[E]^{1/3}$, again $V(\mathbf{r})=V(-\mathbf{r})$. The Fourier transform $V_{\mathbf{p}}$, defined in Eq. (50), has dimension $[V_{\mathbf{p}}]=[E]^{1/3}\cdot[L]^{3}$.

The interaction Hamiltonian is written as

[TABLE]

with

[TABLE]

Using the relation (63) and performing the Bogoliubov approximation, it is found that

[TABLE]

To analyse the previous expression for $\hat{H}_{I}$, we denote by $\hat{A}_{\mathbf{p}}$ a general function of the $\hat{a}$’s operators entering in (64). In general the terms $V_{\mathbf{p}_{12}}\,V_{\mathbf{p}_{12}-\mathbf{p}}^{2}\hat{A}_{\mathbf{p}}$ and $V_{\mathbf{p}_{12}}\,V_{\mathbf{p}_{12}+\mathbf{p}}^{2}\hat{A}_{\mathbf{p}}$ are different. Nonetheless it is possible to show that

[TABLE]

Indeed, by using the fact that $V_{\mathbf{p}}=V_{-\mathbf{p}}$ and by doing the change of variables $\tilde{\mathbf{p}}_{12}=-\mathbf{p}_{12}$ in the left hand side of (65), we can rewrite the latter as

[TABLE]

so that after relabelling the index $\tilde{\mathbf{p}}_{12}\rightarrow\mathbf{p}_{12}$, Eq. (65) is proved.

Then the interaction Hamiltonian can be finally written as

[TABLE]

where

[TABLE]

The complete Hamiltonian is therefore

[TABLE]

where

[TABLE]

The quasi-particle spectrum is therefore

[TABLE]

We derive Eq. (71) in an alternative form in Appendix B, where we extend to the $3$-body interaction the procedure followed in Sec. (III.1). In Appendix B we also show the equivalence of these two approaches.

Passing from sums to integrals in Eq. (71), we obtain

[TABLE]

III.3 A specific example of non-local potential

To see how the finite-range in $2$- and $3$-body potentials modifies the quasi-particle spectrum, we choose a specific form for it, namely a Gaussian form

[TABLE]

applying it to the three cases of $2$-body finite-range [Section III.1], $3$-body finite-range sum of three terms [Section III.2.1] and $3$-body finite-range product of three terms [Section III.2.2].

For a $2$-body finite-range potential, we put

[TABLE]

with $[V]=[E]$ (and clearly $[\kappa]=[L]^{-1}$). The Fourier transform is given by

[TABLE]

To make comparison between the $3$-body finite-range potentials we pass to dimensionless units, denoted by tildes: we set $\tilde{p}=\frac{p}{2\,\hbar\,\kappa}$, $\tilde{\epsilon}_{p}=\frac{\epsilon_{p}}{\varepsilon}$ and $\tilde{V}={\pi^{3/2}}\frac{n_{0}}{\kappa^{3}}\frac{V}{\varepsilon}$, with $\varepsilon=\frac{2\,\hbar^{2}\,\kappa^{2}}{m}$. In this way the quasi-particle spectrum (52) can be written as

[TABLE]

For the $3$-body potential given in Eq. (55) we choose

[TABLE]

where for the case in consideration $[{\cal V}]=[E]^{1/2}$. One then finds

[TABLE]

with \tilde{{\cal V}}=\frac{\pi}{3}\bigl{(}\!\frac{n_{0}}{\kappa^{3}}\!\bigr{)}^{\!{2}}\frac{{\cal V}^{2}}{\varepsilon} (and $\tilde{p}$ and $\tilde{\epsilon}_{p}$ defined as above).

For the $3$-body potential given in Eq. (61) we choose

[TABLE]

where $[{\cal V}_{0}]=[E]^{1/3}$. Eq. (72) assumes then the following simple dimensionless form:

[TABLE]

with \tilde{{\cal V}}_{0}=\frac{\pi}{3\sqrt{3}}\bigl{(}\!\frac{n_{0}}{\kappa^{3}}\!\bigr{)}^{\!{2}}\frac{V^{3}}{\varepsilon} and the same definition for $\tilde{p}$ and $\tilde{\epsilon}_{p}$.

A comparison between Eqs. (76) and (78) shows that the functional form of the quasi-particle spectrum is different between $2$- and $3$-body finite-range interactions, and the two considered $3$-body finite-range interactions give quite different results. To show the differences in the spectra with the same value of the dimensionless coupling strength, i.e. setting $\tilde{V}=\tilde{{\cal V}}=\tilde{{\cal V}}_{0}$, a plot is presented for the sake of comparison in Fig. 1.

IV $2$- and $3$-body $\delta$-interactions

As a first application of the results presented in Sections II and III we consider a model with $2$- and $3$-body contact interactions and we apply the obtained findings for realistic values of the trap parameters.

For a model with $2$- and $3$-body $\delta$-interactions, the Hamiltonian reads

[TABLE]

In the Bogoliubov approximation one finds

[TABLE]

with

[TABLE]

Note that $U_{2}$ and $U_{3}$ have dimensions $[E]\!\cdot\![L]^{3}$ and $[E]\!\cdot\![L]^{6}$, respectively. The ground-state energy is found to be

[TABLE]

with \mathfrak{s}=\sqrt{X^{(3)}/m}=\sqrt{\bigl{(}n_{0}U_{2}+n_{0}^{2}U_{3}\bigr{)}/m}.

The stability as a function of $U_{1}$ and $U_{2}/n_{0}$ is simply given by $X^{(3)}>0$, i.e. by

[TABLE]

Using (85) we can evaluate the values of $U_{3}$ for which one has instability, as depicted in Fig. 2. To make contact with a notation often used in the literature, we set $U_{3}\equiv 2g_{3}$. Theoretical estimates for $g_{3}$ have been given in literature gammal2001 ; kohler2002 ; li2010 . Regarding one-dimensional trapping potentials, we mention that the $3$-body recombination rate in $1D$ Bose gases was experimentally studied laburthe2004 ; haller2011 , and several theoretical studies addressed the problem of determining $g_{3}$ and its effects in the Lieb-Liniger model gangardt2003 ; kheruntsyan2003 ; cheianov2006 ; kormos2009 ; kormos2011 ; piroli2015 .

Finally, the depletion fraction is given by

[TABLE]

where in the right-hand side it has been used the fact that $n_{0}\approx n$.

As an application of the previous results we may consider the case of 85Rb atoms, having a negative scattering length $a_{s}<0$ cornish2000 ; dorney2001 ; cornish2006 ; altin2010 ; mcdonald2014 ; everitt2017 . In the setup described in everitt2015 the breathing frequency of a 85Rb gas was studied after a variation of $a_{s}$. To refer to a realistic setup and reminding that in our study we are not including $3$-body losses (so that $g_{3}=Re[g_{3}]$), we deal with an external potential having the form

[TABLE]

For an anisotropic trap of the form (87) a very simple estimate of the density $n_{0}$ can be obtained by setting $n_{0}\approx N_{T}/\Omega$ and choosing as effective volume $\Omega$ the product $\ell_{x}\ell_{y}\ell_{z}\equiv\Omega_{TF}$ of the Thomas-Fermi quantities $\ell_{\alpha}$ (with $\alpha=x,y,z$) defined in Appendix C, where we study the cubic-quintic Gross-Pitaevskii equation in an isotropic potential parabolic trap. Similarly to what is done for $2$-body interactions, to take into account in a simple way the effect of the anisotropy of the potential $V\left(\vec{r}\right)$, we use the formula (140) of Appendix C by substituting $\bar{\omega}=(\omega_{\perp}^{2}\omega_{z})^{1/3}$ in place of $\omega$ and $\bar{a}=\sqrt{\hbar/m\bar{\omega}}$ in place of $a=\sqrt{\hbar/m\omega}$ [where $\omega$ and $a$ refer to an isotropic potential $V\left(\vec{r}\right)=\frac{1}{2}m\omega^{2}\left(x^{2}+y^{2}+z^{2}\right)$, as studied in Appendix C].

In Fig.2, we report points corresponding to different values of $a_{s}<0$ cornish2000 for a set of realistic values of the parameters, together with result (85). For each value of $a_{s}$ we determine the Thomas-Fermi radius $R$ using Eq. (136) and then the quantities $\ell_{\alpha}$ via Eq. (138). An important point to be observed is that when $|a_{s}|$ increases the ratio $\ell_{\alpha}/R$ increases, but $R$ itself decreases (due to the small but appreciable attractive $2$-body interactions), finally resulting in a decrease of the effective volume and an increase of $U_{3}n_{0}$.

To illustrate the results reported in Fig.2 we introduced the dimensionless variable $u_{2}\equiv\frac{U_{2}}{\hbar\bar{\omega}\bar{a}^{3}}$ and $u_{3}\equiv\frac{U_{3}n_{0}}{\hbar\bar{\omega}\bar{a}^{3}}$, respectively proportional to $U_{2}$ and $U_{3}$. It emerges from the figure that the points lie in the stability region: similar results (even deeper in the stability region) would have been obtained if we had chosen the peak density $n(0)$ at the trap center, which is another reasonable choice. We also considered other values of $U_{3}$ to explore the dependence of the stability plot on the values of $U_{2}$, $U_{3}$, and the result is that, even with rather larger values of $U_{3}$, the system is stable. We observe that with (the very large value of) $a_{s}=-100a_{0}$, one gets $|u_{3}/u_{2}|\sim 10$. Furthermore, for a value $U_{3}$ which is, e.g., $10$ times larger than the one considered in Fig. 2, one would have again $|u_{3}/u_{2}|\sim 10$, again well inside the stability region.

V Other applications to ultracold atom systems

In this Section we discuss some further realistic interaction potentials which are of interest for current experimental setups with ultracold atoms. The discussed applications include dipolar systems, e.g. magnetic atoms and polar molecules, and soft-core potentials that can be simulated with Rydberg dressed atoms. For each of these cases we analyse the energy of the elementary excitations in the homogeneous limit and derive for some interesting parameter regimes the stability diagram of the Bogoliubov spectra.

Before discussing in detail these specific implementations let us consider a general model with a non-local $2$-body interaction potential accompanied by $2$- and $3$-body contact interactions. Note that it is possible to extend this model straightforwardly adding $\delta$-interactions up to $N$-body to a $2$-body non-local interaction potential, but for the sake of simplicity we will not consider this general case.

The Hamiltonian of the system reads

[TABLE]

Proceeding as in Sections II and III and expanding the quantum fields in terms of creation and annihilation operators, one arrives in the Bogoliubov approximation at the following expression for the Hamiltonian:

[TABLE]

where $X_{\mathbf{p}}=n_{0}V_{\mathbf{p}}+n_{0}U_{2}+n_{0}^{2}U_{3}$ is the sum of the Fourier components of each interaction potential times the condensate density $n_{0}$ to the proper power. We then arrive at the Bogoliubov excitation spectrum:

[TABLE]

The condensate density can be obtained by solving the self-consistent equation:

[TABLE]

An example where the stability conditions can be carried out explicitly is the Gaussian 2-body potential Eq. (74) plus a $2$- and $3$-body $\delta$-interaction. The Bogoliubov spectrum then reads

[TABLE]

with $\tilde{U}_{2}={n_{0}}\frac{U_{2}}{\varepsilon}$, $\tilde{U}_{3}={n_{0}^{2}}\frac{U_{3}}{\varepsilon}$. In Eq. (92) dimensionless units are used as in Section III.3, setting $\tilde{p}=\frac{p}{2\,\hbar\,\kappa}$, $\tilde{\epsilon}_{p}=\frac{\epsilon_{p}}{\varepsilon}$ and $\tilde{V}={\pi^{3/2}}\frac{n_{0}}{\kappa^{3}}\frac{V}{\varepsilon}$, with $\varepsilon=\frac{2\,\hbar^{2}\,\kappa^{2}}{m}$.

Setting for convenience $\tilde{U}=\tilde{U}_{2}+\tilde{U}_{3}$, we easily derive the inequality ensuring the stability for each component of the Fourier spectrum:

[TABLE]

We then arrive at the following conditions for the parameters $\tilde{V}$ and $\tilde{U}$:

[TABLE]

which define the stability regions of the Bogoliubov spectrum.

The inequalities (94) are represented in Fig. 3. Of course, for all repulsive non-local interactions and $2$- and $3$-body contact potentials the system is stable against linear perturbations around the uniform solution. When local and non-local potentials are competing (opposite sign) or when both $\tilde{U}$ and $\tilde{V}$ are negative, it is possible to find instabilities which may signal the onset of a structured ground-state configuration. In Fig. 4 we show two cases with negative short-range interaction $\tilde{U}$ where a roton-like minimum occurs. The blue dashed line spectrum is unstable to linear perturbations for finite momentum wave-vector, signalling the onset of a modulated ground-state. This spectrum lies on the separation line of Fig. 3.

V.1 $2$-body Rydberg-dressed potentials and $3$-body contact

interactions

The results presented above can be applied in the case of model potentials like soft-core interactions that can be implemented in the laboratory with Rydberg-dressed potentials to see the effect of $3$-body terms. These potentials were recently created in optical lattices with Rb atoms excited to Rydberg states biedermann16 ; Zeiher16 . The interest for such interactions is general and involves the simulation of novel kinds of spin Hamiltonians vanBijnen2015 ; Glaetzle14 for the creation of exotic phases, like the supersolid boninsegni12 ; cinti14 ; leonard16 ; li16 , and for metrological applications macri16 ; davis16 ; gil14 ; henkel10 . Motivated by these experimental results and theoretical investigations, we focus on the study of the stability diagram of a $2$-body isotropic step-potential:

[TABLE]

Physically realizable potentials generically display long-range tails, decaying generically as a power law ($\sim r^{-3}$ or $\sim r^{-6}$). However, the model potential of Eq. (95) is a good approximation of such more complicated real potentials in the sense that the many-body properties found for such potential do not differ qualitatively from the realistic ones saccani12 ; kunimi12 ; macri13 ; macri14 . It is important to recall that the Gaussian model potential described above does not fall in the same class of soft-core potentials. The reason is that the Fourier transform of the Gaussian potential never changes sign, indicating that such potential alone can never display instabilities as the ones that are found for generic soft-core potentials. In the following we analyse in more detail the $3D$ as well as the $2D$ geometries in free space for the potential (95).

V.1.1 $3$D Case

In $3D$, the Fourier transform of Eq. (95) is:

[TABLE]

where $a=4\pi R_{0}^{3}C$ and $j_{1}$ is the spherical Bessel function of the 1st kind $j_{1}(x)=\frac{\sin{x}}{x^{2}}{-}\frac{\cos{x}}{x}$. From Eq. (90) the excitation spectrum is

[TABLE]

In dimensionless form this excitation spectrum can be written as

[TABLE]

where we defined $\tilde{p}=\frac{R_{0}p}{\hbar}$, $\varepsilon=\frac{\hbar^{4}}{2mR_{0}^{2}}$, $\tilde{\epsilon}_{p}=\frac{\epsilon_{p}}{\varepsilon}$, $\tilde{C}={4\pi n_{0}R_{0}^{3}}\frac{C}{\varepsilon}$, $\tilde{U}_{2}=\frac{n_{0}U_{2}}{\varepsilon}$ and $\tilde{U}_{3}=\frac{n_{0}^{2}U_{3}}{\varepsilon}$. Defining also $\tilde{U}=\tilde{U}_{1}+\tilde{U}_{2}$, the stability condition is represented by

[TABLE]

V.1.2 $2D$ Case

The Fourier transform of Eq. (95) in $2D$ is

[TABLE]

where $J_{1}(x)$ is the Bessel function of the 1st kind. The excitation spectrum is

[TABLE]

with $\sigma_{0}\equiv\frac{N_{0}}{L^{2}}$. In dimensionless form this expression reads

[TABLE]

with $\tilde{C}={2\pi\sigma_{0}R_{0}^{2}}\frac{C}{\varepsilon}$, $\tilde{U}_{2}=\frac{\sigma_{0}U_{2}}{\varepsilon}$ and $\tilde{U}_{3}=\frac{\sigma_{0}^{2}U_{3}}{\varepsilon}$ with $\tilde{p}$, $\varepsilon$ and $\tilde{\epsilon}_{p}$ defined as above. Setting $\tilde{U}=\tilde{U}_{1}+\tilde{U}_{2}$, the stability condition is represented by

[TABLE]

In Fig. 5 we show the stability plots for the $2D$ and $3D$ geometry in terms of the regime parameters $\tilde{C}$ and $\tilde{U}$, picking up two cases for the values $(\tilde{U},\tilde{C})$, represented by a red circle and a blue square; their respective spectra are reported in Fig. 6. Notice that, contrarily to the Gaussian potential, in $2D$ competitive interactions, or attractive potentials are not necessary to have instability.

V.2 Dipolar interactions in magnetic atoms and $3$-body contact

interactions

In this Subsection we analyse the stability diagram of a homogeneous bosonic system interacting via a long-range $2$-body dipolar potential. Such potentials have been investigated in the past years for the study of effects induced by non-local interactions in the physics of BECs both in free space and in optical lattices lahaye09 ; baranov12 . One of the major problems regarding dipolar interaction in free space is their anisotropic character which induces instabilities in $3D$ homogeneous systems koch08 ; lahaye07 . On the other hand, the presence of an asymmetric harmonic trapping in combination with short-range repulsive interactions can eliminate such instabilities, opening the way to the study of interesting many-body physics with long-range interactions santos03 ; fischer06 . Recent experiments with dipolar BECs showed that under certain conditions where instability is expected from a standard Bogoliubov approach, dense clusters with many atoms can occur chomaz16 ; ferrier16 ; kadau16 ; schmitt16 ; ferrier16-1 , which are expected to be superfluid cinti16 . Two interpretations have been proposed to explain the stabilization of this phase, namely the presence of weak $3$-body interactions xi16 ; bisset15 and beyond mean-field effects (Lee-Huang-Yang type corrections) baillie16 ; bisset16 ; waechtler16 .

Motivated by these recent developments, we analyse in further detail the stability of uniform superfluids in the presence of long-range dipolar interactions and $2$- and $3$-body contact potentials. The dipolar potential can be written as:

[TABLE]

where $C_{dd}$ is the strength of the dipolar interaction and $\theta$ is the angle between the direction of polarization and the relative position of the particles. Tuning the relative angle among the particles and the quantization axis the potential can be either attractive and repulsive. The Fourier transform of Eq. (104) is

[TABLE]

where $\alpha$ is the angle between $\mathbf{p}$ and the polarization direction (of the dipole-dipole interaction).

The excitation spectrum is then found to be:

[TABLE]

with condensate density given by the equation

[TABLE]

The spectrum (106) can be expressed in terms of the ratio

[TABLE]

of the dipolar length $a_{dd}=\frac{C_{dd}\,m}{12\pi\hbar^{2}}$ to the $s$-wave scattering length, which compares the relative strength of the dipolar and contact interactions:

[TABLE]

In dimensionless units one gets

[TABLE]

where $\tilde{p}=\frac{p\,a_{dd}}{\hbar}$, $\tilde{\epsilon}_{p}=\frac{\epsilon_{p}}{\varepsilon}$, $\tilde{U}_{2}=\frac{n_{0}U_{2}}{\varepsilon}$ and $\tilde{U}_{3}=\frac{n_{0}^{2}U_{3}}{\varepsilon}$, with $\varepsilon=\frac{\hbar^{2}}{2ma_{dd}^{2}}$.

In Fig. 7 we plot the stability diagram of Eq. (110) as a function of the dimensionless parameters $\tilde{U}=\tilde{U}_{2}+\tilde{U}_{3}$ and $\tilde{C}=\varepsilon_{dd}\tilde{U}_{2}$ for $\alpha=\pi/2$. We compute the specific values of these parameters for a condensate of 164Dy also considering as the $3$-body contact interaction the value given in bisset15 (notice however that by using other values of $U_{3}$, as the ones given in waechtler16 , the effect of such terms is anyway rather small). From Fig. 7 one sees that upon varying the scattering length the uniform phase goes from a stable (blue diamond) to an unstable (red circle) configuration. Note that the $3$-body interaction enhances the stability region to values of $\varepsilon_{dd}$ larger than $1$.

VI Conclusions

In this work we have presented a systematic study of weakly interacting bosonic gases with local and non-local multi-body interactions in the Bogoliubov approximation. We considered conservative multi-body interactions for which the number of particles is conserved. In fact multi-body interactions are associated with the presence of particle losses kagan1985 , that we did not study, rather focusing on the determination of the stability conditions due to the competition between $2$- and higher-body interactions.

A variety of interparticle potentials have been considered. We first considered contact interactions, studying the case in which the interparticle potential can be written as a general sum of $N$-body $\delta$-interactions, providing the quasi-particle spectrum, the ground-state energy and the stability conditions. Results for general effective contact potentials are also presented. Our findings show that the well-known results for the $2$-body $\delta$-interactions in the homogeneous case are generalized in the Bogoliubov approximation by the substitution $U_{2}n_{0}\to X$, where $X$ is a function of the condensate fraction $n_{0}$ given by Eq. (24) for potentials which are sums of $N$-body $\delta$-potentials and by Eq. (45) in the general case. Since the Bogoliubov approximation works well when $n_{0}\approx n$, then one can make the substitution $n_{0}\to n$ in $X$, resulting for a sum of $N$-body contact interactions in the substitution $U_{2}n\to\sum_{\ell}U_{\ell}\ell(\ell-1)n^{\ell-1}/\ell!$. The case of higher-body non-local interactions is instead different from this respect and the final results depend on the specific form of the interactions. We explicitly considered two different cases of $3$-body non-local interactions.

In the last part we discussed a few interaction potentials which are of interest for current experimental setups with ultracold atoms. Implementations include systems with $2$- and $3$-body $\delta$-interactions, where we applied in the homogeneous limit for realistic values of the trap parameters. We also considered the effect of (conservative) $3$-body terms in dipolar systems, e.g. magnetic atoms and polar molecules, and soft-core potentials that can be simulated with Rydberg dressed atoms. For each of these cases we analysed the energy of the elementary excitations and derived the stability diagram of the Bogoliubov spectra for some interesting parameter regimes.

In the present paper we focused on higher-body interactions in the homogeneous limit, having in mind both $3$-body terms and general effective multi-body interactions. Of course ultracold experiments are done in confined traps, and we think that a systematic study of the Bogoliubov equations in inhomogeneous potentials with general multi-body local and non-local interactions is an interesting direction of future research.

Acknowledgements

Discussions with L. Barbiero, G. Gori and L. Salasnich are gratefully acknowledged. Useful correspondence with N. Robins is as well acknowledged. T.M. acknowledges CNPq for support through Bolsa de produtividade em Pesquisa n. 311079/2015-6 and the hospitality of the Physics Department of the University of Padova. Support form the European STREP MatterWave is acknowledged.

Appendix A Bogoliubov approximation for a general contact interaction

In this Appendix we consider the case of a general contact interaction described by a Hamiltonian of the form

[TABLE]

One has

[TABLE]

where we wrote explicitly the operator $\hat{\rho}$ in terms of the operators $\{\hat{a}_{\mathbf{p}}$, $\hat{a}_{\mathbf{p}}^{\dagger}\}$ and we used the Bogoliubov approximation. We assume that the function ${\cal F}={\cal F}(x)$ can be expanded in series up to the second order.

Since $N_{0}\approx N_{T}$ and $N_{T}\gg 1$, it follows that ${N_{0}}\gg\sqrt{N_{0}}\gg 1$, we can write:

[TABLE]

where

[TABLE]

Neglecting products of 3 or more operators, the interaction part reads

[TABLE]

where we integrated out the space variable and used the conservation of momentum, having taken explicitly the normal ordering of the operators. Using the conservation of the total number of particles, we can write

[TABLE]

so that

[TABLE]

Denoting {\cal G}(y)\equiv\frac{\partial{\cal F}\left(x\right)}{\partial x}\Big{|}_{x=y}, we see that

[TABLE]

so that the difference in the second term of the previous equation is of higher order and it can be safely neglected. Thus the complete Hamiltonian can be written as

[TABLE]

from which Eq. (44) follows.

Appendix B Excitation spectrum for a $3$-body non-local factorizable

potential

In this Appendix we consider the $3$-body finite-range potential of Eq. (61):

[TABLE]

where $V(\mathbf{r})=V(r)$. In momentum space, the Hamiltonian can be written as $\hat{H}=\hat{H_{0}}+\hat{H}_{I}$ with the kinetic term $\hat{H_{0}}$ given by Eq. (7). To write the interaction term $\hat{H}_{I}$, we consider here the change of variables {$\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3}$} $\rightarrow$ {$\mathbf{r}_{12},\mathbf{r}_{31},\mathbf{R}$}, where

[TABLE]

In the Bogoliubov approximation, after some algebra we obtain

[TABLE]

where we defined

[TABLE]

Setting

[TABLE]

proceeding as in Section III.1, one can write the Hamiltonian in the form

[TABLE]

which is diagonalizable by the standard procedure. Then the excitation spectrum is:

[TABLE]

and we obtain

[TABLE]

To demonstrate the equivalence between Eqs. (127) and (71) we have to show that their prefactors multiplying $\epsilon_{p}^{0}$ are equal. We start considering the factor of Eq. (127), which we denote by $f_{\mathbf{p}}$:

[TABLE]

Since

[TABLE]

inserting this expression in (128), it is readily seen that $f_{\mathbf{p}}$ is given by

[TABLE]

which is equal to

[TABLE]

Using in the previous expression the properties $\sum_{\mathbf{p}}V_{\mathbf{p}}=\sum_{\mathbf{p}}V_{-\mathbf{p}}$ and

[TABLE]

we retrieve the prefactor of the $\epsilon_{p}^{0}$ term in Eq. (71).

Appendix C Thomas-Fermi approximation for the cubic-quintic

Gross-Pitaevskii in an isotropic parabolic trap

In this Appendix we consider the time-independent cubic-quintic Gross-Pitaevskii equation

[TABLE]

where the parabolic trap $V(\mathbf{r})$ is assumed isotropic:

[TABLE]

To make contact with the notation used in the main text, it is $U_{2}=g_{2}$ and $g_{3}=U_{3}/2$.

In the Thomas-Fermi approximation pethick2002 ; pitaevskii2016 it is found with $g_{3}>0$ that

[TABLE]

Imposing the normalization condition $\int d\vec{r}\,|\psi|^{2}=N$ and defining the Thomas-Fermi radius $R$ such that $V(R)=\mu$, one gets

[TABLE]

where ${\cal R}$ is measured in units of the harmonic oscillator length $a=\sqrt{\hbar/m\omega}$ (with ${\cal R}\equiv R/a$) and the dimensionless quantities $A$, $B$ and $C$ in Eq. (136) are given by $A={2\pi}\sqrt{\frac{2m\omega^{2}}{g_{3}}}{a^{4}}$, $B=\frac{g_{2}^{2}}{2g_{3}\hbar\omega}$ and $C=\frac{2\pi g_{2}a^{3}}{3g_{3}}$. In (136) the function $F(x)$ is defined (with $x>0$) as

[TABLE]

The expectation value of $\langle r^{2}\rangle=\frac{\int d\vec{r}r^{2}|\psi|^{2}}{N}$ is given by

[TABLE]

where $R$ is defined by Eq. (136) and the function $G$ is given by

[TABLE]

Since $\langle r^{2}\rangle=\langle x^{2}\rangle+\langle y^{2}\rangle+\langle z^{2}\rangle$, we get $\langle x^{2}\rangle=\langle y^{2}\rangle=\langle z^{2}\rangle=\langle r^{2}\rangle/3$. To adapt the notation to that of the main text, we set $\sigma_{\alpha}\equiv\sqrt{\langle\alpha^{2}\rangle}$ and $\ell_{\alpha}\equiv 2\sigma_{\alpha}$, with $\alpha=x,y,z$, so that in the isotropic case considered in this Appendix the product $\ell_{x}\ell_{y}\ell_{z}$, which we denote by $\Omega_{TF}$, reads

[TABLE]

where $\langle r^{2}\rangle$ is given by Eq. (138).

The previous formulas simplify for $g_{2}=0$, when only the $3$-body interaction term [i.e., only the quintic term in the Gross-Pitaevskii equation (133)] is present. One gets

[TABLE]

where we introduced the dimensionless parameter ${\cal G}_{3}\equiv\frac{g_{3}}{\hbar\omega a^{6}}$ (remember that $[g_{3}]=[E]\,\cdot\,[L]^{6}$). The expectation value of $\langle r^{2}\rangle$ is simply given by

[TABLE]

Eqs. (141) and (142) have to be compared with the usual results for the cubic Gross-Pitaevskii (with $g_{3}=0$), where one has respectively

[TABLE]

(where ${\cal G}_{2}\equiv\frac{g_{2}}{\hbar\omega a^{3}}$) and

[TABLE]