Off-diagonal long-range order in arrays of dipolar droplets

We report quantum Monte Carlo results of harmonically confined quantum Bose dipoles within a range of interactions covering the evolution from a gas phase to the formation of an array of droplets. Scaling the experimental setup to a computationally accessible domain we characterize that evolution in qualitative agreement with experiments. Our microscopic approach generates ground-state results free from approximations, albeit with some controlled statistical noise. The simultaneous estimation of the static structure factor and the one-body density matrix allows for a better knowledge of the quantum coherence between droplets. Our results show a narrow window of interaction strengths where diagonal and off-diagonal long-range order can coexist. This domain, which is the key signal of a supersolid state, is reduced with respect to the one predicted by the extended Gross–Pitaevskii equation. Differences are probably due to an increase of attraction in our model, observed previously in the calculation of critical atom numbers for single dipolar drops.


I. INTRODUCTION
The experimental realization of cold Bose and Fermi gases with dipolar interactions opened a very rich landscape to explore new phenomena [1][2][3][4][5].What makes these systems exceptional and motivates the interest, from theory and experiment, is the interplay between a long range interaction and their anisotropic character [6].The first experimental realization of a dipolar quantum gas was carried out with 52 Cr [1] atoms, which have a magnetic dipolar moment µ = 6 µ B with µ B the Bohr magneton.By bringing the system close to a Feshbach resonance, it was possible to significantly reduce the scattering length, making it much smaller than the dipolar length and thus enhancing its dipolar character.More recently, the achievement of Bose-Einstein condensate (BEC) states of lanthanides Er [5] and Dy [3,4] atoms, with larger magnetic moments, µ = 7 µ B and 10 µ B , respectively, has given access to a regime where the dipolar interaction dominates.This is more evident by comparing the dipolar length of the three elements a dd = 15 a 0 , 65.5 a 0 , and 131 a 0 (a 0 is the Bohr radius) for Cr, Er, and Dy, respectively.
The dipolar interaction between head-to-tail oriented dipoles is attractive while it becomes repulsive in parallel configurations.In the absence of any other interaction to compensate the attraction, mean-field theories predict a collapsed state [7].The scenario is similar to the one of attractive Bose-Bose mixtures and, in both cases, it has been found that the introduction of quantum fluctuations through the perturbative Lee-Huang-Yang (LHY) correction is able to stabilize the system [8][9][10].The result is the formation of self-bound liquid droplets of extremely low densities, much lower than other liquids in Nature [11].At difference with spherical Bose-Bose drops, which can be produced without confinement, dipolar ones are elongated and require of the presence of an external potential.
Squeezing the harmonic trap in the direction of the oriented magnetic moments, and varying the scattering length a, different phases have been observed [12].Increasing the ratio ϵ dd = a dd /a, the system evolves from a gas phase to the formation of a single drop, followed by the appearance of an array of aligned drops.Interestingly, this array is initially coherent and then becomes insulated, just as in a crystal of droplets.Within the small window where local ordering and coherence are simultaneously observed, this phase constitutes a supersolid of droplets [13][14][15].This ordered superfluid phase arises from the coexistence of the droplets with a small density bath that connects the droplets.In this sense, this is significantly different from the supersolid phase that has been searched for many years in hcp solid 4 He [16].
The most common theoretical tool used in the study of dipolar droplets, and its eventual supersolid character, is the extended Gross-Pitaevskii equation (eGPE), including LHY beyond mean-field effects [9,10,17,18].The range of interaction parameters where supersolidity is observed is quite narrow [6], and has been estimated in two different ways.In the first case, a monopole mode is induced and the dynamic response of the system is measured [14].Within the supersolid regime two modes appear, the lowest one being associated to the coherence between drops.This low-frequency mode is the signature of superfluidity and disappears above a characteristic value of the scattering length [14].In the second approach, the superfluid density is obtained from the response of the system to a translational movement [19].
In the present work, we report results for a harmonically-confined dipolar gas of 162 Dy atoms, obtained from an ab initio approach that relies exclusively on the microscopic Hamiltonian of the system.Our results are valid in the zero temperature limit, as they have been obtained using the path integral ground-state (PIGS) method [20,21], a ground-state version of the path integral Monte Carlo (PIMC) method devised for quantum simulations at finite temperatures.Since PIGS solves exactly the many-body boson problem, up to some statistical noise, one has direct access to two-body distri-arXiv:2310.10222v2[cond-mat.quant-gas]6 Feb 2024 butions, which are not accessible using the eGPE due to its mean-field character.By tuning the s-wave scattering length a, which we obtain from the zero-momentum limit of the two-body scattering T-matrix [22], we study the evolution of the harmonically trapped dipolar gas.We observe the formation of self-bound droplets and calculate the static structure factor and the one-body density matrix as a function of a.The simultaneous observation of Bragg peaks in the former and off-diagonal long-rangeorder in the latter proves the presence of a supersolid state.Our results provide a narrow window of interaction strengths where it is plausible to conclude that this exotic state of matter is observed.
The paper is organized as follows.The PIGS method and the main ingredients of our study are discussed in Sec.II.Results are presented and discussed in Sec.III, paying particular attention to the static structure factors and one-body density matrix.Finally, the main conclusions are reported in Sec.IV.

II. QUANTUM MONTE CARLO METHOD
The starting point in any ab initio microscopic approach to a quantum many-body problem is the Hamiltonian.In the present case, it is given by with {r i } the particle position coordinates, and the dipolar interaction.In Eq. ( 2), C dd is proportional to the square of the dipolar moments, and (r ij , θ ij ) stands for the polar coordinates representation of r ij .The term V s (r ij ) = V s (r ij ) is a central potential that is repulsive at short distances to prevent collapse.It has been chosen to be of the 6-12 Lennard-Jones form, as in Ref. [22], The coefficient C 6 in Eq. ( 3) is known for Dysprosium (C 6 = 2.978 10 −2 in dipolar units) [23].The second parameter σ 12 is fixed by imposing that the full interaction (V s + V d ) has the desired s-wave scattering length.This is accomplished by solving the low momentum limit of the scattering T-matrix [22].The values of σ 12 decrease monotonously with a from 7.852 10 −4 at a = 71.1 to 2.954 10 −4 at a = 45.Finally, V h (r) is a confining harmonic potential.
In our study, we use the first-principles PIGS algorithm which is an exact stochastic projection method used to solve the Bose N -body quantum problem.Basically, it consists in a systematic improvement of a trial wave function Ψ T by repeated application of the propagator in imaginary time, which drives the system to the ground state [20], While the exact G(R ′ , R; τ ) is in general not accessible, short-time approximations are available and allow for the evaluation of averages of diagonal observables, mapping the quantum many-body system onto a classical system of N interacting polymers composed by 2M +1 beads.Each one of these terms represents a time slice of the evolution from the initial trial state Ψ T to a final wave function Ψ PIGS that is sampled at the central bead.By increasing M one is able to reduce the systematic error and therefore to exactly recover the ground-state properties of the system up to residual statistical noise.
A good approximation for the propagator G(R ′ , R; τ ) is capital to improve the numerical efficiency of the method.This greatly reduces the complexity of the calculation but also helps sorting possible ergodicity related issues, allowing to simulate the quantum system with few beads and a larger time step.By using a high-order approximation for the propagator, as the one used in this work, it is possible to obtain an accurate description of the exact ground-state wave function, even when the initial trial state contains no more information than the bosonic statistics, already present in Ψ T = 1 [21].The PIGS method has been used previously in the determination of the critical atom number required for droplet formation [22] and in the study of two-dimensional(2D) dipolar gases [24,25].Its extension to finite temperature, the path integral Monte Carlo method, has also been used to study the BKT phase transition in 2D dipolar gases [26].
The presence of diagonal order in the system is studied from the analysis of the static structure factor S(k), defined as with ρ(k) = N i=1 exp(ik • r i ) the density fluctuation operator.In PIGS, S(k) and other diagonal magnitudes, such as the density profiles ρ(r), are obtained using the coordinates of the middle points of the chains since it is in these points where the exact ground-state wave function is sampled, provided a sufficiently large number of beads is used.
One of the main goals of our work is to estimate the possible off-diagonal long-range order (ODLRO) present in the system, particularly when the formation of droplets is observed.To this end, we calculate the one-body density matrix (OBDM), where the configuration R differs from R ′ only by the position of one atom.In homogeneous systems, translational invariance makes the OBDM a function of the coordinate differences only, so ρ 1 (r 1 , r ′ 1 ) = ρ 1 (r 1 − r ′ 1 ) = ρ 1 (r).If the system is anisotropic, as in the current case, ρ 1 can be expanded in partial waves with (r, θ, ϕ) the polar coordinates of vector r and Y l,m the Spherical Harmonics.The lowest order mode, corresponding to l = m = 0, provides the only isotropic contribution and is the one being analyzed in this work.
In the following, we simplify the notation and denote ρ 0,0 1 (r) as ρ 1 (r).In extended systems, one concludes that ODLRO exists if a plateau is observed in ρ 0,0 1 (r → ∞), which coincides with the condensate fraction value n 0 .In PIGS, the expectation value of non-diagonal observables, like ρ 1 , is computed mapping the quantum system into the same classical system of polymers as in the diagonal case, but cutting one of these polymers in the mid point.Building the histogram of distance frequencies between the cut extremities of the two half polymers, one can compute the isotropic contribution ρ 1 (r).This is effectively carried out using the worm algorithm (WA), a technique previously developed for path integral Monte Carlo simulations at finite temperature [27].The key aspect of the WA is to work in an extended configuration space, containing both diagonal (all polymers with the same length) and off-diagonal (one polymer cut in two separate halves) configurations.

III. RESULTS
We study a system of 162 Dy atoms with a magnetic dipolar moment 10 µ B .All our simulations contain N = 400 harmonically confined particles.We optimize the number of time steps (beads) and the time-step itself, resulting in the values M = 30 and τ = 0.3.We use dipolar units everywhere, with r 0 = mC dd 4πℏ 2 = 3a dd , except when stated otherwise.The experimental setup of Ref. [28] uses N = 40000 particles within an harmonic potential of oscillator lengths l x = 68.908r 0 and l y = l z = 37.744 r0.Approximating the trap by a cylinder the mean experimental density is ρ ≃ 0.06485 r −3 0 .In our simulations, the size of the trap is adjusted to reproduce this experimental density.Notice that we do not have access to the very large number of particles used in realistic experimental setups due to the large computational cost involved in our first-principles calculations.Therefore, our scaled system serves to study the physical behavior of the real dipolar gas in a shifted interval of scattering lengths.We have verified that reducing in our simulations the number of particles to N = 200 we recover the same evolution that the one obtained with N = 400 but reducing adequately the scattering lengths.The progressive reduction of a, and thus an increase of the interatomic attraction, with the size of the system is directly related to the existence of a minimum critical atom number for the formation of a self-bound drop.
In the following, we report results obtained with two different confinement configurations that we characterize by their aspect ratio λ = l x /l z , with l i the oscillator lengths along the three spatial directions (i = x, y, z).The first one (λ = 1.83) corresponds to the scaled experimental trap [28], with oscillator lengths l x = 14.84 r 0 and l y = l z = 8.13 r 0 .The second one (λ = 5.00) is more one-dimensional looking, keeping again the central den-sities close to the experimental ones: l x = 29.06r 0 and l y = l z = 5.81 r 0 .
A. Aspect ratio λ = 1.83 First, we study the evolution of the system under the change of the s-wave scattering length a for λ = 1.83.The value of a comes from a delicate balance between the short-range repulsive interaction and the dipolar term.As a general rule, when a > 0 decreases, the system becomes less repulsive as the hard-core radius decreases.In all our simulations we consider the dipoles to be oriented in the z direction due to the action of an external magnetic field.Figure 1 shows the evolution of the system with decreasing a by projecting the position of particles in the x-z plane.As we are dealing with a quantum system, particles are delocalized.This feature is implicit in the path integral method since every particle is represented by a chain or polymer characterized by a given number of beads.The same figure shows typical configurations once the system has reached equilibrium.
Starting from the largest a = 71.1 a 0 , we observe how the particles tend to align in the z direction since dipoles are attractive in head-to-tail configurations, a magnetostriction effect [29].This alignment is however not complete because the gas is harmonically confined, also in that axis.For the largest a considered, we observe the formation of a single self-bound drop, surrounded by a gas.Reducing a the system starts to split in two drops.For a = 58.5 a 0 the two drops are quite close but tend to separate when the scattering length is further reduced.Finally, for a = 45 a 0 , we observe that the two drops become more mutually repulsive and are well separated.In the last configuration, we find two clearly defined drops, with no other particles around them.The reported evolution with the scattering length reproduces qualitatively what is actually observed in experiments, once properly scaled in number of particles and a.
As commented in the preceding section, the use of a microscopic approach such as the PIGS method allows for the calculation of two of the most relevant groundstate magnitudes of the system, the static structure factor and the one-body density matrix.Figure 2 shows results for both magnitudes for a selected set of interaction strengths.The static structure factor is calculated for momenta along the x direction, where drops align.In panel (a), corresponding to a = 60 a 0 , a single peak at k x = 0 appears since only a single drop is present.Moving to panels (b) to (d) more peaks at k x > 0 appear in a progressive way, in agreement with the formation of new drops.It is worth mentioning that the evolution with a of S(k x ) reproduces qualitatively the experimental data on Dy atoms of Ref. [28].The peaks in S(k x ) resemble the Bragg peaks of a real crystal, but they are clearly different since in this case the solid arrangement is formed by drops rather than by single atoms.
The bottom panel of Fig. 2 displays our results for ρ 1 (r).As the system is harmonically confined in all directions, ρ 1 (r) tends to zero when r grows, no matter the direction.The key point is to check whether this function shows a plateau at intermediate distances, where the effect of confinement is less relevant.When more than one droplet appear, this distance is of the order of the separation between drops.Configuration (b) containing two drops is compatible with the presence of ODLRO, with a small plateau at r = 10 -15 r 0 .For this particular interaction, we show the momentum distribution n(k) obtained by Fourier transforming ρ 1 (r) in Fig. 3.The extent of n(k is quite narrow, with a large peak centered al k = 0.When a is further decreased, as in configurations (c) and (d), this plateau is absent and ρ 1 (r) decreases monotonically, so no ODLRO appears.In this way, therefore, our results seem to indicate that the key ingredients of supersolidity can be present, but only in a quite narrow window of interaction strengths.
To localize the range of scattering length values where a supersolid behavior could be observed, we show in Fig. 4 particular values of both ρ 1 (r) and S(k x ).In panel (a), we report the value of the one-body density matrix for two-values of a large distance r max , as a function of a. Panel (b) shows the value of the first peak at k x > 0 of S(k x ), also as a function of a.The strength of the peak in the static structure factor increases when the drops are more dense, and the amount of gas in between decreases.The peak appears when a < 60 a 0 and saturates at a ≃ 55 a 0 .Looking at panel (a) one can see that in this same interval a = 55 -60 a 0 the long-distance value of ρ 1 (r) is stable and large.Therefore, it is within this range of interactions where a superfluid of droplets could exist in our scaled system.
B. Aspect ratio λ = 5.00 The second configuration considered in this work is set to λ = 5, corresponding to a cigar shape geometry.Compared with the previous case, we expect to find more drops, within the same range of interactions.We show the evolution of the system with decreasing a in Fig. 5.At the largest value of a considered (70 a 0 ) the system looks like a gas, with no trace of drops being formed.However, at a = 60 a 0 two drops, immersed in a cloud of gas, are clearly visible.Decreasing a even more, the two drops separate and the amount of gas surrounding them diminishes.Finally, a third drop appears at the lowest value of a considered.In this case, the separation between them is noticeably reduced when compared to the configurations that contain only two droplets, while there is still a small amount of gas around them.
Figure 6 shows results for the static structure factor and the isotropic component of the one-body density matrix for different interaction strengths.As it can be seen, the results are qualitatively similar to the ones found for λ = 1.83.As before, peaks in S(k x ) at k x > 0 appear when drops are being formed, and their number and strength increase when the drops are more dense and well separated.Furthermore, the one-body density matrices decay to zero at sufficiently long distance but differences between them are appreciated at intermediate r.The emergence of a plateau is clearly visible in panels (b) and (c).Within the range of interaction strengths between these two panels, we can state that the system presents simultaneously diagonal and off-diagonal long-range order, the key signals to supersolidity.

IV. CONCLUSIONS
In the present work we have studied a harmonically trapped gas of 162 Dy atoms, with a large magnetic moment, using a first-principles quantum Monte Carlo method.The ground-state of this system has been studied previously using the extended Gross-Pitaevskii equation that incorporates Lee-Huang-Yang terms [9,10,17,18].At difference with the eGPE approach, our description is fully microscopic, relying directly on the Hamiltonian of the system and on a projecting method based on the stochastic evaluation of path integrals [21].In order to make the calculations feasible in practical (computer time) terms, we have scaled the system and used 400 atoms, keeping the same central density of the trap as in experiment [28].Having access to the actual position of the atoms and its delocalization (number of beads in the path integral representation) our results are somehow closer to experimental views [6].
Having access to the atomic coordinates in a spatial representation of the ground state implies that we can calculate its diagonal and off-diagonal properties in an accurate way.The PIGS method provides exact estimates of ground-state properties within some residual statistical noise.We study the evolution of the confined system with a, which decreases when the interaction becomes more attractive.
The diagonal order is analyzed by calculating the static structure factor in a direction x perpendicular to the dipole moments, which are aligned along the z axis.The possible existence of off-diagonal long-range order is studied by estimating the isotropic mode ρ 1 (r) of the onebody density matrix.The simultaneous existence of diagonal and off-diagonal long-range orders are key signatures of a supersolid state.Our work reports the first results on ρ 1 (r) for the three-dimensional dipolar gas and, with some limitations due to the finite size of the system, shows a regime where the system can be found in a superfluid state.The range of interaction strengths, measured in terms of a, where the droplets are coherent is quite narrow, specially when compared with results derived from the extended Gross-Pitaevskii equation.This can be partially attributed to our scaled system, with a number of particles much smaller than the ones in typical experimental realizations.However, the standard eGPE equation, incorporating the first perturbative LHY correction, has shown relevant shortcuts when applied to study dipolar droplets.It was shown in Ref. [28] that the critical atom number, which is the minimum number of atoms required to form a stable drop, is significantly underestimated within that theory.Direct calculations using PIGS did show better agreement with the experimentally measured critical numbers [22].PIGS results point to a system that is globally more attractive than what is predicted by the eGPE, and thus leads to a smaller window of scattering length values where coherent drops are observed.
In future work, we plan to simulate a similar system but in a cylindrical geometry with periodic boundary conditions in the axial direction [30,31].In this way, we expect to be able to evaluate the precise behavior of the one-body density matrix at long distances, making it possible to find a reliable estimation of the condensate fraction.On the other hand, using PIMC we are going to study the effect of finite temperature in the BEC to supersolid transition, which has shown a counterintuitive behavior where order by heating is observed [32,33].
the imaginary-time propagator, and M the number of intermediate imaginary-time steps.

FIG. 1 .
FIG.1.Evolution with the s-wave scattering length of the column densities of the system along the y direction for λ = 1.83.From (a) to (g) a/a0 = 71.1,60, 58.5, 56.6, 54.9, 51.1, and 45.The size of the system in the x direction is 70 r0 for the largest a and reduces with decreasing a until 25 r0 for a/a0 = 45.In the z direction, it extends 50 r0 and squeezes to 40 r0 for the two lowest a values.

FIG. 2 .
FIG. 2. Static structure factor S(kx) for momenta along the x direction(top row) and isotropic mode of the one-body density matrix ρ1(r)(bottom row) for a set of configurations corresponding to λ = 1.83.From a) to d) the scattering length decreases a/a0 = 60, 56.6, 51.1, and 45.In all cases, the error bars are compatible with the width of the lines.The insets show the distribution of particles in the x-z plane.

FIG. 4 .
FIG. 4. Upper panel: value of ρ1(r) at a distance rmax as a function of the s-wave scattering length a. Blue and red points and lines correspond to rmax = 40 and 30, respectively.Lower panel: height of the first peak of S(kx) as a function of a.

FIG. 5 .
FIG.5.Evolution of the dipolar system with the s-wave scattering length for λ = 5.00, with the drops projected in the x-z plane.From (a) to (g) a/a0 = 71.1,60, 58.5, 56.6, 54.9, 51.1, and 45.The size of the system in the x direction is 160 r0 for the largest a and reduces with decreasing a until 90 r0 for a/a0 = 45.In the z direction, it keeps stable around 40 r0.