Abstract
We calculate analytically the quantum and thermal fluctuations corrections of a dilute quasi-two-dimensional Bose-condensed dipolar gas. We show that these fluctuations may change their character from repulsion to attraction in the density-temperature plane owing to the striking momentum dependence of the dipole–dipole interactions. The dipolar instability is halted by such unconventional beyond mean field corrections leading to the formation of a droplet phase. The equilibrium features and coherence properties exhibited by such droplets are deeply discussed. At finite temperature, we find that the equilibrium density crucially depends on the temperature and on the confinement strength and thus, a stable droplet can exist only at ultralow temperature due to the strong thermal fluctuations.
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
One of the most fascinating phenomena recently observed in dipolar Bose–Einstein condensates (BEC) with Dy and Er atoms is the formation of self-bound droplets [1–4]. Theoretically, a wealth of studies have been spawned for highlighting the behavior of droplet states [5–14]. These so-called liquid droplets are stable even in the absence of external trapping [3, 8, 9] (self-bound) due to the competition between attraction, repulsion and Lee–Huang–Yang (LHY) quantum fluctuations [15–17]. It was also found that dipolar quantum droplets are anisotropic and form a regular array, result in from the anisotropy and the long-range character of dipole–dipole interaction (DDI). These droplets, whose densities are one order of magnitude higher than the density of the ordinary BEC and decay at the droplet critical temperature at which particles undergo a phase transition from a low density gas phase (ordinary BEC) to the high density phase (droplet) [10].
In a 3D case, the LHY corrections provide a term proportional to n3/2, where n is the peak density, that arrests the dipolar instability at high condensed density [2, 5–7]. In quasi-1D dipolar BEC, the LHY quantum corrections present anomalous properties due to the transversal modes and the quantum droplet is appeared only in a low density regime [18]. Quantum fluctuations play also a crucial role in stabilizing droplets in low-dimensional nondipolar Bose–Bose mixtures [19–21].
In this paper we investigate for the first time the formation of a droplet state in quasi-2D weakly interacting dipolar bosons. We show that the system yields many surprising and interesting properties. We find that the LHY quantum corrections change their nature from repulsive to attractive due to the peculiar momentum dependence of the DDI [22], similarly to the quasi-1D droplet [18]. This unconventional behavior not only modifies the density dependence and the quantum stabilization mechanism but also unveils novel phase of matter consisting of a quantum droplet. The nucleation of this state occurs due to the competition between the roton instability results in local collapses [22], and the LHY fluctuations. It has been suggested that the roton softening combined with the quantum stabilization mechanism opens a new avenue for exploring supersolids [23].
The observed 2D self-bound dipolar droplets differ from the bound state of 2D weakly attracting bosons [24] and that of 2D Bose–Bose mixtures with both contact and dipolar interactions [19–21]. The former exists by the increased kinetic energy associated with their nonuniform shape while the latter occurs when the interspecies interaction is weakly attractive and the intraspecies ones are weakly repulsive. In addition, the transition to the droplet state happens at relatively high density compared to the quasi-1D system [18]. The density and the shape of the droplet are analyzed by numerically solving the underlying generalized nonlocal 2D Gross–Pitaevskii equation (GPE). We find that the equilibrium density is markedly affected by transversal modes. It is shown in addition that the LHY corrections invoke non-trivial enhancements in the excitations and the one-body density matrix of the droplet. Finally, we extend our study to finite temperatures and predict effects of thermal fluctuations. We determine the stability condition as well as the condensed density inside the droplet.
The rest of paper is organized as follows. Section 2 introduces the DDI in quasi-2D and the LHY quantum corrections. Section 3 deals with the stability regime, the ground-state properties and the coherence of the droplet. In section 4 we generalize our results to finite temperature. Section 5 contains our conclusions.
2. Model
2.1. Dipolar interactions in quasi-2D geometry
We consider a dilute Bose-condensed gas of dipolar bosons tightly confined in the axial direction z by an external potential
and assume that in the x, y plane the translational motion of atoms is free. The dipole moments d are oriented perpendicularly to the x, y plane. In the ultracold limit
, where
is a characteristic range of the DDI, the momentum representation of the two-body interaction potential
is given as [22]

where C = 2πd2/g,
is the 2D contact interaction coupling constant which strongly depends on the strength of the transverse confinement
, and
with a being the s-wave scattering length (a > 0 throughout the paper). Another model for the effective quasi-2D potential was proposed in [25]
, where Erfc is the complementary error function. Expanding this potential which can be obtained by integration of the full 3D dipolar interaction over the transverse harmonic oscillator at small momenta leads to V(k) = g (1 − C l0k), with l0 adjusts the scale for the strength of the linear term and can be set to unity since the ratio between the dipolar length and the trap length is the most important. Therefore, both potentials require a high momentum cut-off when calculating the beyond-mean field corrections. The large momentum behavior of both potentials is different, the potential of [25] is constant (
) for large k, while the potential (1) is linear in k. This implies different regularization schemes when computing the beyond-mean field LHY corrections.
The Bogoliubov excitation energy is given as
[22], where
and μ0 = ng is the zeroth order chemical potential. For small momenta the excitations are sound waves,
. For C varies as
, where ξ is the healing length, the excitation spectrum exhibits a roton-maxon structure [22]. The observation of such a roton mode has been reported very recently in [26] using momentum-distribution measurements in dipolar quantum gases of highly-magnetic Er atoms. For
, the uniform Bose gas becomes dynamically unstable.
2.2. LHY corrections
Indeed, obtaining reliable estimate for the beyond mean field correction to the equation of state (EoS) in low dimensions is challenging even for Bose systems with contact interactions [27–29]. At zero temperature, the LHY corrections to the EoS can be written [10, 22]

The evaluation of this integral requires special care due to the crucial contribution to the beyond mean field terms of the transverse trap modes of the contact interactions. The large-momentum divergence originating from the dipolar term −gCk (valid only for k ≪ 1/r*) is another issue of the integral (2). One possibility to solve this problem is to work with an arbitrary Λ-cutoff. In the case of contact interactions, the potential (1) takes the form V(k) = g for k < Λ, and 0 otherwise. Then, if Λ is larger than typical momenta in the gas, the obtained LHY corrections are cutoff-independent and in good agreement with the existing literature (see e.g. [27–30]). Now if one applies this method to the dipolar interaction case, it turns out that the resulting corrections to the EoS are cutoff-dependent (the cutoff is not larger than the roton momentum) due to the special character of the DDI (see e.g. [31]). Another possible route to compute the LHY corrections (2) is to take into account the full transverse structure. Obtaining reasonable stable corrections within this technique is also a tedious and time-consuming task (diagonalizing the Bogoliubov–De Gennes equations is extremely difficult both analytically and numerically) [31].
To circumvent this problem, a high-momentum cutoff is considered here which is valid in the ultracold regime k ≪ 1/r* [22]. Despite it gives qualitative correct results, it renders much simpler the calculations and captures the main features of the system at hand [22]. The choice of this momentum cutoff is not only motivated by computational convenience, but also the obtained corrections will be insensitive to the cutoff in contrast to the Λ-cutoff method. After some algebra, we obtain

where
and
. In the absence of the DDI, equation (3) excellently agrees with the usual short-range 2D Bose gas EoS (see e.g. [27, 30]). When the roton minimum is close to zero i.e.
, one has
. The quantum corrections (3) are important to halt the collapse of the system when the roton touches zero (roton instability). They can also substantially impact the collective excitations and the thermodynamics of the system.
Figure 1 clearly shows that for
,
initially increases and after it reaches its maximum at nr*2 = 0.05, starts decreasing. In this ultra-dilute limit, δμLHY provides an additional repulsive term
prohibiting the formation of any droplet in contrast to the 2D Bose–Bose mixture with contact interactions [19]. For
, the effective LHY attraction furnishes an extra term
arresting the dipolar instability, results in the formation of a stable self-bound droplet. This droplet phase has a universal peak density at
where the LHY energy,
, reaches its minimal value (see the inset of figure 1). For nr*2 > 0.65,
grows logarithmically and thus, the system undergoes an instability as the complexity increases.
Figure 1. The LHY corrections to the EoS from equation (3) as a function of nr*2 for l0/a = 40. The inset shows the LHY energy. These parameters are sufficent to reach the roton regime.
Download figure:
Standard image High-resolution image3. Quantum droplets
In this section we discuss the formation and the equilibrium properties of quasi-2D quantum droplets.
The equilibrium density n0 can be obtained by minimizing the energy per particle
with respect to the density n [19], where
. In this manner, we get

where
. Equation (4) clearly shows that the transverse harmonic confinement l0 may strongly change the equilibrium density. Therefore, the weakly interacting regime requires the condition:
or equivalently
.
To gain more insights into these quantum ensembles, we numerically solve the generalized GPE in which
:

where
is the inverse Fourier transform. The wavefunction must satisfy the normalization condition
. In 3D geometry, equation (5) has been intensively used to describe the dynamics of the droplet [5, 7–11] and already validated by quantum Monte Carlo simulations [6]. In the model (5) the effect of higher-momentum modes is just a local density-dependent term and LHY fluctuations are assumed to be large enough to maintain the overall balance with the dipolar instability. The stationary generalized GPE can be obtained from equation (5) using
, where μ is the chemical potential of the system. This yields

The numerical simulation of equation (6) was performed using the split-step Fourier transform and the convolution method to evaluate the DDI term [32]. Our simulations are carried out for
Dy atoms, with typical atom number N = 104 and a = 141a0 [33] (a0 being the Bohr radius) which can be controlled via a magnetic Feshbach. Dy atoms in their ground state have a dipolar length r* = 130a0 [23, 34]. In this configuration a roton mode is expected to appear. Figure 2 shows that a stable droplet is formed due to the competition between the roton instability and the LHY quantum fluctuations and aquires an equilibrium density n0r*2 ≃ 0.63 at which the energy develops a local minimum as is foreseen above (see the inset of figure 1). For
Dy atoms, the stability is reached at densities n0 ∼ 1016 m−2 and confinement strength l0 = 2.6 × 10−7 m. By further reducing l0, the droplet contracts to a small size (dashed line in figure 2).
Figure 2. Density profiles of a self-bound droplet. Parameters are : N = 104 of
Dy atoms, a = 141a0 [33] (a0 being the Bohr radius), and r* = 130a0 [23, 34]. Solid line: b = 20. Dashed line: b = 10.
Download figure:
Standard image High-resolution imageTo quantitatively check the existence of the droplet, we additionally analyze the behavior of the one-body density matrix which can be determined within the realm of the phase-density representation [35]. Writting the field operator in the form
, where
and
are the phase and density operators, which obey the commutation relation
. Expanding the density and the phase in the basis of the excitations:
and
(see e.g [27, 30, 36]). Assuming small density fluctuations, we then obtain for the excitation spectrum of homogeneous gas

where
. Equation (7) shows that the presence of the LHY quantum corrections in the dispersion relation may lead to modify the full spectrum of the system. If C varies in the narrow interval

the system emulates roton-maxon excitation spectrum. The position and the gap of the roton are shifted owing to the LHY quantum corrections (see figure 3(a)) in agreement with recent numerical and experimental predictions [26]. For
, the droplet becomes unstable and thus, the ground state completely disappears. In the limit
, the dispersion law (7) is linear in k and well approximated by the phonon-like linear dispersion form
, where the sound velocity is given by
with
being the standard sound velocity. The collective excitations of the droplet can be determined by numerically solving the full Bogoliubov–de-Gennes equations which are however, beyond the scope of the present work.
Figure 3. (a) Excitations energy
of the quasi-2D dipolar droplet as a function of momentum k. Solid line: with LHY corrections and dashed line: without LHY corrections. (b) One-body correlation function for nr*2 = 0.2 (dotted line) and nr*2 = 0.75 (solid line).
Download figure:
Standard image High-resolution imageThe one-body density matrix is defined as
. We see from figure 3(b) that the one-body density matrix is tending to its asymtotic value n at
signaling the existence of the long-range order allowing formation of a droplet in quasi-2D geometry at zero temperature. This confirms the scenario anticipated above whereby the combined effect of the quantum fluctuations and the dipolar instability may lead to a stable droplet. Close to the roton region, g1(r) is increased by ∼14% and exhibits pronounced oscillations when r is approaching to zero. These oscillations are most likely a signature of the destruction of the long-range order, unlocking the possibility of a novel quantum phase transition.
4. Effects of thermal fluctuations
In this section we deal with the finite-temperature behavior of the droplet. In uniform Bose gas, at any nonzero temperature, thermal fluctuations distroy the condensate [37, 38]. However, according to Berezinskii–Kosterlitz–Thouless (BKT) [39, 40], quasicondensate takes place at low temperature, characterized by a power-law decay of the one-body spatial correlation function [27, 30]. In such a quasicondensate, the phase coherence governs only regime of a size smaller than the size of the condensate, marked by the coherence length [30]. For a distance smaller than the coherence length, one can work with the true BEC theory [30, 41]
At finite temperature, the LHY thermal fluctuations reads

In contrast to the zero temperature case, integral (9) is finite. At low T, the main contribution to equation (9) comes from the phonon branch. This yields

where ζ (3) is the Riemann Zeta function. The most striking feature of the thermal fluctuations (10) which introduce a new extra term
, is that they change their nature from repulsive at lower T to attractive interactions at higher T. Notice that at T > μ0, the leading term for the chemical potential coincides with that of an ideal gas.
The thermal contribution to the equilibrium density can be given by minimizing the free energy
[21, 22]. This yields

its behavior as a function of l0 and T is displayed in figure 4(a). We see that the thermal equilibrium density
is important only for
and
, revealing that thermal fluctuations may substantially affact the stability of the droplet. At
and
,
and hence, the condensate is weakly depleted, results in the droplet remains in its equilibrium state.
Figure 4. (a) Equilibrium density from equation (11) as function of temperature T/E0 and confinement strength l0. (b) Condensed density for several values of temperatures T/E0. Solid line: T/E0 = 2.5. Dahsed line: T/E0 = 3.5. Dotted line: T/E0 = 8.5. Parameters are the same as in figure 2.
Download figure:
Standard image High-resolution imageLet us now look at how the condensed density inside the droplet behaves by varying the temperature. To this end, we insert the quantum (3) and thermal (10) fluctuations into the nonlocal GPE (5) and pursue typically the same numerical method. Figure 4(b) depicts that the condensed density decreases with increasing temperature. At higher T, the droplet evaporates into an expanding gas owing to the strong thermal fluctuations. Note that nc could be also shifted by changing l0/a at fixed temperature.
5. Conclusions
In conclusion, we predicted the formation of a self-bound droplet in quasi-2D dipolar Bose gas at both zero and low temperatures. Interestingly, the roton instability inducing a local collapse instability can be stabilized by the LHY corrections. Unlike the Bose mixtures [19], the LHY quantum and thermal fluctuations which present an intriguing density dependence, found to be pivotally influenced by the transversal modes. Such modes may change the nature of the LHY corrections from attractive to repulsive at certain density. Their impacts on the structure of the droplet are also considerable. The one-body correlation function of the droplet is decaying over distance and displays a remarkable behavior near the roton instability. At finite temperature, we pointed out that the droplet state can survive only at ultralow temperatures (T < E0), that should be smaller than the BKT transition temperature. Experimentally, the realization of the quantum droplet remains challenging in particular at finite temperatures due to its self-evaporation. One can expect, on the other hand, that the unusual density dependence of the quantum corrections persists also in the presence of the three-body correlations [10, 42, 43]. Future experimental investigations and Monte Carlo simulation are required in order to fully understand the confidentiality of the droplet state in 2D configuration.
Acknowledgments
We are grateful to Dmitry Petrov, Lauriane Chomaz, Krzysztof Jachymski, Grigori Astrakharchik and Pawel Zin for fruitful discussions and comments on the manuscript.
: Appendix. Low-energy S-wave scattering of dipolar bosons in quasi-2D
In this appendix we discuss low-energy two-body scattering of identical particles undergoing the 2D translational motion and interacting with each other at large separations via the potential

where r* is the characteristic dipole–dipole distance (see the main text). The term low-energy means that their momenta satisfy the inequality
.
The on-shell scattering amplitude is defined as

where Jl (kr) is the Bessel function, and
is the true wavefunction of the relative motion with momentum k. It is governed by the Schrödinger equation

where m is the reduced mass.
For the solution of the scattering problem, it is more convenient to normalize the wavefunction of the radial relative motion with orbital angular momentum l in such a way that it is real and, for
, one has

where Nl is the Neumann function and
. In order to calculate the s-wave part of the scattering amplitude, we divide the range of distances into two parts r < r0 (region I) and r > r0 (region II), where the distance r0 is selected such that
[44, 45]. Here we distinguish two contributions to the scattering amplitude : the short-range contribution coming from distances r < r*, and the so-called anomalous contribution coming from distances of the order of the de Broglie wavelength of particles, r ∼ k−1 [46].
In region I, the s-wave relative motion of two particles is governed by the Schrödinger equation with zero kinetic energy:

which admits the solution

where K0 and I0 are the modified Bessel functions and the constant A is determined by the behavior of V(r) at shorter distances. If V(r) = d2/r3 at all distances (pure dipole–dipole potential), then A = 0. In the case when V(r) behaves as d2/r3 at distances r > r1 ∼ r* (square well potential at short distance), then the coefficient A can be determined by equalizing the logarithmic derivatives of the wavefunction obtained in the region r < r1 and the one of ψI (r) at r = r1. In this model, we should have
, so that
.
In region II, the relative motion is practically free and the potential V(r) can be considered as perturbation. To zero-order we then have for the relative wavefunction

where the scattering phase shift
is due to the interaction between particles in region I.
Matching the logarithmic derivatives of ψI(r) and ψII(0)(r) at r = r0, and taking into account only terms up to r*, we get

where

with
being the Euler constant and
.
On the other hand, the contributions to the s-wave scattering phase shift from distance r > r0 should be included perturbatively. In this region, to first-order in V(r), the relative wavefunction is given by

where the Green function for the free s-wave motion is given by

Substituting the Green function (21) into (20) and taking the limit
, we have [44, 45]

where the first-order contribution to the phase shift is given by:

The sum
and the first-order contribution gives

Now for identical bosons, the full scattering amplitude
, can be obtained by making summation over all partial amplitudes with even l

We then omit the term proportional to kr*/Q2 and notice that in quasi-2D where the inequality
is satisfied, the parameter rd depends on the confinement length
in the z-direction as
[30, 45]. Substituting this into equation (24), and keeping in mind that in the quasi-2D geometry the short-range constant
, we finally obtain the result (1) employed in the main text,
.



