Effective expression of the Lee-Huang-Yang energy functional for heteronuclear mixtures

TL;DR

This paper derives an accurate approximation of the Lee-Huang-Yang energy functional for heteronuclear Bose mixtures at mean-field collapse, showing it closely resembles the homonuclear case with specific parameter conditions.

Contribution

It introduces a simplified, accurate expression for the LHY energy functional in heteronuclear mixtures, extending the homonuclear form with two fitting parameters.

Findings

The approximated LHY functional matches the homonuclear form under permutation invariance.

Parameter values for heteronuclear mixtures are identical to the homonuclear case.

Quantitative analysis of deviations from the exact LHY functional.

Abstract

We consider a homogeneous heteronuclear Bose mixture with contact interactions at the mean-field collapse, i.e. with interspecies attraction equal to the mean geometrical intraspecies repulsion. We show that the Lee-Huang-Yang (LHY) energy functional is accurately approximated by an expression that has the same functional form as in the homonuclear case. The approximated energy functional is characterized by two exponents, which can be treated as fitting parameters. We demonstrate that the values of these parameters which preserve the invariance under permutation of the two atomic species are exactly those of the homonuclear case. Deviations from the exact expression of LHY energy functional are discussed quantitatively and a specific application is described.

Tables

Table 1. Table 1: Bose mixtures experimentally realized and amenable to collapse at g 12 = − g 11 ​ g 22 subscript 𝑔 12 subscript 𝑔 11 subscript 𝑔 22 g_{12}=-\sqrt{g_{11}g_{22}} . For the purpose of comparison, we assume the same harmonic potential for all species, corresponding to an harmonic frequency of 100 100 100 Hz for 87 Rb, i.e. ω j / ( 2 ​ π ) = 100 ​ m R ​ b / m j subscript 𝜔 𝑗 2 𝜋 100 subscript 𝑚 𝑅 𝑏 subscript 𝑚 𝑗 \omega_{j}/(2\pi)=100\sqrt{m_{Rb}/m_{j}} Hz for species j 𝑗 j .

Species 1,2		$a_{11}/a_0$	$a_{22}/a_0$	$z$	${10}^{7}B$	$n_0/N^{4/13}$[cm-3]	Ref.
39K	39K	86	34	1	2.4	$4.0\cdot {10}^{13}$	Cabrera et al. (2018),Semeghini et al. (2018)
23Na	39K	52	13	1.7	0.32	$8.9\cdot {10}^{13}$	Schulze et al. (2018)
41K	87Rb	65	100.4	2.1	8.2	$3.1\cdot {10}^{13}$	D’Errico et al.
39K	87Rb	10	100.4	2.2	0.76	$7.9\cdot {10}^{13}$	Wacker et al. (2015)
23Na	87Rb	54.5	100.4	3.8	5.0	$3.0\cdot {10}^{13}$	Wang et al. (2016)

Equations

$$E_{\mathrm{LHY}}$$

$${\mathcal{E}}_{\mathrm{LHY}}$$

$$=\kappa m_1^{3/2}(g_{11}{n}_1{)}^{5/2}f(z,u,x),$$

$$f(1,u,x)=\sum _{\pm }{\left(1+x\pm \sqrt{(1-x{)}^2+4ux}\right)}^{5/2}/(4\sqrt{2})$$

$$f(z,u,x)=\frac{15}{32}{\int }_0^{\infty }{k}^2\mathcal{F}(k,z,u,x)dk,$$

$$\mathcal{F}(k,z,u,x)$$

$${\mathcal{E}}_{\mathrm{LHY}}$$

$$f_L(z,x)$$

$$f_L(z,x)\simeq 1+\frac{5}{2}x+\frac{3}{2}x(z-1)+O[z-1{]}^2.$$

$${\mathcal{E}}_{\mathrm{LHY}}=\kappa m^{3/2}{\left(gn\right)}^{5/2}.$$

$${\mathcal{E}}_{\mathrm{LHY}}=\kappa (m_1{m}_2{)}^{3/4}(g_{11}{n}_1{g}_{22}{n}_2{)}^{5/4}{z}^{-3/4}{x}^{-5/4}f(z,1,x).$$

$$f(z,u,x)=z^{3/2}{x}^{5/2}f\left(\frac{1}{z},u,\frac{1}{x}\right).$$

$$F_q(x,z)={\left[1+{\left(z^{3/2}{x}^{5/2}\right)}^{1/q}\right]}^q,$$

$$f\left(z,1,x\right)\simeq F_{5/2}(z,x)\equiv (1+z^{3/5}x{)}^{5/2}.$$

$${\mathcal{E}}_{\mathrm{LHY}}\simeq \frac{8m_1^{3/2}(g_{11}{n}_1{)}^{5/2}}{15{\pi }^2{\hslash }^3}{\left[1+{\left(\frac{m_2}{m_1}\right)}^{3/5}\frac{g_{22}{n}_2}{g_{11}{n}_1}\right]}^{5/2},$$

$$g(z,\xi )\equiv f(z,1,z^{-3/5}\xi );$$

$$g(z,\xi )=\sum _{n=0}^{+\infty }{G}_n(\xi )(z-1{)}^n,$$

$$G_0(\xi )$$

$$G_1(\xi )$$

$$G_2(\xi )$$

$$G_3(\xi )$$

$$f(z,1,x)$$

$$+\frac{4z^{3/5}x(47z^{3/5}x+43)}{525\sqrt{z^{3/5}x+1}}(z-1{)}^3+O\left((1-z{)}^4\right)$$

$$\equiv f^{(3)}(z,x)+O\left((z-1{)}^4\right)$$

$$E=\int d\mathbf{r}\left\{\sum _{i=1,2}\left[\frac{{\hslash }^2}{2m_i}|\nabla {\psi }_i{|}^2+V_i|{\psi }_i{|}^2+\frac{1}{2}{g}_{ii}|{\psi }_i{|}^4\right]+g_{12}|{\psi }_1{|}^2|{\psi }_2{|}^2+{\mathcal{E}}_{\mathrm{LHY}}\right\}$$

$$E$$

$$m^{\ast }=m_1\left[\frac{1+\sqrt{z{a}_{11}/a_{22}}}{1+\sqrt{a_{11}/z{a}_{22}}}\right].$$

$${\mathcal{E}}_{\mathrm{LHY}}$$

$$=C_{\mathrm{LHY}}^{\ast }|\psi {|}^5$$

$$C_{\mathrm{LHY}}^{\ast }\equiv \frac{256\sqrt{\pi }}{15}\frac{{\hslash }^2}{m_1}{\left[a_{11}\frac{1+z^{1/10}\sqrt{a_{22}/a_{11}}}{1+\sqrt{z{a}_{11}/a_{22}}}\right]}^{5/2}.$$

$${\xi }_{\mathrm{LHY}}^{-2}$$

$${\left[\frac{1+z^{1/10}\sqrt{a_{22}/a_{11}}}{1+\sqrt{z{a}_{11}/a_{22}}}\right]}^{5/2},$$

$$\omega /{\omega }^{\ast }=2+\frac{64\sqrt{2}{N}^{3/2}}{5\sqrt{5}{\pi }^{7/4}}{\left(\frac{a_{11}}{a_{ho}}\right)}^{5/2}{\left[\frac{1+\sqrt{z{a}_{11}/a_{22}}}{1+\sqrt{a_{11}/z{a}_{22}}}\right]}^{13/8}{\left[\frac{1+z^{1/10}\sqrt{a_{22}/a_{11}}}{1+\sqrt{z{a}_{11}/a_{22}}}\right]}^{5/2}.$$

$$n_0$$

Full text

I Introduction

Bose-Einstein condensates (BEC) with ultracold gases are so dilute that interatomic interactions are generally weak, making the consolidated framework of mean-field (MF) theory a reliable and accurate approach in many cases. However, even with ultracold gases, there are instances where quantum fluctuations become relevant, thus producing sizeable deviations from the MF theory. Going beyond the MF analysis, N. Bogoliubov derived the many-body ground state of a weakly interacting gas of bosons Bogolubov (1947), assuming that most particles occupy the same quantum state, i.e. belong to the condensed fraction, except for a small fraction (“quantum depletion”). The leading correction to the MF energy for identical bosons with hard-sphere interactions was later obtained by Lee-Huang-Yang (LHY) Lee et al. (1957) in the homogeneous case. These results were then extended to the case of two-component bosons by D. Larsen Larsen (1963).

A well-known example where quantum fluctuations are important in ultracold atomic gases is provided by optical lattices that suppress the single-particle kinetic energy and increase the interatomic interactions Jaksch et al. (1998); Greiner et al. (2002); Xu et al. (2006). Another case is when the interatomic interactions are enhanced by increasing the scattering length in proximity of a Feshbach resonance. This route was firstly followed in fermionic condensates, where accurate measurements of the frequency shifts in collective oscillations Altmeyer et al. (2007) and of density profiles Shin et al. (2008); Navon et al. (2010) showed deviations from the MF predictions. Similar effects were later observed in strongly interacting bosonic gases in the excitation spectrum of 85Rb Papp et al. (2008) and in the equation of state of 7Li Navon et al. (2011). The momentum distribution of the quantum (and thermal) depletion has been measured in a BEC of metastable He Chang et al. (2016) and the total depletion density in a BEC of 39K Lopes et al. (2017).

Recently, D. Petrov proposed a new setting where quantum fluctuations could have striking effects, i.e. a BEC binary mixture Petrov (2015): when both components have repulsive intraspecies interactions but the interspecies interactions are attractive and sufficiently large, the MF energy may vanish or even become negative. In this case, the LHY energy correction dictates the behavior of the system, allowing for the existence of “quantum droplets”, i.e. stable states self-bound by interatomic attractive interactions. Such states exhibit the properties of a liquid even if their density is several orders of magnitude lower than in ordinary liquids Petrov (2015); Cappellaro et al. (2017); Cabrera et al. (2018); Cheiney et al. (2018); Semeghini et al. (2018); Ferioli et al. (2019); D’Errico et al. . Similarly, LHY corrections also stabilize droplets with attractive dipolar forces Ferrier-Barbut et al. (2016); Chomaz et al. (2016).

To describe this novel system, the MF Gross-Pitaevskii equation (GPE) has been modified to incorporate the extra energy due to quantum fluctuations. For mixtures, the two-component LHY energy calculated in the homogeneous case is added to the MF energy functional to obtain a generalized GPE (g-GPE), that can be also used in the inhomogeneous case within the local density approximation Cabrera et al. (2018); Cheiney et al. (2018); Semeghini et al. (2018); Ancilotto et al. (2018); D’Errico et al. . While for an homonuclear mixture of equal mass components the LHY energy functional is analytically known in a closed and simple form Petrov (2015), to date no analytic expression is known for a heteronuclear mixture and the LHY term must be calculated numerically Ancilotto et al. (2018). This implies that the LHY energy must be tabulated for a large number of values of the ratio of two components densities $n_{2}/n_{1}$, prior to the numerical integration of the g-GPE.

In this paper, guided by symmetry considerations and asymptotic behavior, we derive a simple and appealing functional form for the heteronuclear mixtures at the MF collapse, as detailed below. We verify that this function accurately approximates the numerical integral and, as such, can be efficiently used in the numerical solution of the g-GPE and to extend analytical results existing for homonuclear mixtures to the heteronuclear case.

In Section II, after summarizing known results to introduce the notation, we review the known analytical LHY energy in the special case of a largely imbalanced mixture, where the lighter species is the majority density component. In other words, this corresponds to the case of heavy impurities Larsen (1963). In Section III we propose an approximated formula for the LHY energy, which is the main result of this work, and we show its limits of applicability by comparing it with the exact numerical results of the LHY integral. In Secton IV we justify our proposed approximation by means of a Taylor expansion of the LHY integral around the homonuclear case. In Section V we show the usefulness of an analytic formula with a concrete example: we apply our analytic formula to the case of heteronuclear mixtures in the so-called LHY “quantum fluids” regime, i.e. when the MF interaction is null and the mixture is governed by the LHY term alone.

II LHY energy functional

The LHY energy of a Bose mixture can be written as

[TABLE]

where ${\cal E}_{\rm LHY}$ is the energy density given by Petrov (2015); Ancilotto et al. (2018)

[TABLE]

with $n_{i}(\bm{r})$($i=1,2$) being the density of each component, and $\kappa\equiv 8/(15\pi^{2}\hbar^{3})$, $z\equiv m_{2}/m_{1}$, $u\equiv g_{12}^{2}/(g_{11}g_{22})$, $x\equiv g_{22}n_{2}/(g_{11}n_{1})$. The interaction strengths are proportional to the scattering lengths $a_{11},a_{22},a_{12}$: $g_{ii}=4\pi\hbar^{2}a_{ii}/m_{i},g_{12}=2\pi\hbar^{2}a_{12}/m_{12}$, with $m_{12}$ the reduced mass. We note that the LHY energy density is independent of the sign of $g_{12}$, but the formation of quantum droplets occurs only for negative $g_{12}$. Therefore, as done in Ref. Petrov (2015), here we consider the case of a mixture with interspecies attractive interaction, namely $g_{12}<0$, with $g_{11}>0$ and $g_{22}>0$.

For the homonuclear case, $z=1$, the function is given in Ref. Petrov (2015):

[TABLE]

so that $f(1,1,x)=(1+x)^{5/2}$. For the heteronuclear case, $z\neq 1$, a simple expression similar to the one above is not available and approximate or numerical methods are needed.

The expression of $f(z,u,x)$ for heteronuclear mixtures is obtained from Ref. Petrov (2015); Ancilotto et al. (2018):

[TABLE]

[TABLE]

Since the dimensionless function $f(z,u,x)$ is weakly dependent on $u$, as shown in Ref. Petrov (2015), the latter is usually replaced by its value at MF collapse, where $g_{12}+\sqrt{g_{11}g_{22}}=0$, i.e. $u=1$:

[TABLE]

The expression of $f(z,1,x)$ coincides with that given in Ref. Ancilotto et al. (2018), where it is remarked that, since the integral converges due to the cancellation of divergent terms, its evaluations requires some care.

II.1 Limit of heavy impurities

Lacking a closed form like Eq. (3), an analytic expression of the LHY energy is given in Ref. Larsen (1963) for the special case of a repulsive mixture with $a_{11}=a_{12}=a_{22}>0$, where the density of the heavy component is much lower than that of the light component, i.e. $z>1$ and $x\ll 1$:

[TABLE]

Even if the above equation was derived for all positive scattering lengths, we can also apply it to the case of $-a_{12}=a_{11}=a_{22}>0$ since the LHY energy is independent of the sign of $a_{12}$, as shown by Eq. (2). Note, however, that the condition of all equal scattering lengths corresponds to $u=(1+z)^{2}/(4z)$, thus $u=1$ only for $z=1$.

In the general case, however, the numerical integration of the g-GPE requires the knowledge of the LHY term for generic values of $x$, and the above Eq. (7) is obviously insufficient. However, we will use it to confirm the validity of our approximated formula proposed below. For the purpose, it is useful to consider $f_{L}(z,x)$ around $z=1$, where it can be compared to $f(z,1,x)$:

[TABLE]

III Effective Ansatz

We introduce here the main result of our work, namely an approximate analytic expression for $f(z,1,x)$. The exact function, defined by Eqs. (4, II), can only be evaluated by points, through numerically integration of the above equations. As an example, in Fig. 1 we plot $f(41/87,1,x)$, corresponding to the case of a 41K-87Rb mixture.

When seeking for an approximate expression for $f(z,1,x)$, it is useful to recall that in the case of a single component the LHY correction takes the form Lee et al. (1957); Dalfovo et al. (1999)

[TABLE]

Thus, since for $x\to 0$ the LHY energy has to converge to the value of the single species 1, in this limit $f(z,1,x)\to 1$. In the opposite limit, i.e. for $x\gg 1$, the LHY energy converges to the value of the single species 2, ${\cal E}_{{\rm LHY}}\to\kappa m_{2}^{3/2}(g_{22}n_{2})^{5/2}$, implying $f(z,1,x)\to z^{3/2}x^{5/2}$.

Having established the above asymptotic behavior, we also recall that ${\cal E}_{{\rm LHY}}$ must be invariant under permutation of the two species. It is therefore convenient to rewrite Eq. (2) in a more symmetric form with respect to the species index:

[TABLE]

Now it is easy to see that invariance under permutation of indices $i=1,2$ implies

[TABLE]

The above properties suggest that one could approximate $f(z,1,x)$ with the following function

[TABLE]

which has a functional form that resembles that of the homonuclear case Eq. (3), except for the presence of the factor $z^{3/2}$, accounting for the mass ratio, and of the exponent $q(1/q)$, that is here considered as a fitting parameter. In the specific case shown of Fig. 1, a best fit to $f(z,1,x)$ returns $q=2.4612(2)$ (with the relative residuals being below $2\%$ for all values of $0.01<x<100$). In general the $q$ parameter returned by the fit depends on $z$, but for a broad range of $z$ values it is very close to the value $q=5/2$, which characterizes the function $f$ in the homonuclear mixture, Eq. (3) fit . Thus, we find convenient to approximate $f(z,1,x)$ with the following, simple generalization of the expression for the homonuclear case:

[TABLE]

As shown in Fig. 2, this form, in spite of its simplicity, provides a very accurate representation of the exact (numerical) values of $f(z,1,x)$, in a wide range of mass ratios. In addition, as requested, Eq. (14) is consistent with above Eq. (9) when we consider $x,|z-1|\ll 1$.

Finally, by combining Eq. (6) and Eq. (14), we arrive at the following effective expression for the LHY energy functional

[TABLE]

which constitutes one of the main results of this work.

IV Perturbative expansion

In order to quantify the deviation of the exact $f(z,1,x)$ from $F_{5/2}(z,x)$, we note that the latter is only function of $\xi\equiv z^{3/5}x$ (and not of $x$ and $z$ separately). Then we define the auxiliary function

[TABLE]

if $f(z,1,x)$ was a function of $\xi$ alone, then $g(z,\xi)$ would be independent on $z$. Although this is not the case in general, close to $z=1$, $g(z,\xi)$ is almost independent of z. To show this, we write the Taylor expansion:

[TABLE]

and compute analytically the coefficients $G_{n}$ from the integral Eq. (II)

[TABLE]

The vanishing first derivative, $G_{1}(\xi)=0$, confirms that $g(z,\xi)$ depends weakly on $z$ around $z=1$. Therefore,

[TABLE]

The above equation shows that our ansatz function $F_{5/2}(z,x)$ is the zeroth-order term of the Taylor expansion around $z=1$. A comparison between $F_{5/2}(z,x)$ and the third-order truncated Taylor series, $f^{(3)}(z,x)$, is shown in Fig. 3 for $z=41/87$. The lower panel shows that the first term, $F_{5/2}(z,x)$, approximates $f(z,1,x)$ better for small and large values of $x$, as expected since it is determined to match the asymptotic behaviours. Anyway, the maximum relative deviation is below 4%, proving that $F_{5/2}(z,x)$ suffices for a reliable estimation of the LHY integral.

V A Lee-Huang-Yang fluid with heteronuclear mixtures

Recently the so called “Lee-Huang-Yang fluid” has been proposed, consisting in a dilute quantum mixture where mutual interactions and atom numbers are tuned such that the MF interactions cancel out, and the resulting system is governed only by quantum fluctuations Jørgensen et al. (2018); bey . This system has been studied in the case of a homonuclear mixture ($m_{1}=m_{2}$) in Ref. Jørgensen et al. (2018), where explicit expressions for important quantities such as, e.g., the healing length and the frequency of monopole oscillation frequency, are obtained.

To prove the usefulness of our analytical, albeit approximate, form for the LHY energy discussed in Section III, we provide here a generalization of results derived in Ref. Jørgensen et al. (2018) to the experimentally relevant case of heteronuclear mixtures, by explicitly using our ansatz formula Eq. (14).

The total energy of a Bose-Bose binary mixture, with densities $n_{1},n_{2}$ and interacting via contact potentials, is given by

[TABLE]

where $\psi_{i}$ and $V_{i}$ are the wavefunction and the external potential of the $i-$th component ($i=1,2$), respectively, and $\cal{E}_{\rm LHY}$ is the LHY energy density contribution Eq. (2). We assume in the following a symmetric harmonic trap, $V_{1}({\bf r})=V_{2}({\bf r})=m_{1}\omega_{1}^{2}r^{2}/2\equiv V_{ho}$. The associated harmonic length is $a_{ho}=\sqrt{\hbar/m_{1}\omega_{1}}$.

The atom numbers and interactions can be tuned to cancel completely the mean-field terms in the above energy functional. In fact, by choosing $g_{12}=-\sqrt{g_{11}g_{22}}$ (namely, the MF collapse condition), and atom numbers such that $N_{2}/N_{1}=\sqrt{g_{11}/g_{22}}$, only the LHY energy term remains in Eq. (23), where the two component wavefunctions satisfy the condition $\psi_{2}=\psi_{1}(g_{11}/g_{22})^{1/4}$.

By defining $|\psi|^{2}=|\psi_{1}|^{2}+|\psi_{2}|^{2}$ and $\alpha=\sqrt{g_{11}/g_{22}}=z^{1/2}\sqrt{a_{11}/a_{22}}$, one has $|\psi_{1}|^{2}=(1+\alpha)^{-1}|\psi|^{2}$ and $|\psi_{2}|^{2}=\alpha(1+\alpha)^{-1}|\psi|^{2}$. Thus the total energy functional Eq. (23) reduces to:

[TABLE]

where

[TABLE]

Then, by approximating $\cal{E}_{\rm LHY}$ by the effective expression in Eq. (15), we can write

[TABLE]

with

[TABLE]

The expression in Eq. (26) makes the above Eq. (24) formally identical to the total energy functional for the homonuclear LHY fluid of Ref. Jørgensen et al. (2018), the only difference being that the coefficient $256\sqrt{\pi}\hbar^{2}|a_{12}|^{5/2}/(15m)\equiv C_{{\rm LHY}}$ is replaced by $C^{*}_{{\rm LHY}}$, $m$ by $m^{*}$, and $\omega_{1}$ by $\omega^{*}\equiv\sqrt{m_{1}/m^{*}}\,\omega_{1}$.

With the above replacements the results of Ref. Jørgensen et al. (2018) can be easily generalized to the case of different masses, $z\neq 1$. For instance, we readily obtain the healing length in the homogeneous case, $V_{ho}(\mathbf{r})=0$,

[TABLE]

and the frequency $\omega$ of monopole collective oscillation in the weakly interacting limit, i.e. $N^{3/2}|a_{12}/a_{ho}|^{5/2}\ll 1$ (where $N=N_{1}+N_{2}$ is the total atom number):

[TABLE]

As expected, in the case of equal masses $z=1$ the above expressions reduce to the corresponding expressions (10) and (15) of Ref. Jørgensen et al. (2018).

In order to estimate the effect of mass imbalances on the properties of the LHY fluid, we evaluate the deviation of the monopole oscillation frequency from the case of ideal gases, defined as $(\omega/\omega^{*}-2)/N^{3/2}\equiv B$, for various heteronuclear Bose mixtures which are currently realized in experiments of tunable dual-species condensates, as summarized in Table 1. Experimentally, the frequencies of collective modes oscillations are measured with high accuracy and have been used throughout as a powerful tool to study many-body dynamics. With typical atom numbers of $N\sim 10^{4}$, the deviations $\omega/\omega^{*}-2$ should be readily measurable.

In the Thomas-Fermi regime, i.e. the opposite limit $N^{3/2}|a_{12}/a_{ho}|^{5/2}\gg 1$, we obtain the total peak density from Eq. (9) of Ref. Jørgensen et al. (2018):

[TABLE]

where $A\simeq 1.815$. This quantity is relevant for the lifetime of the atomic LHY quantum fluid that is limited by 3-body inelastic collisions, whose rate is proportional to $n_{0}^{2}$. In Table 1 we report also the peak densities evaluated for the different Bose mixtures.

VI Conclusions

We have shown that the Lee-Huang-Yang energy functional for a heteronuclear Bose mixture at the MF collapse can be accurately approximated by an expression that has a functional form similar to the simple expression of the homonuclear case, see Eq. (15). Two different ansatz functions, with one or two exponents as fitting parameters, have been considered fit . Remarkably, we found that the values of these fitting parameters which preserve the invariance under permutation of the two species are exactly those of the homonuclear case. A quantitative analysis of the deviations from the exact expression of LHY energy functional indicates that the simple expression we propose is effective in a wide range of mass ratios, and may be useful for describing current and future experiments. Our formula greatly simplifies the numerical integration of the generalized GPE’s and allows further analytic study, as demonstrated in the case of a Lee-Huang-Yang fluid with heteronuclear mixtures.

Acknowledgements.

We are grateful to M. Prevedelli and L. Salasnich for comments and suggestions, and to L. Marmugi for carefully reading the manuscript. We acknowledge support by the Spanish Ministry of Science, Innovation and Universities and the European Regional Development Fund FEDER through Grant No. PGC2018-101355-B-I00 (MCIU/AEI/FEDER,UE), the Basque Government through Grant No. IT986-16, and the Fondazione CR Firenze through project “SUPERACI-Superfluid Atomic Circuits”.