Path-integral Monte Carlo study on a droplet of a dipolar Bose-Einstein condensate stabilized by quantum fluctuation

Motivated by the recent experiments [H. Kadau et al., Nature (London) 530, 194 (2016); I. Ferrier-Barbut et al., arXiv:1601.03318] and theoretical prediction (F. W\"achtler and L. Santos, arXiv:1601.04501), the ground state of a dysprosium Bose-Einstein condensate with strong dipole-dipole interaction is studied using the path-integral Monte Carlo method. It is shown that quantum fluctuation can stabilize the condensate against dipolar collapse.

many-body properties of ultracold atoms. [23][24][25][26][27][28][29][30][31] We will show that a stable droplet state is obtained by the PIMC method, even when the ground state does not exist and the collapse occurs in the simple mean-field theory. The density profiles of the atomic clouds obtained by the PIMC method are compared with those by the GP equation with LHY correction proposed in Ref. 20.
We consider a system of 164 Dy atoms with mass m confined in a trap potential V(r), in which the direction of the magnetic dipole moment of the atoms is fixed to the z axis. The Hamiltonian for the system is given by where N is the number of atoms, and r j and p j are the position and momentum operators of the jth atom. The system is confined in a harmonic potential V(r) = m(ω 2 x x 2 +ω 2 y y 2 +ω 2 z z 2 )/2, where ω x , ω y , and ω z are the trap frequencies. The interaction U between atoms consists of the hard-sphere potential with a radius a and the magnetic DDI as where U hard (r) = ∞ for r < a and U hard (r) = 0 for r > a, µ 0 is the magnetic permeability of the vacuum, µ = 9.93µ B is the magnetic dipole moment of a 164 Dy atom with µ B being the Bohr magneton, and χ is the angle between r and the z axis. It is known that the s-wave scattering length coincides with the radius a of the hard-sphere potential.
In thermal equilibrium at temperature T , the probability that the atoms are located at R = {r 1 , r 2 , · · · , r N } is proportional to P R|e −βH |PR , where P represents permutation of indices to assure the Bose symmetry and β = 1/(k B T ) with k B being the Boltzmann constant. The bracket R|e −βH |PR is divided into the path-integral form, where M is the number of "slices". Each bracket in Eq. (3) is approximated by The part of noninteracting particles in a harmonic potential has the form, 32 where τ = β/M. For the hard-sphere interaction part in Eq. (4), we adopt the expression derived in Ref. 33, where and H is the Heaviside step function, i.e., P hard vanishes when the distance between any two particles is less than a. The DDI part in Eq. (4) is approximated as where f cutoff (r) = r 3 for r > R cutoff and f cutoff (r) = R 3 cutoff for r < R cutoff . The cutoff radius R cutoff is introduced to avoid the steep increase in the DDI potential near r = 0, which ruins the calculation. The validity of the cutoff will be discussed later. Using these expressions of P noint , P hard , and P ddi , P R|e −βH |PR is evaluated by the multilevel Metropolis sampling, 22 where R 1 , R 2 , · · · , R M−1 , and R are sampled with an appropriate probability. Taking the average of R, one obtains the density distribution n(r) of the atomic cloud in thermal equilibrium. Typically, after 10 3 -10 4 Monte Carlo sweeps are performed to relax the system, 10 3 -10 4 samples are taken for the average.
Before showing the PIMC results, we briefly review the mean-field theory with the LHY correction proposed in Ref. 20. The LHY correction of the chemical potential in a homogeneous dipolar BEC with density n is given by 34 where g = 4π 2 a/m, ǫ dd = µ 0 µ 2 /(3g), and The integral in Eq. (9) is taken for the range in which the integrand is real. Using the local density approximation, the LHY correction in Eq. (8) is incorporated into the GP equation, giving where χ ′ is the angle between r− r ′ and the z axis. The macroscopic wave function ψ is normalized as |ψ| 2 dr = N. The DDI energy and the LHY correction ∆µ are roughly proportional to |ψ| 2 and |ψ| 3 , respectively. Therefore, when the peak density is increased by the DDI, the energy is dominated by the LHY correction term, which stops the collapse. The LHY quantum fluctuation can thus prevent the collapse and stabilize droplets. In the following results, the stationary states of the GP equation are obtained by the imaginary-time propagation method, in which i on the left-hand side is replaced with −1 and the wave function is normalized in every time step. We first check that the PIMC method reproduces the meanfield theory, when the LHY correction is small. We consider a system of N = 1024 atoms confined in a harmonic trap with frequencies (ω x , ω y , ω z ) = 2π × (46, 44, 133) Hz. 14 We take ωβ ≡ (ω x ω y ω z ) 1/3 β = 0.3, which corresponds to T ≃ 7.4 nK. The critical temperature for Bose-Einstein condensation of an ideal Bose gas is T c ≃ 0.94 ωN 1/3 /k B ≃ 29 nK. Figure 1(a) shows the results without DDI, where a = 100a 0 with a 0 being the Bohr radius. We define the integrated density distributions as where n(r) is the atom density. In Fig. 1(a), the density profiles obtained by the PIMC method almost agree with those by the GP equation. For these parameters, the GP results with and without the LHY correction cannot be discerned. Figure 1(b) shows the result with DDI for a = 70a 0 . For this value of a, the relative strength of the DDI is ǫ dd ≃ 1.87 and the GP equation without the LHY correction has a metastable state, where the energy barrier originates from the quantum pressure. 6 Due to the anisotropic nature of the DDI, the atomic cloud is slightly elongated in the z direction and shrunk in the x-y direction, compared with that without DDI. We see that the density distribution obtained by the PIMC method is in good agreement with those by the GP equation with DDI. This indicates that the PIMC method can be used to obtain not only the ground state but also a metastable state. The cutoff radius used in Fig. 1 where a x = [ /(mω x )] 1/2 . Almost the same result is obtained for R cutoff = 0.1a x . We next examine whether the quantum fluctuation can stop the dipolar collapse. The state in Fig. 1(b) is the metastable state, and beyond the energy barrier, the energy of the atomic cloud decreases as it shrinks. If the LHY correction is absent, the GP equation has no lower energy bound and the peak density diverges; there is no ground state. The LHY correction in Eq. (10) suppresses the divergence of the peak density and allows the ground state. 20 To cross the energy barrier in the numerical calculations, the radial harmonic frequencies ω x and ω y are temporarily increased during Monte Carlo sweeps in the PIMC and during imaginary-time propagation in the GP equation, which shrinks the atomic cloud in the x-y directions. Starting from these states, the system goes to the ground state beyond the energy barrier. Figure 2 shows the density distributions of the state that has crossed the energy barrier, where the parameters are the same as those in Fig. 1(b). From Fig. 2(a), we find that both PIMC method and GP equation with LHY correction provide stable states, in which the dipolar collapse is suppressed and the peak density is kept finite. The density distribution obtained by the PIMC method slightly deviates from that by the GP equation with LHY correction, mainly due to the errors in the PIMC, which will be explained later. The GP equation with LHY correction may also be inaccurate due to the local density approximation. Figure 2(b) shows the isodensity surfaces of the three-dimensional density distribution obtained by the PIMC method. The atomic cloud is highly deformed to the cigar shape by the anisotropic DDI, while the trap potential is pancake shaped. The peak density in Fig. 2 is ∼ 3 × 10 15 cm −3 and the gas parameter is na 3 ∼ 10 −4 . The three-body recombination is expected to occur predominantly at the density peak, which is the reason for the atomic loss observed in the experiment. 14 We check the validity of the cutoff made in the DDI potential in Eq. (7). Figure 3 shows the dependence of the density distribution n x (x) on the cutoff radius R cutoff and the number of slices M in the PIMC method. For R cutoff = 0.25a x , the distribution n x (x) is substantially wider than others, and we presume that R cutoff = 0.25a x is too large to give the accurate result. For R cutoff = 0.2a x , n x (x) of the PIMC is close to that of the GP equation with LHY correction. Since the results of M = 256 and M = 512 are almost the same, the number of slices is enough. However, for R cutoff = 0.15a x , n x (x) significantly depends on M, when M is inadequate; M must be 2048 or larger. Therefore, the number of slices M must be increased with a decrease in R cutoff . The computational amount is proportional to M, and the calculation for R cutoff < 0.15a x is extremely difficult. It seems that the density distribution converges to around that of the GP equation with LHY correction.
The accuracy of the present PIMC calculation is thus re-  stricted by the the primitive approximation of P ddi in Eq. (7), whose r −3 steepness requires large numbers of slices M and Monte Carlo samplings. A more suitable expression for P ddi is needed. Another bottleneck is the long-range nature of the DDI, which costs O(N 2 ) calculations per Monte Carlo sweep. The O(N) Monte Carlo technique 35 may circumvent this problem. With these improvements, it may be possible not only to perform more accurate calculation, but also to simulate the droplet pattern formation observed in the experiment with N ∼ 10 5 atoms. 14 The results obtained by the PIMC method should be compared with other methods, such as the diffusion Monte Carlo method.
In conclusion, we have investigated the stability of a strong dipolar BEC against collapse, motivated by the recent experiments 14,15 and the theoretical prediction. 20 Using the PIMC method, we showed that the system has a stable ground state even in the parameters for which the simple GP equation cannot sustain the system against dipolar collapse, which implies that the quantum fluctuation stabilizes the system. We compared the PIMC results with those obtained by the GP equation with the LHY correction proposed in Ref. 20, and found that they are in qualitative agreement. The present results indicate that the quantum fluctuation plays an important role in the droplet stabilization observed in the experiments. 14,15