Droplets of trapped quantum dipolar bosons
A. Macia, J. S´anchez-Baena, J. Boronat, and F. Mazzanti
Departament de F´ısica, Campus Nord B4-B5, Universitat Polit`ecnica de Catalunya, E-08034 Barcelona, Spain
Strongly interacting systems of dipolar bosons in three dimensions confined by harmonic traps are
analyzed using the exact Path Integral Ground State Monte Carlo method. By adding a repulsive
two-body potential, we find a narrow window of interaction parameters leading to stable ground-
state configurations of droplets in a crystalline arrangement. We find that this effect is entirely due
to the interaction present in the Hamiltonian without resorting to additional stabilizing mechanisms
or specific three-body forces. We analyze the number of droplets formed in terms of the Hamiltonian
parameters, relate them to the corresponding s-wave scattering length, and discuss a simple scaling
model for the density profiles. Our results are in qualitative agreement with recent experiments
showing a quantum Rosensweig instability in trapped Dy atoms.
Dipolar effects in quantum gases have been considered
of major experimental and theoretical interest in the last
decade since the initial studies of dilute clouds of Cr
atoms, which present a relatively large magnetic dipo-
lar moment. In the pioneering experiments of Ref. [1],
the two-body scattering length of a cloud of
52
Cr atoms
was drastically reduced by bringing it close to a Fesh-
bach resonance. In this way, dipolar effects were en-
hanced and interesting new features, not observed be-
fore in other species like Rb or Cs, appeared. The long
range and anisotropic character of the dipolar interac-
tion has been largely explored since then, leading to in-
teresting new phenomena such as d-wave superfluidity
or d-wave collapse [2, 3]. All these experiments have
opened new perspectives on the field of dipolar quantum
physics, and new systems with stronger dipolar interac-
tions have since then been explored. The most promis-
ing ones, consisting initially in ultracold polar molecules
of K and Rb or Cs and Rb, are unfortunately problem-
atic due to the inherent difficulty to bring them down to
the quantum degeneracy limit, although recent progress
have been achieved with NaK [4] and NaRb [5] molecules.
Alternatively, Bose-Einstein condensate (BEC) states of
Er [6] and Dy [7] have recently been produced, enabling
for instance to observe the deformation of the Fermi sur-
face [8], or the influence of the anisotropy on the super-
fluid to Mott insulator transition in the extended Hub-
bard model made with dipoles [9].
The anisotropy of the dipolar interaction plays a fun-
damental role on the behavior of the system, with differ-
ent regimes and phases depending on the geometry and
dimensionality. The particular form of the dipole-dipole
potential makes the interaction attractive or repulsive
depending on the relative orientation of the dipoles, ac-
cording to the expression
V
dd
(r) =
C
dd
4π
ˆp
1
·ˆp
2
− 3(ˆp
1
·ˆr)(ˆp
2
·ˆr)
r
3
, (1)
where C
dd
sets the strength of the interaction that is
proportional to the square of the (magnetic or electric)
dipolar moment, p
j
is the dipolar moment itself, and
r is the relative position vector of the two interacting
dipoles. The particular form of this interaction leads to
surprising new features not present in other systems, like
stripe phases in two-dimensional (2D) Bose systems [10].
Similar phases in Fermi systems have also been pre-
dicted [11, 12], although these are more controversial [13].
One of the most interesting phenomenon recently re-
ported in the field of dipolar quantum gases is the for-
mation of self-bound droplets when a gas of trapped
164
Dy atoms is brought to the regime of mean field col-
lapse [14, 15]. More interestingly, the resulting set of
droplets are reported to arrange themselves in a crys-
talline structure that becomes more clearly visible when
the number of droplets increases. This phenomenon re-
sembles the Rosensweig instability of a classical ferrofluid
appearing when the magnetization increases, and thus
the dipolar gas can be considered as a first realization
of a quantum ferrofluid. From the theoretical side, this
is a clear example of a situation where beyond mean-
field effects drastically determine the physics of the sys-
tem. Initial mechanisms based on the inclusion of three-
body forces were reported to produce both effects [16–
18] (droplet formation and crystallization), although this
mechanism does not seem to be fully compatible with
experimental data [15]. Beyond mean field effects at the
level of the Lee-Huang-Yang correction (LHY) [19], as
proposed in [20–22], produce equivalent effects. A sim-
ilar stabilization mechanism has been recently proposed
to make possible the formation of liquid droplets in at-
tractive Bose-Bose mixtures [23].
In this work, we address the problem of droplet for-
mation of trapped dipolar bosons from a microscopic
point of view, using the stochastic Path Integral Ground
State (PIGS) method [24–26]. Starting from a variational
ansatz for the ground-state wave function, propagation
in imaginary time leads to a statistically exact represen-
tation of the actual ground state of the system that is
used to sample relevant observables. Differently from the
perturbative approaches discussed before, PIGS includes
correlations induced by the interactions to all levels. Ac-
curate fourth-order approximations of the many-body
propagator are employed in order to guarantee conver-
gence with the lowest possible number of imaginary time
2
steps. In particular, we have employed the 4A fourth-
order short time expansion propagator proposed by Chin
and Chen in Ref. [27], and have set the trial wave func-
tion to a constant. Recently, Saito [28] performed a path
integral Monte Carlo simulation of a single droplet but
the formation of the ordered array of droplets found in
experiments was not reproduced.
In the following we describe a system of N trapped
dipolar bosons of mass m, fully polarized along the z-
axis, by the Hamiltonian
H = −
¯h
2
2m
X
j
∇
2
j
+
X
j
V
ho
(r
j
)+
X
i<j
V
σ
(r
ij
)+
X
i<j
V
dd
(r
ij
) ,
(2)
where
V
ho
(r) =
1
2
m
ω
2
x
x
2
+ ω
2
y
y
2
+ ω
2
z
z
2
(3)
is the three-dimensional harmonic trapping potential of
frequencies ω
x
, ω
y
and ω
z
, and V
σ
(r) is a short-range two-
body repulsive interaction of the form V
σ
(r) = (σ/r)
12
,
where the parameter σ can be varied to tune the s-
wave scattering length. This interaction is so steep at
short distances that in practical terms it is hardly dis-
tinguishable from a pure hard core. In the following we
use dimensionless quantities, introducing a characteris-
tic length scale r
0
= mC
dd
/(4π¯h
2
) and a characteristic
energy scale E
0
= ¯h
2
/(mr
2
0
).
Inspired by recent experimental measurements [14], we
have considered different pancake-shaped harmonic traps
with oscillator lengths a
x
= a
y
> a
z
(a
α
=
p
¯h/(mω
α
)).
Two different choices, leading to more than one stable
droplet, have been used: Trap1 with a
x
= 1.20, a
z
= 0.6,
and Trap2 with a
x
= 1.38, a
z
= 0.45. Additionally, a
Trap3 model yielding a single stable droplet with a
x
=
1.73 and a
z
= 1.00 has also been considered. Notice
that the first choice is close to the experimental setting
of Ref. [14] where a
x
∼ a
z
√
3, while in the second model
the aspect ratio is different. With these settings, different
values of σ lead to different ground-state configurations.
In general, the resulting ground state is a gas, but we
have checked that for the number of particles in the few-
hundreds range used, there is a narrow window of values,
spanning the range σ ∈ [0.24, 0.28], where stable config-
urations of droplets are obtained.
Results of the ground state configurations obtained
from the PIGS simulations are reported in Fig. 1. Plots
(a), (b) and (c) correspond to Trap1, and (d), (e) and (f)
to Trap2. Each plot corresponds to σ = 0.28 but different
number of particles, N = 120, 150, 270, 90, 120, 150 (from
(a) to (f)), leading to an increasing number of droplets
as can be seen in the figure. These representative config-
urations show well-formed droplets that are easily distin-
guishable from each other, although not every configura-
tion after thermalization is that neat. This indicates that
at some point the system may reach metastable configu-
rations. Notice that, despite the number of particles in
FIG. 1: Spatial density distributions of the stable droplet
configurations obtained from the quantum simulations. Axes
in (a), (b) and (c) range from −4r
0
to 4r
0
while those in (d),
(e) and (f) span the range from −5r
0
to 5r
0
.
all cases is similar, Trap1 tends to form fewer but larger
droplets than Trap2. This fact is more clearly seen in
Fig. 2, where the average number of droplets as a function
of the number of particles for the Trap1 and Trap2 mod-
els is reported. The error bars indicate the fluctuations
that results from the averaging procedure performed over
the different configurations used in each case. The am-
plitude of the error bars can then be taken as a measure
of the degree of metastability of the analyzed states.
As it can be seen, in the Trap2 case the observed de-
pendence is linear, in agreement with the behavior ob-
served in the experiments of Ref. [14]. However, in the
Trap1 case an approximately linear dependence at the
beginning of the curve is followed by a change of ten-
dency, where one can not decide if the average number of
droplets saturates or increases very slowly with the total
number of particles in the simulation. A similar behavior
has been reported in Ref. [20] using the LHY correction,
while in that case the number of atoms in the system is
much larger.
The stochastic sampling of the different configurations
allows for the statistical analysis of the density profiles
of the droplets formed in each case. In general, chang-
3
FIG. 2: Number of droplets as a function of the total number
of particles in the simulation. Upper and lower plots corre-
spond to the Trap1 and Trap2 set of parameters explained in
the text, respectively.
ing the model parameters (σ, number of particles and
trapping frequencies) leads to different droplet configu-
rations and density profiles. We analyze the specific cases
where a reduced number of droplets appear, due to the
restrictions on the total number of particles employed.
We report results for the normalized and marginalized
n(x), n(y) and n(z) density profiles. The normalization
has been set such that the integral of each separate pro-
file is one, i.e.,
R
dα n(α) = 1 for α = {x, y, z}. Notice
that due to axial symmetry, the n(x) and n(y) profiles
are identical up to statistical errors.
The obtained density profiles are shown in Fig. 3. The
upper plot shows n(x), n(y) and n(z) for the Trap3 model
with N = 150 particles. All three profiles are approxi-
mately Gaussian but with marked differences in width
when the z component is compared to the other two.
This is due to the combined effect of the anisotropy of
the dipolar interaction and the shape of the confining
trap. Despite the fact that the trap is tighter in the z di-
rection, the droplet is larger along the z axis. In fact, the
dipolar interaction is attractive along this line, and thus
particles prefer to arrange themselves in head-to-tail con-
figurations with a minimum distance imposed by the core
of V
σ
(r). Therefore, the height of the droplet is mainly
constrained by ω
z
. However, the width r
x
' r
y
along the
radial directions is set by a more complex combination
of parameters, as can be seen from the fact that a single
droplet does not cover the whole available space. Actu-
ally, according to the obtained profiles, the self-induced
FIG. 3: Upper panel: normalized marginal density profiles
n(x), n(y) and n(z) of the single droplet obtained in the Trap3
setup with N = 150. Middle and lower panels: scaled n(z)
and n(x) density profiles for the droplets found in the Trap1
and Trap3 setups, with 120 and 150 particles per droplet,
respectively.
confinement along the radial direction is much stronger,
forming prolate stable droplets.
One simple model that captures these features assumes
that the average density of the droplet does not change
much in the range of parameters analyzed in this work.
Assuming also that the height of the droplet is essen-
tially fixed by ω
z
, a simple 1D harmonic oscillator model
leads to a height directly proportional to the a
z
oscilla-
tor length. These two assumptions also imply a simple
dependence of r
x
(∼ r
y
) on
p
N
d
/a
z
, with N
d
the total
number of particles in the droplet. In order to check the
accuracy of this simple model, the middle and lower pan-
els in Fig. 3 depict scaled z and x profiles for the Trap1
and Trap3 cases and a total number of particles N = 120
and N = 150, respectively. As it can be seen, this simple
model seems to capture the main trends reasonably well,
although the fine detail due to more complex contribu-
tions is missing.
4
FIG. 4: Upper left panel: spatial density distribution of the
averaged stable Trap2 configuration with nine droplets. Up-
per right panel: density distribution in momentum space of
the configuration shown in the upper left panel. Lower panel:
angularly averaged density in momentum space.
Following the analysis of Ref. [14], we report in Fig. 4
the Fourier spectrum of the averaged stable Trap2 config-
uration with nine droplets shown in the upper left panel.
The resulting Fourier transform density is depicted in the
upper right panel. Notice that the latter shows clear reg-
ularities that are due to the crystalline spatial distribu-
tion of the droplet configuration. The lower panel shows
the radially averaged density ρ
q
k
2
x
+ k
2
y
in momen-
tum space. This quantity presents a large value around
the origin proportional to the average number of parti-
cles and, more importantly, a clear peak around k ∼ 3
which reflects the spatial periodicity of the underlying
triangular lattice formed by the droplets.
A better comparison with experiments requires the
analysis of the range of scattering lengths covered in our
simulations. As previously mentioned, the interaction
parameters where stable configurations of droplets are
found span the range from σ = 0.24 to 0.28 in dipo-
lar units. Relating these to the s-wave scattering length
a
0
is not immediately trivial, as it is known that the
presence of the dipolar interaction in the trap modifies
its value [29, 30]. Following a standard procedure, we
have determined the low-energy behavior of the coupled-
channel T-matrix associated to the two-body potential
employed, including the full dipole-dipole interaction.
Table I reports the values of a
0
obtained for the differ-
ent σ’s used in this work. In contrast to the experi-
σ a
0
hN
droplet
i ∆hN
droplet
i
0.24 -1.417 1.60 0.44
0.25 -0.868 1.85 0.44
0.26 -0.618 2.79 0.31
0.27 -0.472 2.85 0.69
0.28 -0.374 3.25 0.19
TABLE I: s-wave scattering length, average number of
droplets hN
drop
i and standard deviation ∆hN
drop
i for the
Trap1 model with N = 120 particles, as a function of σ.
ments [14, 15], in our case all situations leading to stable
droplet configurations correspond to negative values of
a
0
, appearing before a first resonance is found. This is
due to the attractive components of the dipolar interac-
tion. In any case, the range of σ’s (and therefore scat-
tering lengths) leading to droplet formation may change
when a much larger number of particles and/or different
trapping frequencies are used, as in the experiments. The
dependence of the number of droplets on a
0
is also indi-
cated in the third and fourth columns of Table I for the
Trap1 set of parameters. As it can be seen, increasing a
0
leads to a larger number of smaller droplets.
In summary, we have used the Path Integral Ground
State algorithm with a fourth-order propagator to deter-
mine the ground state structure of a system of trapped
dipolar bosons with an additional repulsive (σ/r)
12
core.
Contrarily to previous perturbative estimates, PIGS is
exact and relies only on the Hamiltonian of the system
and the chosen geometry. The formation of self-bound
droplets and its arrangement in a crystal lattice is on
the basis of the microscopic Hamiltonian and, therefore,
there is no need to resort to three-body interactions We
find that, for a few hundred dipoles, there is a window of
σ’s corresponding to negative s-wave scattering lengths
where the ground state forms a crystal of stable droplets,
rather than collapsing to a single one or remaining as a
gas. We have analyzed the density profiles of the droplets
to find that a very simple model is able to qualitatively
explain how these scale when the number of particles
and trapping parameters are changed. The droplets are
prolate in the z direction with a height which increases
with a
z
, resembling the formation of filaments or stripes
already predicted to exist in quantum tilted dipoles in
two dimensions [31]. Finally, the analysis of the density
profiles in momentum space supports the general picture
of crystalline order in the droplet configurations. Work
is in progress now to determine the superfluidity of the
system and its dependence on the temperature using the
path integral Monte Carlo method. A supersolid scenario
seems plausible at very low temperature.
We acknowledge partial financial support from the
MICINN (Spain) Grant No. FIS2014-56257-C2-1-P.
5
[1] A. Griesmaier, J. Werner, S. Hensler, J. Stuhler, and T.
Pfau, Phys. Rev. Lett.94, 160401 (2005).
[2] T. Lahaye, T. Koch, B. Fr¨ohlich, M. Fattori, J. Metz, A.
Griesmaier, S. Giovanazzi, and T. Pfau, Nature 448, 672
(2007).
[3] T. Lahaye, J. Metz, B. Fr¨ohlich, T. Koch, M. Meister,
A. Griesmaier, T. Pfau, H. Saito, Y. Kawaguchi and M.
Ueda, Phys. Rev. Lett.101, 080401 (2008).
[4] J. W. Park, S. A. Will, and M. Zwierlein, Phys. Rev.
Lett. 114, 205302 (2015).
[5] M. Guo, B. Zhu, B. Lu, X. Ye, F. Wang, R. Vexiau, N.
Bouloufa-Maafa, G. Qu´em´ener, O. Dulieu, and D. Wang,
Phys. Rev. Lett. 116, 205303 (2016).
[6] K. Aikawa, A. Frisch, M. Mark, S. Baier, A. Rietzler, R.
Grimm, and F. Ferlaino, Phys. Rev. Lett. 108, 210401
(2012).
[7] M. Lu, N. Q. Burdick, S. H. Youn, and B. L. Lev, Phys.
Rev. Lett 107, 190401 (2011).
[8] K. Aikawa, S. Baier, A. Frisch, M. Mark, C. Ravensber-
gen and F. Ferlaino, Science 345, 1484 (2014).
[9] S. Baier, M. J. Mark1, D. Petter, K. Aikawa1, L.
Chomaz1, Z. Cai, M. Baranov, P. Zoller and F. Ferlaino,
Science 352, 201 (2016).
[10] A. Macia, D. Hufnagel, F. Mazzanti, J. Boronat and R.
E. Zillich, Phys. Rev. Lett. 109, 235307 (2012).
[11] M. M. Parish and F, M. Marchetti, Phys. Rev. Lett. 108,
145304 (2012).
[12] Y. Yamaguchi, T. Sogo, T. Ito and T. Miyakawa, Phys.
Rev. A 82, 013643 (2010).
[13] N. Matveeva and S. Giorgini, Phys. Rev. Lett. 109,
200401 (2012).
[14] H. Kadau, M. Schmitt, M. Wenzel, C. Wink, T. Maier,
I. Ferrier-Barbut and T. Pfau, Nature 530, 194 (2016).
[15] I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel and
T. Pfau, Phys. Rev. Lett. 116, 215301 (2016).
[16] R. N. Bisset and P. B. Blakie, Phys. Rev. A 92, 061603
(2015).
[17] P. B. Blakie, Phys. Rev. A 93, 033644 (2016).
[18] K.-T. Xi and H. Saito, Phys. Rev. A 93, 011604(R)
(2016).
[19] T. D. Lee, K. Huang and C. N. Yang, Phys. Rev. 106,
1135 (1957).
[20] F. W¨achtler and L. Santos, Phys. Rev. A 93, 061603(R)
(2016).
[21] F. W¨achtler and L. Santos, arXiv:1605.08676 (2016).
[22] R. N. Bisset, R. M. Wilson, D. Baillie, and P. B. Blakie,
arXiv:1605.04964 (2016).
[23] D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
[24] A. Sarsa, K. E. Schmidt and W. R. Magro, J. Chem.
Phys. 113, 1366 (2000).
[25] D. E. Galli and L. Reatto, Mol. Phys. 101, 1697 (2003).
[26] R. Rota, J. Casulleras, F. Mazzanti, and J. Boronat,
Phys. Rev. E 81, 016707 (2010).
[27] S. A. Chin and C. R. Chen, J. Chem. Phys. 117, 1409
(2002).
[28] H. Saito, J. Phys. Soc. Jap. 85, 053001 (2016).
[29] D. C. E. Bortolotti, S. Ronen, J. L. Bohn and D. Blume,
Phys. Rev. Lett. 97, 160402 (2006)
[30] S. Ronen, D. C. E. Bortolotti, D. Blume and J. L. Bohn,
Phys. Rev. A 74, 033611 (2006)
[31] A. Macia, J. Boronat, and F. Mazzanti, Phys. Rev. A 90,
061601(R) (2014).
