Many-Body Physics of Low-Density Dipolar Bosons in Box Potentials

Abstract

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 reviewed 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 characterized the many-body phases investigating correlations and superfluidity.

1. 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, here we propose an original work, which aims to understand the regime of low density, where a quantum phase transition from superfluid to cluster is expected. For this reason, we carries out simulations in order to support recent experimental findings about the stability of droplet phases in free space.

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.

2. Two-Body Physics of Dipole–Dipole Interactions

The two body-potential describing a short-range and an anisotropic long-range dipolar interaction reads

$$V\left(\mathbf{r}\right)=\frac{4\pi {\hslash }^2}{m}{a}_s\delta \left(\mathbf{r}\right)+\frac{C_{dd}}{4\pi }\frac{\left(1-3{cos}^2\theta \right)}{r^3}$$

(1)

where m is the particle mass, $a_s$ is the s-wave scattering length and $a_d=m{C}_{dd}/4\pi {\hslash }^2$ is the dipolar length representing the characteristic length of the dipole interactions $C_{dd}/4\pi$ [23]. $\theta$ denotes the angle between the vector $\mathbf{r}$ and the polarization axis of the dipoles is the z-axis. Equation (1) applies to vertically aligned dipoles interacting via magnetic or electric dipole moments [23]. The pseudo-potential $V\left(\mathbf{r}\right)$ was discussed by Yi and Yu [24,25] making use of the Born approximation. More realistic model potentials, which include also a van der Walls short range potentials, are recently considered in [26] for a more detailed description of Dy two-body collisions. These potentials eventually lead to the same effective potential as in [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$ and $a_s$ were defined in terms of a unit length $r_0$, ${\epsilon }_0={\hslash }^2/m{r}_0^2$ being the corresponding energy with the effective pairwise potential $\mathcal{V}\left(\mathbf{r}\right)$:

$$\mathcal{V}\left(\mathbf{r}\right)=\left\{\begin{array}{cc}\infty \hfill & \mathrm{if}r<r_0,\hfill \\ \frac{C_{dd}}{4\pi }\frac{1-3{cos}^2\theta }{r^3}\hfill & \mathrm{if}r\ge r_0.\hfill \end{array}\right.$$

(2)

$\mathcal{V}\left(\mathbf{r}\right)$ describes a short-range hard core with cut-off at $r_0$ and the usual anisotropic dipolar component. The authors of [27,28] 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.2a_0$, $a_0$ being the Bohr radius) and ${}^{164}$Dy ($a_d=130.8a_0$), relevant to the experimental observations by Kadau et al. [15]. In addition, the background scattering lengths were recently measured for both isotopes [15]: $a_s=122a_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 $\mathcal{V}\left(\mathbf{r}\right)$, as shown in Figure 1. For $a_d/r_0\approx 2.9$ and $a_d/r_0\approx 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 Figure 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$.

3. Many-Body Properties

The present section is devoted to discussing many-body properties of bosons interacting via a dipolar interaction. The quantum-mechanical Hamiltonian of an ensemble of N interacting identical bosons reads

$$\mathcal{H}=-\frac{1}{2}\sum _{i=1}^N{\nabla }_i^2+\sum _{i<j}^N\mathcal{V}({\mathbf{r}}_i-{\mathbf{r}}_j).$$

(3)

The Hamiltonian is reported in units of $r_0$ and ${\epsilon }_0$. It is important to note that the zero-temperature physics is exclusively controlled by the dimensionless interaction strength ${\mathcal{V}}_0=m{C}_{dd}/4\pi {\hslash }^2{r}_0$ and the dimensionless density $n{r}_0^3$. Recently, the zero-temperature phase diagram of Equation (3) has been throughly investigated applying QMC techniques [20]. In particular, considering low ${\mathcal{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 ($n{r}_0^3\lesssim 5\times {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\gtrsim 1$), corresponding to the experimentally relevant regime in [15], to the strongly coupled limit of large dipolar interactions.

Theoretical results obtained solving the Hamiltonian in Equation (3) can be compared with recent experiments with ${}^{162}$Dy and ${}^{164}$Dy making use of the following condition:

$$n{r}_0^3=\frac{n{a}_d^3}{{(a_d/r_0)}^3}.$$

(4)

As an example, setting an average density $n=5\times {10}^{20}$ m${}^{-3}$, and with $n{a}_d^3=2\times {10}^{-4}$, taking $a_d=130a_0$ [15], we showed in our previous work [20] that experimental results lie in the transition region (within the error bars) from superfluid to cluster phase by tuning the scattering length.

4. Results

4.1. Methodology

The low-temperature equilibrium properties of the systems described via the Hamiltonian in Equation (3) have been investigated employing first-principle computer simulations based on the worm algorithm in the continuous-space path integral representation [29], 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 [30], Rydberg atoms [4] and dipolar systems [31,32,33,34,35]. For the details of the implementation of the worm algorithm, the reader may refer to Boninsegni et al. [29]. Due to the intrinsic anisotropy of the potential in Equation (2), we 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 [36] or the fourth-order one [37]. All the results reported below were extrapolated to the limit of zero temperature. In this work, we studied the equilibrium properties of scaled average density $n{r}_0^3<0.1$ and different interaction strengths ${\mathcal{V}}_0$. We worked with N atoms in a cubic box of linear dimension L (fixing the density $n=N/L^3$) and using periodic boundary conditions. We performed simulations with N between 100 and 400. Keeping the system’s density constant when N increases, simulations showed a corresponding growth of the number of filaments or clusters. From these simulations, we obtained density profiles, radial correlation functions, and the superfluid fraction. The limit of zero temperature was obtained lowering the temperature until observables (e.g., energy, or superfluid fraction) do not change on further decreasing T.

4.2. 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 $n{r}_0^3={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. Figure 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 delocalization of particles in the box, and therefore of their probability distribution.

The system shapes in four different well defined filaments, which we now characterize. We fit the single-filament density with a Gaussian $n_G\left(r\right)$

$$n_G\left(r\right)=\frac{N_p}{\pi a_{\perp }^{2}L}{e}^{-\frac{r^2}{a_{\perp }^2}}$$

(5)

and extract the average number of particles within each filament. The inset in Figure 3b shows the width of the Gaussian $a_{\perp }$ in units of $r_0$, as a function of the particle number $N_p$. The four different colored squares refer to the distribution for each filament of Figure 3a. 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

$$〈r〉\approx \frac{1}{n_f^{1/3}}={\left(\frac{\pi a_{\perp }^{2}L}{N_p}\right)}^{1/3}=2.1r_0$$

(6)

for the data in Figure 3a, where $n_f$ is the density of the filament. From the fit of the parameters above, we get $N_p=25\pm 7$ and $a_{\perp }/r_0=1.7\pm 0.4$. For sake of clearness, one Gaussian fit is only reported in Figure 3b (see blue line).

4.3. Weakly Interacting Regime

With the many-body Hamiltonian in Equation (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\ll 1$ and $n{a}_d^3\ll 1$ simultaneously fulfilled. Quantitatively, as already discussed in Reference [20], both inequalities hold keeping $n|a_s{|}^3$ and $n{a}_d^3\le 5\times {10}^{-3}$.

We present QMC data at ${ϵ}_d=0.6,1.2,1.8$. Figure 4 (top) 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\pi {\left(n{r}_0^3\right)}^{2/3}/\zeta {(3/2)}^{2/3}$; whereas the bottom panels show their corresponding radial distribution functions $g\left(r\right)$. In addition, in this case, simulations were obtained with $N=100$, $P=1000$, keeping a fixed density of $n{r}_0^3=1\times {10}^{-4}$.

For ${ϵ}_d=0.6$ (Figure 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., Figure 4b,c), we observed a stabilization of clusters. As expected, their emergence exhibits a $g\left(r\right)$ with stronger correlations at short r if compared to Figure 4a. However, the same observable in Figure 4b,c shows a gentle modulation at long distances, excluding long range crystalline order at larger r, and therefore any cluster crystal phase. In addition, contrary to the prediction of mean field theory for the simulation parameters, this phase is not globally superfluid. To illustrate this point, we show the computed the frequency $P\left(n\right)$ of the occurrence of permutation cycles involving n particles ($1\le n\le N$) for the system with ${ϵ}_d=1.2$ (see Figure 5). The important point here is that permutations extend only along cycles that remain confined within the single cluster or most close by clusters, supporting the interpretation that locally particles exchange and quantum droplets are locally superfluid. Finally, increasing further ${ϵ}_d$ (Figure 4c), we observed nothing but the stabilization of the cluster phase. This completely excludes the existence of the filament regime discussed at higher densities.

5. Conclusions

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.