Strongly correlated quantum droplets in quasi-1D dipolar Bose gas

TL;DR

This paper investigates strongly interacting dipolar bosons in a quasi-1D system, revealing the existence of quantum droplets and proposing a new density equation that captures fermionization effects, expanding understanding of quantum phase stability.

Contribution

It introduces a novel density amplitude equation incorporating Lieb-Liniger energy for strongly repelling bosons, and demonstrates the formation of quantum droplets in a quasi-1D dipolar Bose gas.

Findings

Discovery of quantum droplets in strongly coupled dipolar bosons.

Development of a new density equation accounting for fermionization.

Identification of regimes where Lee-Huang-Yang correction is not applicable.

Abstract

We exploit a few- to many-body approach to study strongly interacting dipolar bosons in the quasi-one-dimensional system. The dipoles attract each other while the short range interactions are repulsive. Solving numerically the multi-atom Schr\"{o}dinger equation, we discover that such systems can exhibit not only the well known bright soliton solutions but also novel quantum droplets for a strongly coupled case. For larger systems, basing on microscopic properties of the found few-body solution, we propose a new equation for a density amplitude of atoms. It accounts for fermionization for strongly repelling bosons by incorporating the Lieb-Liniger energy in a local density approximation and approaches the standard Gross-Pitaevskii equation (GPE) in the weakly interacting limit. Not only does such a framework provide an alternative mechanism of the droplet stability, but it also introduces means to further analyze this previously unexplored quantum phase. In the limiting strong repulsion case, yet another simple multi-atom model is proposed. We stress that the celebrated Lee-Huang-Yang term in the GPE is not applicable in this case.

Equations

$$\begin{array}{rl}\hat{H}=& -\frac{{\hslash }^2}{2m}\int dx{\hat{\psi }}^{†}(x){\nabla }^2\hat{\psi }(x)\\ & +\frac{1}{2}\int dxd{x}^{\prime }{\hat{\psi }}^{†}(x){\hat{\psi }}^{†}(x^{\prime })V_{\mathrm{eff}}(x-x^{\prime })\hat{\psi }(x)\hat{\psi }(x^{\prime })\end{array}$$

$$\begin{array}{rl}⟨\hat{H}⟩=& -\frac{1}{2}\int dx⟨{\hat{\psi }}^{†}(x){\nabla }^2\hat{\psi }(x)⟩\\ & +\frac{1}{2}\int dxd{x}^{\prime }{G}_2(x,x^{\prime })V_{\mathrm{eff}}(x-x^{\prime })\end{array}$$

$$\begin{array}{rl}E& =\int dx\left[\frac{N}{2}{\left|\nabla \psi \right|}^2+\frac{gN(N-1)}{2}\frac{|\psi {|}^6}{|\psi {|}^2+\frac{3g}{N{\pi }^2}}\right]\\ & -\frac{g_{\mathrm{dd}}N(N-1)}{2}\int dxd{x}^{\prime }{\left|\psi (x)\right|}^2{v}_{\mathrm{dd}}^{\sigma }(x-x^{\prime }){\left|\psi (x^{\prime })\right|}^2\end{array}$$

$$\begin{split}\mu&\psi(x)=-\frac{N}{2}\frac{\partial^{2}}{\partial\,x^{2}}\psi(x)+f_{\mathrm{LL}}\big{[}\psi(x)\big{]}\\ &-g_{\mathrm{dd}}N(N-1)\int dx^{\prime}v_{\mathrm{dd}}^{\sigma}(x-x^{\prime})\left|\psi(x^{\prime})\right|^{2}\psi(x),\end{split}$$

$$\hat{H}=\hat{T}+{\hat{V}}_{\mathrm{sr}}+{\hat{V}}_{\mathrm{dd}}$$

$${\hat{H}}_{ϵ}=\hat{H}-{\hat{T}}_{CM}-\widehat{ϵ}$$

$${\forall }_{|x-y|\gg d/N}⟨{\hat{\psi }}^{†}(x){\hat{\psi }}^{†}(y)\hat{\psi }(x)\hat{\psi }(y)⟩\approx ⟨{\hat{\psi }}^{†}(x)\hat{\psi }(x)⟩⟨{\hat{\psi }}^{†}(y)\hat{\psi }(y)⟩$$

$$\begin{array}{rl}\bullet & {\forall }_{x^{\prime }\in [x,x+l]}{V}_{\mathrm{dd}}(x^{\prime })\approx V_{\mathrm{dd}}(x)\\ \bullet & {\forall }_{x^{\prime }\in [x,x+l]}⟨{\hat{\psi }}^{†}(x^{\prime })\hat{\psi }(x^{\prime })⟩\approx ⟨{\hat{\psi }}^{†}(x)\hat{\psi }(x)⟩\\ \bullet & {\int }_x^{x+l}d{x}^{\prime }⟨{\hat{\psi }}^{†}(x^{\prime })\hat{\psi }(x^{\prime })⟩\gg 1\end{array}$$

$$\frac{1}{2}\int dxdy⟨{\hat{\psi }}^{†}(x){\hat{\psi }}^{†}(y)V_{\mathrm{dd}}(x-y)\hat{\psi }(x)\hat{\psi }(y)⟩\approx \frac{1}{2}\int dxdy⟨{\hat{\psi }}^{†}(x)\hat{\psi }(x)⟩V_{\mathrm{dd}}(x-y)⟨{\hat{\psi }}^{†}(y)\hat{\psi }(y)⟩$$

$$⟨\hat{H}⟩=\sum _i⟨{\hat{H}}^i⟩=\sum _i⟨{\hat{H}}_{\mathrm{loc}}^i⟩+\sum _{ij}⟨{\hat{H}}_{\mathrm{dd}}^{ij}⟩$$

$$⟨{\hat{H}}^i⟩\approx ⟨({\hat{T}}^i+{\hat{V}}_{\mathrm{sr}}^i)⟩+\frac{1}{2}\sum _j{V}_{\mathrm{dd}}((i-j)l)⟨{\hat{\psi }}^{†}(j\cdot l)\hat{\psi }(j\cdot l)⟩l.$$

$$e_{\mathrm{LL}}\left(|\psi {|}^2\right)=\frac{1}{2}{N}^2(N-1)|\psi {|}^{6}e\left(\frac{g}{N|\psi {|}^2}\right),$$

$$⟨\hat{H}⟩\approx \sum _i{e}_{\mathrm{LL}}\left(|\psi (i\cdot l){|}^2\right)l+\frac{N(N-1)}{2}\sum _{ij}|\psi (i\cdot l){|}^2{V}_{\mathrm{dd}}((i-j)l)|\psi (j\cdot l){|}^2{l}^2,\sum _i|\psi (i\cdot l){|}^{2}l\approx 1$$

$$⟨\hat{H}⟩\approx$$

$$\int |\psi (x){|}^{2}dx=1$$

$$E\left(\psi (x)\right)=$$

$$\int |\psi (x){|}^{2}dx=1$$

$$e(\gamma )=\frac{{\pi }^2}{3}\frac{\gamma }{\frac{{\pi }^2}{3}+\gamma }.$$

$$e_{\mathrm{LL}}\left(|\psi {|}^2\right)=\frac{gN(N-1)}{2}\frac{|\psi {|}^6}{|\psi {|}^2+\frac{3g}{N{\pi }^2}}..$$

$${\phi }_n(x)=\sqrt{\frac{2}{d}}\sin \left(\frac{\pi n}{d}x\right)E_n=\frac{{\hslash }^2{\pi }^2{n}^2}{2m{d}^2},\mathrm{where}n=1,2,\dots$$

$$\psi (x,y)=|{\phi }_n(x){\phi }_m(y)-{\phi }_n(y){\phi }_m(x)|/\sqrt{2}$$

$$E_{\mathrm{int}}={\int }_0^d\mathrm{d}x{\int }_0^d\mathrm{d}y{\left(\psi (x,y)\right)}^2{V}_{\mathrm{eff}}(x-y)$$

$$u$$

$$v$$

$$x=\frac{1}{2}\left(2v+u\right)y=\frac{1}{2}\left(2v-u\right)$$

$$E_{\mathrm{int}}={\int }_0^d\mathrm{d}x{\int }_0^d\mathrm{d}y{\left(\psi (x,y)\right)}^2{V}_{\mathrm{eff}}(x-y)=\int \mathrm{d}u\mathrm{d}v\tilde{\psi }(u,v{)}^2{V}_{\mathrm{eff}}(u),$$

$$\tilde{\psi }(u,v)=\psi (v+\frac{1}{2}u,v-\frac{1}{2}u)=\psi (v-\frac{1}{2}u,v+\frac{1}{2}u)=\tilde{\psi }(-u,v)$$

$$E_{\mathrm{int}}=2{\int }_0^d\mathrm{d}u{\int }_{u/2}^{d-u/2}\mathrm{d}v\tilde{\psi }(u,v{)}^2{V}_{\mathrm{eff}}(u)$$

$$E_{\mathrm{int}}$$

$$+\frac{4}{\pi dmn(m^2-n^2)}{\int }_0^d\mathrm{d}uV(u)\left(\cos \left(\frac{n\pi u}{d}\right)-\cos \left(\frac{m\pi u}{d}\right)\right)\left(m^3\sin \left(\frac{n\pi u}{d}\right)+n^3\sin \left(\frac{m\pi u}{d}\right)\right)$$

$$E_{\mathrm{KIN}}=N\sum _{k}P(k)\frac{{\hslash }^2{k}^2}{2m}=N\sum _{n=1}^N\frac{1}{N}\frac{{\hslash }^2{\pi }^2{n}^2}{2m{d}^2}=\frac{{\pi }^2{\hslash }^{2}N(N+1)(2N+1)}{12m{d}^2}$$

$$E_{\mathrm{INT}}=\sum _{ij}⟨{\hat{V}}_{ij}⟩,$$

Full text

The advent of degenerate quantum gases consisting of atoms interacting via strong long-range dipolar forces has brought a fascinating perspective of the study of new quantum states of matter. During the last three years, self-bound droplets have been unexpectedly observed Schmitt et al. (2016), so was a long-awaited roton excitation Chomaz et al. (2018) leading to the subsequent realization of dipolar supersolids Tanzi et al. (2019); Böttcher et al. (2019); Chomaz et al. (2019).

Quantum droplets in three dimensions in ultracold dipolar gases and Bose-Bose mixtures are stable due to the many-body corrections to the ground state energy that introduce an additional repulsion in the system given by a seminal Lee-Huang-Yang (LHY) term Petrov (2015). Usually negligible, it determines the properties of the system when repulsive and attractive mean-field contributions almost cancel out each other. Since the first observation of quantum droplets Kadau et al. (2016), intense theoretical Wächtler and Santos (2016); Baillie et al. (2016); Saito (2016); Ołdziejewski and Jachymski (2016); Baillie et al. (2017) and experimental Ferrier-Barbut et al. (2016); Schmitt et al. (2016); Chomaz et al. (2016); Ferrier-Barbut et al. (2018); Böttcher et al. (2019); Cabrera et al. (2018); Cheiney et al. (2018); Semeghini et al. (2018) efforts have been devoted to characterizing their properties in both dipolar gases and Bose-Bose mixtures.

Quantum liquids are also believed to be present in lower dimensions. For Bose-Bose mixtures, beyond mean-field effects are enhanced due to geometrical confinement, and as a consequence, the liquid state becomes even more ubiquitous and remarkable Petrov and Astrakharchik (2016); Astrakharchik and Malomed (2018). In this case, LHY contribution causes the additional attraction present in the system. Unfortunately, calculations of similar corrections in the dipolar gases proved to be more demanding. For a quasi-two-dimensional system, the result is non-universal and depends on features of the confinement Jachymski and Ołdziejewski (2018). For a quasi-one-dimensional dipolar gas, the attractive beyond mean-field contribution was calculated, supporting the existence of dipolar quantum droplets in such a geometry Edler et al. (2017). However, in this case, the correction comes from the possible transverse excitations of the system, which strongly suggests that this result is also non-universal.

In this Letter we present a quasi-1D quantum droplet of a different kind. Our new self-bound state emerges in the system of $N$ polarized cold dipoles with the net repulsive interaction. We focus on a regime which is far beyond the Bogoliubov approximation, and therefore in which LHY correction is not applicable. We start with a few-body system, sufficiently small to perform the numerical diagonalization. We show, that due to the interplay between short range van der Waals interaction and non-local dipolar forces, the ground state of the system is either a dipolar bright soliton (net attractive interaction) or another self-bound state, called here the quantum droplet (net repulsive interaction).

Next, we move to larger systems. As exact numerics is impossible for thousands of atoms, we change completely our framework. Referring to the energy functional discussed by Lieb Lieb and Liniger (1963) and using the local-density approximation, we propose a new equation for a density amplitude. This equation is called here the Lieb-Liniger Gross-Pitaevskii equation (LLGPE). Although our approach does not assume a universal orbital for particles, it reduces to the GPE in the limit of weak interactions or large densities. In another extreme, of the infinitely strong repulsive interaction, it was studied in Kolomeisky et al. (2000); Baizakov et al. (2009). Using our equation, we calculate the quantum phase diagram and show the characteristic flat-top shape of the density of the gas in the quantum droplet state. Finally, we use a many-body variational Ansatz to show that in the infinitely strong repulsive interaction limit the quantum droplet can be understood as the self-confined Tonks-Girardeau gas.

We consider $N$ dipolar bosons confined in both transverse directions $\hat{y}$ and $\hat{z}$ with a tight harmonic trap of a frequency $\omega_{\perp}$. Multi-particle wave-function is approximately Gaussian in tight directions for all variables. In the longitudinal direction $\hat{x}$ the space is assumed to be finite, with the length $L$. Atoms are polarized along the $\hat{x}$ axis in head-to-tail configuration Lahaye et al. (2009). Our quasi-1D system is governed by the Hamiltonian:

[TABLE]

with $\hat{\psi}(x)$ being a bosonic field operator. The effective potential consists of the long-range dipolar part and the short-range part, namely ${V_{{\rm eff}}(x)=V_{{\rm dd}}(x)+V_{{\rm sr}}(x)}$. The quasi-1D dipolar potential reads $V_{\mathrm{dd}}(x)=-\frac{\mu_{0}D^{2}}{2\pi l_{\perp}^{2}}\frac{v_{\mathrm{dd}}\left(x/l_{\perp}\right)}{l_{\perp}}$, where ${v_{\mathrm{dd}}(u)=\frac{1}{4}\left(-2|u|+\sqrt{2\pi}(1+u^{2})e^{u^{2}/2}{\rm{Erfc}}\left(|u|/\sqrt{2}\right)\right)}$ is normalized by $\int v_{\mathrm{dd}}(u)du=1$ and $l_{\perp}=\sqrt{\hbar/m\omega_{\perp}}$. Here, $D$ is a value of atomic dipole moment and $\mu_{0}$ is the vacuum permeability. This effective quasi-1D potential comes from integration of the full 3D dipolar interaction over both transverse variables Deuretzbacher et al. (2010). The singular part coming from this integration Deuretzbacher et al. (2013) is incorporated within the short range interaction.

In this Letter, the atoms repel each other on the short distance. Thus, we use the usual model of short-range interactions $V_{{\mathrm{sr}}}(x)=\frac{\hbar^{2}a}{ml_{\perp}^{2}}\delta(x)$ with $a\geq 0$ mimicking the scattering length, which can by tuned in experiments by Feschbach resonances. Below we use box units where $L$, $\hbar/L$ and $\hbar^{2}/mL^{2}$ are the units of length, momentum and energy respectively. In addition, we also define coefficients $g=\frac{\hbar^{2}a}{ml_{\perp}^{2}}$, $g_{\mathrm{dd}}=\frac{\mu_{0}D^{2}}{2\pi l_{\perp}^{2}}$, the aspect ratio $\sigma=l_{\perp}/L$ and the rescaled function $v^{\sigma}_{\mathrm{dd}}(x):=\frac{1}{\sigma}v_{\mathrm{dd}}(x/\sigma)$, so that finally the effective potential takes a compact form ${V_{\mathrm{eff}}(x)=-g_{\mathrm{dd}}v^{\sigma}_{\mathrm{dd}}(x)+g\delta(x)}$.

We are interested in how the properties of the ground state of the system depend on the strength of interactions. We manipulate two of the parameters of interactions: the ratio between dipolar and contact interactions $f_{\mathrm{dd}}=g_{\mathrm{dd}}/g$ and $g_{\mathrm{dd}}$ itself. We expect that for $f_{\mathrm{dd}}>1$ (net attraction) the ground state has negative energy and atoms form a self-bound state similar to the bright solitons studied earlier in ultracold gases. Edmonds et al. (2017). On the other hand, the bound ground state for $f_{\mathrm{dd}}<1$ (net repulsion) was also presented in quasi-1D systems within the modified mean-field analysis in the extreme case of $f_{\mathrm{dd}}=0$ (with $g_{\mathrm{dd}}>0$ and $g\rightarrow\infty$) Baizakov et al. (2009). With our methods, we can investigate features of the system in a many-body manner across the whole range of parameters. Chiefly, we are interested in the potentially transitional change of the ground state while crossing $f_{\mathrm{dd}}\sim 1$, keeping $g_{\mathrm{dd}}$ constant.

We attempt to find a few-body ground state for both $f_{\mathrm{dd}}<1$ and $f_{\mathrm{dd}}>1$ with negative energy and characteristic width smaller than $L$, so that we can study the spatial properties of such a bound-state. We access the many-body eigenstates of the Hamiltonian by diagonalization using the Lanczos algorithm Lanczos (1950). We assume the periodic boundary conditions for a numerical convenience Ołdziejewski et al. (2018, 2018). The test of convergence of our numerical method is discussed in the Supplemental Material.

As an example, in Fig. 1, we present probability density histograms for the net repulsive interactions (top panel, $f_{\mathrm{dd}}=0.9$) and the net attractive interactions (bottom panel, $f_{\mathrm{dd}}=20$). We select the interaction parameters so that the histograms have a similar width for attractive and repulsive scenarios. Both histograms are obtained by drawing particles’ positions from the many-body probability distribution with the Metropolis algorithm and aligned by rotating them such that their center of mass point in the same direction Ołdziejewski et al. (2018). We observe two spatially localized bound-states with completely different properties. Firstly, for $f_{\mathrm{dd}}=0.9$, the width increases as the number of atoms grows, which is the opposite for $f_{\mathrm{dd}}=20$. In the first case we observe local peaks whose number agrees with the number of particles, whereas in the latter only a single central peak is observed. The above features strongly resemble the quantum droplets and bright solitons differences discussed in the recent papers about dipolar systems and Bose-Bose mixture Edler et al. (2017); Astrakharchik and Malomed (2018). Therefore, we name the first case as a droplet-like and the second one as a soliton-like solution.

The origin of the differences between the two quantum phases is evident when we look at quantum bunching properties, characterized by the normally ordered second-order correlation function $G_{2}(x,x^{\prime}):=\langle\Psi^{\dagger}(x)\Psi^{\dagger}(x^{\prime})\Psi(x^{\prime})\Psi(x)\rangle$. This function is crucial for the average energy of any system described by Eq. (1):

[TABLE]

In particular, we see that in order to develop a bound-state for $f_{\mathrm{dd}}<1$, the second-order correlation function needs exhibiting a local minimum for atoms at the same place, to decrease the contribution from the short-range repulsion which requires sufficiently strong interactions.

With our diagonalization technique, we have an access to the function $G_{2}(x,x^{\prime})$ shown in the insets of Fig. 1. A dramatic difference between both situations can be found: atoms avoid each other in quantum droplets ($\frac{G_{2}(0,0)}{N(N-1)}\ll 1$) and tend to bunch in solitons ($\frac{G_{2}(0,0)}{N(N-1)}>1$) 111Note, that in our example the system is translationally invariant. Somewhat unusual normalization factor comes from the fact that we compare $G_{2}$ function with the density probability given by the histograms.).

To study larger systems around $f_{\mathrm{dd}}\sim 1$, we propose a new framework, independent from the diagonalization scheme but resulting from our previous microscopic considerations. As in the case of droplets $G_{2}(x,x^{\prime})$ exhibits the anti-bunching, one could not assume a universal orbital describing all particles as one does in the GPE derivation. Instead, the density functional approach similar to the Thomas-Fermi approximation (TFA) for fermions can be constructed in the following way SM .If only dipolar interactions change much slower than a range of anticorrelations in $G_{2}(x,x^{\prime})$, namely that $N\sigma\gg d$ with $d$ being a width (FWHM) of a solution density $\rho(x)=N\left|\psi(x)\right|^{2}$, one can treat them classically (like a trapping potential in TFA). Accordingly, we assume, that locally the atoms obey the ground state from the Lieb-Liniger model (like the kinetic energy of fermions in TFA). We approximate the Lieb-Liniger energy as a density function in a very simplified way, namely as $e_{\mathrm{LL}}=\frac{gN(N-1)}{2}\frac{|\psi|^{6}}{|\psi|^{2}+\frac{3g}{N\pi^{2}}}$ SM . With that, the energy functional for our system reads:

[TABLE]

where $\int dx\left|\psi(x)\right|^{2}=1$. Then, using variational calculus, we arrive at a new equation for $\psi(x)$:

[TABLE]

where f_{\mathrm{LL}}\big{[}\psi(x)\big{]}=\frac{\delta e_{\mathrm{LL}}}{\delta\psi^{*}} and $\mu$ is a Lagrange multiplayer.

Note that the above equation reduces to the two equations previously studied in the literature. In the limit of very weak contact interactions $g\rightarrow 0$ we restore the standard GPE, while for $g\rightarrow\infty$ (the Tonks-Girardeau limit) we obtain the equation investigated earlier in Kolomeisky et al. (2000); Baizakov et al. (2009). Therefore, we name Eq. (4) the Lieb-Liniger GPE (LLGPE). However, note that our reasoning does not originate from the usual product state approximation for the many-body wave function for ultracold bosons and does describe the density amplitude of $N$ particles rather than a universal orbital occupied by all particles. We emphasize that we use a simplified energy density functional for the ground state, which is an approximation of the full Lieb-Liniger expression, see for instance Lang et al. (2017); Yao et al. (2018).

We aim to solve the Eq. (4) for exactly the same parameters as for the few-body case. For this, we use the imaginary time evolution (ITE) technique. In Fig. 1, we compare the outcomes from the ITE with the results from the few-body calculations ($\rho(x)$ from the ITE are marked by dashed lines). Most importantly, LLGPE captures the same $N$ dependence as the few-body considerations. We see that especially for N=5 both approaches correspond to each other in a satisfying way both for soliton-like states and droplet-like states. In the first case, we also confirm sech-shape of the solution Edmonds et al. (2017), so we will call it hereafter a bright dipolar soliton. Remembering that Eq. (4) should work for a large number of atoms where $N\sigma\gg d$ rather than $N\sigma>d$ as in the few-body case, this is an excellent agreement between two different methods. It seems that the LLGPE describing macroscopic density amplitude of the many-body system reproduces few-body behaviour of the system well enough. The agreement between these two approaches gives yet another example where already not very many particles mimic the many-body limit as for instance in Wenz et al. (2013); Grining et al. (2015). Our comparison provides good reasons to focus on the features of LLGPE solutions themselves.

As we discussed before for the small system analysis, droplet-like states get wider as $N$ grows, but bright solitons shrink. Then, it is instructive to consider the first derivative of a width $d$ of the solution of Eq. (4) over $N$ as a function of the parameters $f_{\mathrm{dd}}$ and $g_{\mathrm{dd}}$. In Fig. 2 we present such an analysis. As we have expected, for $f_{\mathrm{dd}}<1$ the derivative is positive. Then, it decreases abruptly for $f_{\mathrm{dd}}>1$ 222Note that there is no sign flip of the derivative as one would expect for a bright-solitons for a purely attractive system. That would be the case for $f_{\mathrm{dd}}\rightarrow\infty$ (so that $g\rightarrow 0$).. The rapid drop of the derivative of the width is associated with a change of a density profile $\rho(x)=N\left|\psi(x)\right|^{2}$ as one crosses $f_{\mathrm{dd}}\approx 1$. We observe a transformation of the density from a flat-top to sech-shape profile.

Then, we examine the dependence of the solution on the particle number $N$. In the top panel of Fig. 3, we show the chemical potential $\mu=\frac{\partial E}{\partial N}$ as a function of $N$ for the same coupling strengths as in Fig. 2. For all cases, the total energy of the system and the chemical potential are negative, entailing that the states are self-bound. For small $N$ the chemical potential decreases for all $f_{\mathrm{dd}}$, for higher $N$ it becomes constant for $f_{\mathrm{dd}}\leq 1$ and linear for $f_{\mathrm{dd}}>1$. The constant chemical potential is a hallmark of a droplet. We observe also that as $N$ increases, the density $\rho(x)$ of droplets changes to a flat-top shape (inset in Fig. 3). Note, that even though there is no critical number of particles below which a state is not stable in the empty 1D space, one has to remember that for a very small $N$ it may happen that the solution of LLGPE does not satisfy $N\sigma\gg d$, and thus existence of the corresponding physical state of the system is not guaranteed.

Previously, we have deduced from microscopic properties of few-body systems the density amplitude equation describing self-bound quantum droplets. This equation assumes locally the Lieb-Liniger ground state. For $g\rightarrow\infty$, the ground state becomes the Tonks-Girardeau solution. Therefore, we propose a following many-body Ansatz for the droplet state for $g\rightarrow\infty$. We assume that the droplet state is a ground state of $N$ infinitely repulsive bosons in a hard-wall box of length $d$ that is a variational parameter SM . We compare our numerical minimalization findings with a solution of Eq. (4) obtained with ITE for $g\rightarrow\infty$, $\sigma=0.05$ and $g_{\mathrm{dd}}=180$ in Fig. 4. We see a good agreement between both methods for the width of a droplet. Our Ansatz reflects the structure of the droplet — the atoms stay in a compact group thanks to the dipolar attractive forces, but locally they can strongly repel (see inset of Fig. 4). The comparison between the many-body Ansatz for droplets in the Tonks-Girardeau limit and corresponding solutions of LLGPE suggests that for $g<\infty$ and $f_{\mathrm{dd}}<1$ a proper many-body Ansatz for the droplet ground state in 1D can be constructed similarly by using the Lieb-Liniger model.

In conclusion, we have found a novel dipolar quantum droplet in the regime of strongly interacting bosons for the net repulsive two-body interactions. We have also studied the droplet - bright soliton transition when changing the effective interactions in the system from being repulsive to attractive. By considering general properties of the Hamiltonian and microscopic features of the few-body system, we have argued that the existence of the 1D droplet state requires a local minimum of the normally ordered second-order correlation function $G_{2}(x,x^{\prime})$ for $x=x^{\prime}$. These requirements lead to a new equation for the density amplitude (LLGPE) that incorporates local correlations of the Lieb-Liniger model. Therefore, quantum droplets found in this Letter are far away from the usual mean-field theory even with seminal first order corrections as their effective description contains the fully quantum solution of the Lieb-Liniger model even for strong interactions. The underlying microscopic mechanism of their presence comes from the competition of short-range repulsion and non-local attraction in the system. In contrast to the recent findings with LHY-like term in Edler et al. (2017), we have found a new self-bound state for any number of atoms both in the few-body system and in LLGPE.

Our findings also motivates a similar few-body approach to strongly interacting systems in higher dimensions. It may stimulate the research to look for quantum droplets for strongly interacting systems, like electric dipoles, where the standard beyond mean-field approximations are not valid.

Supplemental material for ’Strongly correlated quantum droplets in quasi-1D dipolar Bose gas’

Appendix A Derivation of Lieb-Liniger Gross-Pitaevskii equation

We consider $N$ bosons interacting via potential $V_{{\rm eff}}(x)=V_{{\rm sr}}(x)+V_{{\rm dd}}(x)$, where $V_{{\rm sr}}(x)$ is a short-range repulsive interaction (modeled by $g\delta(x)$) and $V_{{\rm dd}}(x)$ is a long-range attractive interaction (here the dipol-dipol interaction) with characteristic range $\sigma$. The Hamiltonian can be written in general as:

[TABLE]

with $\hat{T}$ being the kinetic energy. We are interested in the regime for which $V_{{\rm dd}}(x)$ is sufficiently strong, such that the ground state of the above Hamiltonian is the bound state. In this section, we propose and justify the simplified formalism, based on applying the exact solution of the Lieb-Liniger model to the local-density approximation approach, which allows us to investigate the energy and density profile of this ground state.

First, we need to clarify the center of mass issue. As the Hamiltonian (5) is translationally invariant, therefore its ground state has no well-localized center of mass and the ground state density has to be completely flat $\langle\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\rangle={\rm const}$. It does not mean that there are no self-bound states though. This is the same problem as in the simpler case of $N$ atoms interacting only via short range attractive potential. In that case, the ground state is a bright soliton, but, as here, with an ambiguous center of mass position. As our goal is to find and study the self-bound state, we would like, in line with a typical experiment, break the symmetry. The typical way of breaking the symmetry in theoretical considerations is to mimic a measurement of few particles (with random positions $\{x_{1},...,x_{K}\}$), after which the bound state may appear in the higher order correlation function $\langle\hat{\psi}^{\dagger}(x_{K})...\hat{\psi}^{\dagger}(x_{1})\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\hat{\psi}(x_{1})...\hat{\psi}(x_{K})\rangle$ instead of $\langle\hat{\psi}^{\dagger}(x)\hat{\psi}(x)\rangle$ (see for instance Javanainen and Yoo (1996)). However, for the purpose of this derivation, we use another approach (leading to the analogous results). Instead of (5), we consider

[TABLE]

where $\hat{T}_{CM}$ is the kinetic energy of the center of mass and $\left(-\hat{\epsilon}\right)$ is an infinitesimally small, well-localized external attractive potential. Therefore, also its ground sate is well-localized and from now we will discuss the features of the ground state of $\hat{H}_{\epsilon}$ and by $\langle\cdot\rangle$ we will understand the expectation value in this state.

Let us name by $d$ the full width at the half maximum of the bound ground state. As we consider repulsive contact interaction, density distributions of any two particles are anticorrelated at small distances. It corresponds to the hole in $G_{2}(x,0)$ function around $x=0$ (see FIG. 1 in main text), whose typical width is of the order of $d/N$. For much bigger distances we have:

[TABLE]

Now, if the typical distance between neighboring particles is much smaller than the characteristic range of the non-local potential $V_{\rm dd}$, namely if $d$ satisfies $d/N\ll\sigma$, then one may find an intermediate length scale $l$, i.e. $d/N\ll l\ll\sigma$, such that:

[TABLE]

Firstly, this implies that dipolar interactions may be treated classically:

[TABLE]

Then, we discretize space into intervals of length $l$ indexed by $i$. We can rewrite Hamiltonian as a sum of a local term (short-range interactions and the kinetic energy of atoms) and a non-local term (dipolar interactions):

[TABLE]

Using our approximations we get:

[TABLE]

The last term of the above equation is constant within the $i$-th interval because we previously assumed that the dipolar potential and the density profile varies slowly. Therefore, locally, the Hamiltonian admits the Lieb-Liniger form. The only departure from the Lieb-Liniger model is the lack of the periodic boundary conditions, which does not affect the spectrum significantly. In fact, the recent experimental study of the Lieb-Liniger spectrum was performed on a gas in a waveguide Meinert et al. (2015).

All of the above provides sensible reasons to approximate the sum of the kinetic and short-range interaction energies in the $i$-th interval by the energy of the ground state of the Lieb-Liniger model Lieb and Liniger (1963) calculated for the density $N|\psi(i\cdot l)|^{2}$:

[TABLE]

where $e(\gamma)$ is a function defined in the Lieb-Liniger work Lieb and Liniger (1963). Taking sum over all $i$’s we get:

[TABLE]

which can be approximated by an integral

[TABLE]

Note, that the last approximation does not correspond to $l\rightarrow 0$; we only state that the value of sum (13) is close to the value of the integral (14) (which comes directly from assumptions (8)).

Note also, that when making the approximation $|\psi(x)|^{2}=\rm{const.}$ in every interval $x\in[i\,l,(i+1)\,l]$, we have neglected the kinetic energy connected with changing its value between neighbouring intervals. To recompense this fact, we add the kinetic energy of the envelope of the state. We end with an energy functional:

[TABLE]

The final equation resembles Gross-Pitaevskii energy functional (with additional term $e_{\mathrm{LL}}\left(|\psi(x)|^{2}\right)$). It is also worth to note, that for small $g$ it reduces to this functional (as for $g\rightarrow 0$ in local Lieb-Liniger model, the ground state tends to the product state where every particle is in the zero-momentum mode). In any regime of parameters, the functional (15) may be minimized by using analogous methods as for the Gross-Pitaevskii energy functional. Because of this mathematical similarity and to pay tribute to great physicists, we name corresponding differential equation (Eq. (4) in the main text) by the Lieb-Liniger Gross-Pitaevskii equation (LL-GPE).

We want to stress, that the only arbitrary assumption during the derivation was that the size of the bound state $d$ satisfies $d/N\ll\sigma$. Therefore, for given interaction parameters, one has to minimize (15) with respect to the density amplitude $\psi(x)$ and check if the solution satisfies this assumption. If it does, the application of the equation (15) will be justified in that case. Otherwise, the set of parameters is outside of the range of applicability of our approximation and the results should be treated with very limited confidence.

Approximated formula for Lieb-Liniger energy

To make derived functional useful for calculations, we need to have simple approximation for Lieb-Liniger energy, as there is no analytical formula for $e(\gamma)$. Therefore taking into account the following features:

• •

in the limit $\gamma\rightarrow 0$ it should agree with the first-order perturbation theory ($\frac{{\rm d}e}{{\rm d}\gamma}|_{\gamma=0}=1$)

• •

in the limit $\gamma\rightarrow\infty$ it should reproduce Tonks-Girardeau energy $\lim_{\gamma\rightarrow\infty}e(\gamma)=\frac{\pi^{2}}{3}$

• •

it should be smooth and convex in the intermediate region.

we roughly approximate it by

[TABLE]

With that, we rewrite Eq. (12) into

[TABLE]

Appendix B Many-body Ansatz in the Tonks-Girardeau limit

In the limit of infinitely strong interactions $g\rightarrow\infty$, we propose a variational Ansatz for the many-body ground state. The Ansatz assumes, that the real ground state of our system, $N$ atoms with dipolar attraction and short-range repulsion is close to the ground state of the Tonks gas, but in the box of width d. In such state one can compute the kinetic energy and interaction energy (short range and dipolar) even for hundreds of atoms and perform minimization with respect to the box width $d$. The details of the computation of the energies are presented step by step below.

Our system is characterized by:

• •

Open box - hard walls at $x=0$ and $x=d$

• •

Inifnite contact interaction strength $g\rightarrow\infty$

• •

Finite $g_{\mathrm{dd}}>0$

Our Ansatz assumes that the droplet state of such system is the ground state of infinitely repulsive bosons (fermionization) where the variational parameter $d$ depends on the dipolar interactions between the atoms. In the following subsections, we are going to introduce it in a more precise form. We use the box units as in the main text.

B.1 Single particle eigenstates in the box with the open boundary conditions

Single particle eigenstates and energies are given by:

[TABLE]

B.2 Ansatz for the two-particle system

Tonks-Girardeau Ansatz for an eigenstate with the first particle in orbital $n$ and the second in $m$:

[TABLE]

The interaction energy is given by:

[TABLE]

For arbitrary interaction potential $V(x-y)$ (in our case $V_{\mathrm{eff}}(x-y)$ from the main text), it is possible to simplify the expression by performing one integration. To do this, we will change variables:

[TABLE]

The new variables are chosen such that the surface element is unchanged, i.e. ${\rm d}x\,{\rm d}y={\rm d}u\,{\rm d}v$. The inverse transform is simply

[TABLE]

The interaction energy after the variable change will be

[TABLE]

The integration region can be simplified using the symmetries $V_{\rm{eff}}(x-y)=V_{\rm{eff}}(y-x)=V_{\rm{eff}}(u)$ and $\psi(x,y)=\psi(y,x)$. The latter implies, that the wavefunction in the variables $(u,v)$ is even with respect to $u$:

[TABLE]

Therefore the integration with respect to the whole diamond is twice the integration with respect to its part for $u>0$:

[TABLE]

The result is

[TABLE]

B.3 Tonks-Girardeau ground state with $N$ particles

Kinetic energy of $N$ atoms in the Tonks-Girardeau limit:

[TABLE]

The total interaction energy is given by

[TABLE]

where $\hat{V}_{ij}=V_{\rm{eff}}(\hat{x}_{i}-\hat{x}_{j})$ indicates the interaction potential between $i$th and $j$th particles. Due to indistinguishability the binary interaction energy is the interaction energy between the "first" and the "second" atoms multiplied by the number of pairs:

[TABLE]

The many body Ansatz for the ground state can be written as:

[TABLE]

where the $\mathcal{A}$ is antisymmetrization operator. Note, that the ground state for the box with open boundary condition is just the absolute value of the fermionic ground state. This relation does not hold for excited states. It is also not strictly true for the ground state in the box with periodic boundary conditions. The interaction energy between the two particles is expressed by:

[TABLE]

where $E_{\rm int}(m,n)$ is the contribution to interaction energy coming from two particles, one on the $m$th energy level and one on $n$th energy level. The formula for $E_{\rm int}(m,n)$ is given in Eq. (27). Combining Eq. (36) with (31) one gets:

[TABLE]

Therefore, the total energy is given by:

[TABLE]

Then, we numerically minimize Eq. (38) over the width $d$.

Appendix C Numerical diagonalization

We access the ground state of the Hamiltonian (1) by diagonalization using the Lanczos algorithm. As we are interested in investigating bound states, boundary conditions should not have significant impact on the results and therefore, we assume the periodic boundary conditions for numerical convenience. We use the Fock space spanned by the plane-wave basis. We choose the maximum single-particle momentum value which potentially may occur in calculations $2\pi k_{\rm max}$ and then we set the energy cutoff $E_{\rm cut}=(2\pi k_{\rm max})^{2}/2$, e.i. we consider only the Fock states with total the kinetic energy smaller or equal $E_{\rm cut}$. Using the fact, that the total momentum of the system is conserved (due to the periodic boundary conditions), in calculating the ground state we only need to consider the states satisfying $\sum_{k=-k_{\rm max}}^{k_{\rm max}}2\pi kN_{k}=0$, where $N_{k}$ is a number of particles in a given mode $k$.

We want to stress the fact, that to provide the same accuracy with an increasing number of particles, it is not enough to keep $k_{\rm max}$ constant – one should rather ensure that the amount of available kinetic energy per particle is the same, which means $k_{\rm max}^{2}/N=\rm{const.}$

To get some information how far from the actual ground state we are, we present in Fig. 5 the ground state energy as a function of increasing $k_{\rm max}^{2}/N$ for the droplet-like states from the main text (which was much more numerically challenging than soliton-like states). Assuming that a relative error decreases like some power of $k_{\rm max}^{2}/N$, e.i.

[TABLE]

we are able to estimate the exact value of this energy. The fact that coefficients $c$ and $p$ estimated independently for $n=3,4,5$ does not varies significantly, supports reasonability of this assumption. From Fig. 5 we see that the energy of the ground state calculated with the cutoff $k_{\rm max}^{2}/N\approx 1000$ is accurate up to $\approx 3\%$. For $N=5$ particles the number of Fock states needed for such accuracy is $4.8\cdot 10^{5}$. As a required dimension of Hilbert space increases exponentially with $N$, we are restricted to few particles calculations only.

As we did not achieve full convergence for $N=5$, one may wonder if our discussion about this state (presented in FIG. 1 of main text) is reliable. To legitimate it, in FIG. 6 we show the changes of shape of second order correlation function $G_{2}(x,0)$ with the increasing cutoff. As these changes are relatively small, we find our results accurate enough to perform qualitative analysis presented in the main text.