Many-body physics of low-density dipolar bosons in box potentials

Crystallization is a generic phenomenon in classical and quantum mechanics arising in a variety of physical systems. In this work we focus on a specific platform, ultracold dipolar bosons, which can be realized in experiments with dilute gases. We review the relevant ingredients leading to crystallization, namely the interplay of contact and dipole-dipole interactions and system density, as well as the numerical algorithm employed. We characterize the many-body phases investigating correlations and superfluidity.


Introduction
Cluster phases are ubiquitous in diverse systems, from biological structures [1] to soft matter [2,3] to ultra-cold atomic and molecular condensates [4]. In many respects, particle aggregates display fascinating properties, whose features can be mainly described microscopically by means of effective two-body potentials. In a simple mean field picture, a cluster crystalline phase (or, more formally, a spontaneous breaking of translational symmetry) appears when the Fourier transform (when it exists) of the particle-particle interaction exhibits a negative region with a local minima. The cluster phase has then a density modulation with wavelength corresponding to those minima. The validity of this condition has been investigated also in the quantum regime, and it holds also if one considers fermionic or bosonic quantum gases [5]. In particular, having potentials displaying a soft-shoulder-like shape at short distance is a crucial ingredient for designing novel quantum phases. As an example, alkaline atoms off-resonantly excited to Rydberg states display an effective interaction with a soft-core at short distances that causes a supersolid phase [6][7][8][9]. A supersolid is a phase of matter that simultaneously accommodates diagonal as well as off-diagonal long-range order, which means that particles self-assemble into a rigid, regular crystal but at the same time can they can rotate with a corresponding reduced moment of inertia, which is an indication of the presence of a finite superfluid fraction [5,[10][11][12][13][14] Recently, ground-breaking experiments with bosonic dipolar condensates (Dysprosium and Erbium) demonstrated the existence of self-bound droplets in trapped configurations as well as in free space [15,16]. These experiments provides a useful isolated system to probe general quantum-mechanical properties related to cluster phases. A mean-field treatment mainly based on a generalized nonlocal nonlinear Schrödinger equation reveal a good agreement with the density distribution and the excitation spectra observed in the laboratory [17]. The same methodology has been successfully implemented to Bose-Bose mixtures [18,19]. In addition, quantum Monte Carlo simulations (QMC) have investigated the zero temperature phase diagram of bosons interacting via dipolar interactions in three dimensions in free space [20]. Other theoretical works have shown how the physical nature of self-bound droplets smoothly evolves from classical to quantum mechanical as the range of the repulsive two-body potential increases [21,22]. In the present paper, by means of first-principles numerical simulations, we characterize the many-body physics of bosons interacting via anisotropic dipole-dipole potentials with different densities and interaction strengths. In particular we carry out simulation to better understand the quantum behaviour of the system that supports recent experimental findings about the stability of droplets.
The rest of the paper is organized as follows: In section 2 we introduce the microscopic potential that describes the physics of dipolar systems in three-dimension, whereas in section 3 we present the details of the many-body properties known in literature.In section 4, we illustrate our results, outlining our conclusions in section 5.

Two-body physics of dipole-dipole interactions
The two body-potential describing a short-range and an anisotropic long-range dipolar interaction reads where m is the particle mass, while a s is the s-wave scattering length and a d = mC dd /4πh 2 the dipolar length representing the characteristic length of the dipole interactions C dd /4π [23]. θ denotes the angle between the vector r and the polarization axis of the dipoles, the z-axis. Equation (1) applies both to vertically aligned dipoles interacting via magnetic or electric dipole moments [23]. The pseudo-potential V(r) has been already discussed by Yi and Yu [24,25] making use of the Born approximation. More realistic model potentials, that include also a Van der Walls short range potentials, were recently considered in [26] for a more detailed description of Dy two-body collisions. These potentials eventually lead to the same effective potential of [24,25] , with a dipolar length depending on the short range scattering length. In our first-principle simulations based on a QMC technique (see section 3), a d as well as a s are defined in terms of a unit length r 0 , ε 0 =h 2 /mr 2 0 being the corresponding energy with the effective pairwise potential V (r): V (r) describes a short-range hard core with cut-off at r 0 and the usual anisotropic dipolar component. The authors of [27,28] have computed numerically the full low-energy scattering amplitude for the same model potential in equation (2). A comparison with equation (1) found a quite good agreement. As a useful example, we now consider the Dysprosium's isotopes 162 Dy (a d = 129.2 a 0 , a 0 being the Bohr radius) and 164 Dy (a d = 130.8 a 0 ), relevant to the experimental observations by Kadau et al.  [15]. In addition, the background scattering lengths were recently measured for both isotopes (6-8): a s = 122 a 0 and a s = 92a 0 , for 162 Dy and 164 Dy, respectively.
In [27] the authors extracted the relation between the dipolar length a d and the scattering length a s for the model potential V (r) as shown in fig. 1. For a d /r 0 ≈ 2.9 and a d /r 0 ≈ 6.5 we observe two resonances at which the scattering length a s diverges. Moreover, one also extracts the value of d = a d /a s as a function of the corresponding ratios a s /r 0 and a d /r 0 (see Fig. 2). From the calculation of d for the two isotopes we can also compute the value of the parameter r 0 . In particular, for 162 Dy we have r 0 = 177a 0 , while considering 164 Dy we get r 0 = 154a 0 .

Many-body properties
The present section is devoted to discuss many-body properties of bosons interacting via a dipolar interaction. The quantum-mechanical Hamiltonian of an ensemble of N interacting identical bosons reads The Hamiltonian is reported in units of r 0 and ε 0 . It is important to note that the zero-temperature physics is exclusively controlled by the dimensionless interaction strength V 0 = mC dd /4πh 2 r 0 and the dimensionless density nr 3 0 . Recently the zero-temperature phase diagram of equation (3) has been throughly investigated applying QMC techniques [20]. In particular, considering low V 0 (that is 2.1 in the present units), simulations agree with the mean-field phase boundary predicted by standard Bogoliubov analysis. Upon increasing the dipole-dipole interaction we can distinguish two different phases. At low densities (nr 3 0 5 · 10 −3 ) the system displays a cluster phase with vanishing superfluidity and characterized by droplet structures with few particles. For higher densities, one observes a phase marked by elongated filaments with an anisotropic superfluid fraction; that is, whereas the superfluidity is observed along the direction of the filaments, it results greatly suppressed on the corresponding orthogonal plane, excluding then the occurrence of a global supersolid phase. In addition, finite-temperature calculations confirmed the stability of the filament phase against thermal fluctuations and provide an estimate of the superfluid fraction in the weak coupling limit in the framework of the Landau two-fluid model. It is important to stress that both cluster phase and filament phase extend from a small positive induced contact potential ( d 1), corresponding to the experimentally relevant regime of [15], to the strongly coupled limit of large dipolar interactions. Theoretical results obtained solving Hamiltonian (3) can be compared with recent experiments with 162 Dy and 164 Dy making use of the following condition: As an example, setting an average density n = 5 · 10 20 m −3 , and with na 3 d = 2 · 10 −4 , taking a d = 130a 0 [15,29], we have shown in Ref. [20] that experimental results lie in the transition region (within the error bars) from superfluid to cluster phase by tuning the scattering length.

Methodology
The low-temperature equilibrium properties of the systems described via the Hamiltonian (3) have been investigated employing first-principle computer simulations based on the worm algorithm in the continuous-space path integral representation [30], which allows one to essentially obtain the exact thermodynamics properties of Bose systems. Over the last ten years this computational technique has been successfully tested on a large variety of bosonic systems, including, for instance, 4 He [31], Rydberg atoms [4] as well as dipolar systems [32][33][34][35][36]. For the details of the implementation of the worm algorithm the reader may refer to Boninsegni et al. [30]. Due to the intrinsic anisotropy of the potential (2), we have implemented in this work the so-called primitive approximation for the imaginary time propagator, which requires a greater number of time slices with respect to more complex propagators, such as the pair-product ansatz [37] or the fourth-order one [38]. All the results reported below are extrapolated to the limit of zero temperature. In this work we study the equilibrium properties of scaled average density nr 3 0 < 0.1 and different interaction strengths V 0 . We work with N atoms in a cubic box of linear dimension L (fixing the density n = N/L 3 ) and using periodic boundary conditions. We have performed simulations with N between 100 and 400; in such a limit one can exclude finite-size issues. From these simulations we obtain density profiles, radial correlation functions, and the superfluid fraction. The limit of zero temperature is obtained lowering the temperature until observables (e.g. energy, or superfluid fraction) do not change on further decreasing T.

Density profile of the Filament Phase
In this section we present results applying the QMC method discussed in section 4.1. We consider for now a regime of high density such as nr 3 0 = 10 −2 . For this peculiar case our calculations obtained a stable filament phase at temperature much lower then the transition temperature from normal to superfluid liquid. Fig. 3a shows to a filament phase with N = 100 and P = 1000, P being the number of imaginary-time slices. This snapshot displays the projection of the single-particle imaginary-time evolution onto the real space. This representation of the condensate provides information of the delocalisation of particles in the box, and therefore of their probability distribution.
The system shapes in four different well defined filaments, that we now characterize. We fit the single-filament density with a gaussian n G (r) and extract the average number of particles within each filament. The inset in Fig. 3b shows the width of the gaussian a ⊥ in units of r 0 , as a function of the particle number N p . The estimate of the average inter-particle distance with a gaussian wave packet in the radial direction and a cylinder of length L along z yields for the data in Fig. 3a, where n f is the density of the filament. From the fit of the parameters above we get N p = 25 ± 7 and a ⊥ /r 0 = 1.7 ± 0.4.

Weakly interacting regime
With the many-body Hamiltonian (3) one can also investigate the weakly interacting regime. To do that, we consider a limit where a s and a d length scales are much smaller than the inter-particle distance. This corresponds to maintain the general conditions n|a s | 3 1 and na 3 d 1 simultaneously fulfilled. Quantitatively, as already discussed in Ref. [20], both inequalities hold keeping n|a s | 3 and na 3 d ≤ 5 · 10 −3 . We present QMC data at d = 0.6, 1.2, 1.8. The upper part of Fig. 4 reports snapshot configurations at a rescaled low temperature equal to T/T 0 = 0.25, T 0 being the critical temperature of the ideal Bose gas k B T 0 = 2π(nr 3 0 ) 2/3 /ζ(3/2) 2/3 ; whereas the bottom panels show their corresponding radial distribution functions g(r). Also in this case, simulations have been obtained with N = 100, P = 1000, keeping a fixed density of nr 3 0 = 1 · 10 −4 . For d = 0.6 ( Fig. 4a) the condensate appears modified by the dipolar strength only marginally. Radial distribution functions display strong correlations at short distances, presenting a typical fluid behavior at large r. Upon increasing the dipolar interaction, for d > 1 (e.g. Fig. 4b-c) we observe a stabilization of clusters. As expected, their emergence exhibits a g(r) with stronger correlations at short r if compared to Fig. 4a. Yet, the same observable in Fig. 4b-c shows a gentle modulation at long distances, excluding long range crystalline order at larger r, and therefore any cluster crystal phase. Also, contrarily to the prediction of mean field theory for the simulation parameters, this phase is not globally superfluid. In order to illustrate this point, we show the computed the frequency P(n) of the occurrence of permutation cycles involving n particles (1 ≤ n ≤ N) for the system with d = 1.2 (see Fig. 5). The important point here is that permutations extend only along cycles that remain confined within the single cluster, supporting the interpretation that locally particles exchange and quantum droplets are locally superfluid. Finally, increasing further d (Fig. 4c), we observe nothing but the stabilization of the cluster phase. This completely excludes the existence of the filament regime discussed at higher densities.

Conclusion
We discussed the appearance of cluster phases in ultracold dipolar bosons. We reviewed the relevant pairwise interactions among dipoles, leading to the formation of clusters and filaments, as well as the numerical algorithm employed. We characterized the many-body phases varying the strength of the potential as well as the density of the system via extensive QMC calculations.