Droplets of trapped quantum dipolar bosons

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 dipolar 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 Feshbach resonance. In this way, dipolar effects were enhanced and interesting new features, not observed before in other species like Rb or Cs, appeared. The long range and anisotropic character of the dipolar interaction has been largely explored since then, leading to interesting 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 interactions have since then been explored. The most promising ones, consisting initially in ultracold polar molecules of K and Rb or Cs and Rb, are unfortunately problematic 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 surface [8], or the influence of the anisotropy on the superfluid to Mott insulator transition in the extended Hubbard model made with dipoles [9].
The anisotropy of the dipolar interaction plays a fundamental role on the behavior of the system, with different regimes and phases depending on the geometry and dimensionality. The particular form of the dipole-dipole potential makes the interaction be attractive or repulsive depending on the relative orientation of the dipoles, according to the expression 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 predicted [11,12], although these are more controversial [13].
One of the most interesting phenomenon recently reported in the field of dipolar quantum gases is the formation of self-bound droplets when a gas of trapped 164 Dy atoms is brought to the regime of mean field collapse [14,15]. More interestingly, the resulting set of droplets are reported to arrange themselves in a crystalline structure that becomes more clearly visible when the number of droplets increases. This phenomenon resembles 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 meanfield effects drastically determine the physics of the system. Initial mechanisms based on the inclusion of threebody forces were reported to produce both effects [16][17][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][21][22], produce equivalent effects. A similar stabilization mechanism has been recently proposed to make possible the formation of liquid droplets in attractive Bose-Bose mixtures [23].
In this work, we address the problem of droplet formation of trapped dipolar bosons from a microscopic point of view, using the stochastic Path Integral Ground State (PIGS) method [24][25][26]. Starting from a variational ansatz for the ground-state wave function, propagation in imaginary time leads to a statistically exact representation 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. Accurate fourth-order approximations of the many-body propagator are employed in order to guarantee convergence with the lowest possible number of imaginary time arXiv:1607.07184v1 [cond-mat.quant-gas] 25 Jul 2016 steps. In particular, we have employed the 4A fourthorder short time expansion propagator proposed by Chin and Chen in Ref. [27]. 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 by the Hamiltonian is the three-dimensional harmonic trapping potential of frequencies ω x , ω y and ω z , and V σ (r) is a short-range twobody repulsive interaction of the form V σ (r) = (σ/r) 12 , where the parameter σ can be varied to tune the swave scattering length. This interaction is so steep at short distances that in practical terms it is hardly distinguishable from a pure hard core. In the following we use dimensionless quantities, introducing a characteristic 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 α = 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 fewhundreds range used, there is a narrow window of values, spanning the range σ ∈ [0.24, 0.28], where stable configurations 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 configurations show well-formed droplets that are easily distinguishable from each other, although not every configuration after thermalization is that neat. This indicates that at some point the system may reach metastable configurations. Notice that, despite the number of particles in 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 models is reported. The error bars indicate the fluctuations that results from the averaging procedure performed over the different configurations used in each case. The amplitude 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 dependence is linear, in agreement with the behavior observed 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 tendency, 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, changing the model parameters (σ, number of particles and trapping frequencies) leads to different droplet configurations 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 profile is one, i.e., 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 approximately 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 direction, 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 configurations 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. Actually, according to the obtained profiles, the self-induced 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 essentially fixed by ω z , a simple 1D harmonic oscillator model leads to a height directly proportional to the a z oscillator length. These two assumptions also imply a simple dependence of r x (∼ r y ) on 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 panels 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 contributions is missing.
Following the analysis of Ref. [14], we report in the Fourier spectrum of the averaged stable Trap2 configuration 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 regularities that are due to the crystalline spatial distribution of the droplet configuration. The lower panel shows the radially averaged density ρ k 2 x + k 2 y in momentum space. This quantity presents a large value around the origin proportional to the average number of particles 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 dipolar 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 coupledchannel 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 different σ's used in this work. In contrast to the experiments [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 interaction. In any case, the range of σ's (and therefore scattering 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 indicated 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 determine 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.