Patterned Supersolids in Dipolar Bose Systems

TL;DR

This paper uses Quantum Monte Carlo simulations to explore the phase diagram of dipolar bosons, revealing a broad supersolid phase with various crystalline patterns and a roton minimum near crystallization.

Contribution

It provides the first detailed phase diagram of dipolar bosons including a broad supersolid region influenced by interaction range and density.

Findings

Supersolid phase exists over a broad parameter range.

Crystalline patterns depend on density and interaction range.

A roton minimum appears near the transition to crystallization.

Abstract

We study by means of first principle Quantum Monte Carlo simulations the ground state phase diagram of a system of dipolar bosons with aligned dipole moments, and with the inclusion of a two-body repulsive potential of varying range. The system is shown to display a supersolid phase in a relatively broad region of the phase diagram, featuring different crystalline patterns depending on the density and on the range of the repulsive part of the interaction (scattering length). The supersolid phase is sandwiched between a classical crystal of parallel filaments and a homogeneous superfluid phase. We show that a "roton" minimum appears in the elementary excitation spectrum of the superfluid as the system approaches crystallization.

Equations

$$\hat{H}=-\frac{1}{2}\sum _i{\nabla }_i^2+\sum _{i<j}U({\mathbf{r}}_i,{\mathbf{r}}_j)$$

$$U_{sr}(r)=\frac{4\pi {\hslash }^2{a}_s}{m}\delta (\mathbf{r})$$

$$U_{sr}(r)=(\sigma /r{)}^{12}$$

$$\frac{a_s}{\sigma }\approx 0.76{\sigma }^{1/5}$$

$$U_d(\mathbf{r},{\mathbf{r}}^{\prime })=\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^3}\left(1-\frac{3(z-z^{\prime }{)}^2}{|\mathbf{r}-{\mathbf{r}}^{\prime }{|}^2}\right)$$

$$e(q)=\frac{q^2}{2S(q)}$$

Full text

1 Introduction

Dipolar gases are the focus of extensive experimental and theoretical research, motivated by the possibility that yet unobserved, exotic phases of matter may be underlain by the distinctive character of the inter-particle interaction, both long ranged and anisotropic menotti . In particular, the experimental achievement of Bose-Einstein Condensation of atomic systems with large magnetic moments griesmeier ; lu ; aikawa ; depaz ; ni ; yan ; takekoshi ; park ; balewski suggests that one might be able to predict and observe phases featuring more than one type of order. Of particular interest is the supersolid phase, which spontaneously breaks both continuous translational and global U(1) symmetries, thus featuring simultaneously crystalline order and flow without dissipation (see, for instance, Ref. rmp ).

As supersolid behavior in a crystal of 4He, once deemed the most promising candidate, has so far eluded unambiguous observation, dipolar bosons have been suggested as a likely physical setting for its detection in various experimental wentzel ; ferlaino and theoretical cinti10 ; jltp ; cinti works. In Ref. cinti , for instance, it is contended that, if dipole moments are aligned, an ordered array of filaments (or, prolate droplets) constitutes the ground state of the system in rather broad conditions, and preliminary evidence of global phase coherence among such droplets was offered for specific values of density and interaction parameters.

It seems now accepted that the appearance of filaments in this system is the result of the competition between the attractive part of the dipolar interaction, and the presence of a short range repulsion, which prevents the system from collapsing. Such a repulsive part is often modeled theoretically through the so-called scattering length approximation. Because the scattering length is experimentally controllable and can be varied by means of the Feshbach resonance (see, for instance, Ref. feshbach ), the possibility arises of exploring the quantum phase diagram of the system in its entirety, at least within known limitations (e.g., three-body recombination) fb ; pfau2 ; ferlaino3 ; ferlaino2 . A comprehensive theoretical study of the bulk phase diagram of the system at low temperature based on reliable, first principle computational methods, aimed at helping in the design and interpretation of present and future experiments, seems therefore timely lastnote .

In this work, we carry out state-of-the-art Quantum Monte Carlo (QMC) simulations to investigate the low temperature phase diagram of dipolar bosons with aligned dipole moments. We model the repulsive part of the interaction by means of an inverse power law potential, as done in previous work cinti ; a straightforward connection exists between the characteristic range $\sigma$ of this interaction, and the scattering length $a_{s}$. We map out the phase diagram by computing relevant correlations, as well as the superfluid density, all directly accessible in our numerical approach, as a function of particle density and $\sigma$.

The system displays several distinct phases, ranging from an essentially classical crystal of parallel, particle-thin filaments in one limit ($\sigma\to 0$), to a hard-sphere-like superfluid, reminiscent of liquid 4He in the opposite ($\sigma\to\infty$) limit, at any particle density. An intermediate region exists between these two phases, in which the ground state displays both crystalline order and a finite superfluid response. This supersolid phase exists within a range of scattering length that depends monotonically on the density.

The crystal structure of the supersolid is quite distinct from that of the classical crystal, determined by quantum-mechanical effects, both zero-point motion as well as exchanges of particles (as already noted in Ref. cinti ). Remarkably, the supersolid phase features different arrangements of particles (patterns) in different regions of the phase diagram. At low density, the supersolid phase consists of a crystal of prolate droplets, with quantum-mechanical exchanges of particles across droplets, as suggested in Ref. cinti . As the density is increased, on the other hand, the patterns resemble some of the periodic structures originally predicted to occur in two-dimensional dipolar systems vanderbilt ; spivak .

We study the elementary excitation spectrum for the superfluid phase, which features the experimentally observed “roton” minimum ferlaino , as crystallization is approached from the superfluid side, both on reducing $\sigma$ (i.e., moving toward the supersolid phase) as well as on increasing it, in which case the system behaves essentially as a hard sphere fluid, transitioning into a conventional (i.e., non-superfluid) crystal.

The remainder of this paper is organized as follows: in section 2 we describe the model of the system; in Sec. 3 we briefly describe our methodology; we present and discuss our results in Sec. 4 and finally outline our conclusions in Sec. 5, where we discuss the relevance of this study to recent experimental work.

2 Model

The system is modeled as an ensemble of $N$ identical particles of mass $m$ and spin zero, hence obeying Bose statistics. These particles possess a magnetic moment $d$, pointing in the $z$-direction. We are interested in studying the phase diagram of the bulk; thus, unlike in a typical experiment, we do not confine the simulated system by means of an external potential. Rather, our system is enclosed in a three dimensional box, shaped like a cuboid of volume $V$, with periodic boundary conditions in the three directions. The shape of the cell was varied, depending on the particular structures and patterns forming at the various physical conditions. We take the characteristic length of the dipolar interaction, $a\equiv md^{2}/\hbar^{2}$ as our unit of length, and $\epsilon\equiv\hbar^{2}/(ma^{2})$, as that of energy and temperature. The quantum-mechanical many-body Hamiltonian in dimensionless form reads as follows:

[TABLE]

where the first (second) sum runs over all particles (pairs of particles), and the pair potential consists of two parts, $U({\bf r},{\bf r}^{\prime})=U_{sr}(|{\bf r}-{\bf r}^{\prime}|)+U_{d}({\bf r},{\bf r}^{\prime})$, $U_{sr}$ being the repulsive part. As mentioned above, in most theoretical studies the repulsive part of the interaction is modeled by means of the so-called scattering length approximation, namely

[TABLE]

To the extent that such a representation is valid, expression (2) can be replaced by any other potential that has the same scattering length $a_{s}$. In this work, we use for $U_{sr}$ the repulsive part of the standard Lennard-Jones potential, i.e.,

[TABLE]

whose use is more convenient in numerical simulations. The parameter $\sigma$ of the potential $U_{sr}$ used here is directly related to the scattering length $a_{s}$, through

[TABLE]

(see, for instance, Ref. flugge ). $U_{d}$ is the classical dipolar interaction between two aligned dipole moments, namely

[TABLE]

At zero temperature there are two parameters that govern the physics of the system, namely the particle density $\rho\equiv N/V$, and the characteristic range $\sigma$ of the repulsive interaction. The classical limit is approached as $\sigma\to 0$, whereupon the attractive well of the anisotropic interaction becomes deep enough to dominate the physics, and quantum mechanical effects are small. On the other hand, as $\sigma$ is increased, quantum mechanical effects are increasingly prominent, eventually destabilizing the classical ground state and giving rise to interesting physics.

3 Methodology

We carry out QMC simulations of the system described in section 2 using the continuous-space worm algorithm worm ; worm2 . We shall not review the details of this method, referring instead the reader to the original references. We utilized a canonical variant of the algorithm in which the total number of particles $N$ is held constant, in order to simulate the system at fixed density mezz1 ; mezz2 . Although this is a finite temperature technique, we perform simulations at sufficiently low temperatures, so that computed physical properties are essentially those of the ground state.

We survey the phase diagram of the system by performing simulations at a fixed density for different values of $\sigma$, and then repeating the process for different values of the density. We performed simulations of systems of sizes ranging from $N=160$ to $N=648$ particles, and densities between $\rho=0.125$, which is close to the value of some of the current experiments ferlaino , up to three orders of magnitude higher, i.e., $\rho=100$, quite likely not attainable in present time experiments but nonetheless of fundamental interest. To ensure that the physics is independent of the initial configuration of the particles, most simulations were started from high-temperature, disordered configurations.

Details of the simulation are standard; we made use of the primitive approximation for the high-temperature density matrix. Although more accurate forms exist, we found that in practice with this particular interaction the primitive approximation is the most efficient. All of the results quoted here are extrapolated to the limit of time step $\tau\to 0$; quite generally we found numerical estimates for structural and superfluid properties of interest here obtained with a value of the time step $\tau\sim 10^{-3}\epsilon^{-1}$ to be indistinguishable from the extrapolated ones, within the statistical uncertainties of the calculation.

Occurrence of crystalline order in the system can be detected through the calculation of the static structure factor $S({\bf q})$; because of the anisotropy of the interaction and the ensuing tendency of the system to form filaments along $z$ (the direction of alignment of the dipole moments), we computed $S({\bf q})$ for q-vectors lying in the $xy$ plane. As mentioned above, distinct, characteristic patterns form for different values of the parameters $\rho,\sigma$. The identification of the various patterns is achieved through visual inspection of the configurations generated by the algorithm in the course of sufficiently long computer runs.

The superfluid response of the system is assessed through the direct calculation of the superfluid fraction $\rho_{S}(T)$, by means of the well-known winding number estimator pollock . Due to the anisotropic nature of the interaction, and consequently of the superfluid (and supersolid) phases observed here, we offer results for the in-plane ($xy$) superfluid response only, throughout this paper. A typical result for the in-plane superfluid fraction $\rho_{S}$ is shown in Fig. 1b.

4 Results

4.1 Phase diagram

The ground state phase diagram of the system, as it emerges from our extensive QMC simulations, is shown in Fig. 1a in the $\rho-\sigma$ plane. Here we outline some of its generic features, offering a more detailed discussion of the two most interesting phases (the supersolid and the superfluid) in subsections 4.2 and 4.3.

For any fixed value of the density, two clear physical limits can be easily identified, more intuitively discussed in terms of the average interparticle distance $b\equiv\rho^{-1/3}$. Specifically, if $\sigma<<b$ the ground state of the system consists of an ordered array (a triangular lattice) of thin filaments oriented along the $z$-direction. An example is shown in Fig. 2a, displaying the density of the system, integrated along the $z$-direction morepr . In this physical limit, the depth of the attractive well of the dipolar interaction causes the potential energy to dominate the behavior of the system, which can be understood and predicted quantitatively along classical lines, as shown in Ref. cinti . Specifically, the system forms an ordered array of particle-thin parallel filaments, arranged on a triangular lattice. In this case, exchanges of identical particles are suppressed, both within a filament, as well as across filaments.

Conversely, when $\sigma$ becomes of the order of the interparticle distance $b$, the attractive well of the dipolar interaction weakens, and the physics of the system morphs into that of a hard-sphere fluid, as the repulsive part of the interaction becomes the dominant feature. The behavior of the system in this limit is very similar to that of superfluid 4He, which undergoes “conventional” crystallization as $\sigma\sim b$ (this part of the phase diagram is not shown in Fig. 1a).

As shown in Fig. 1a, there is in an intermediate range of $\sigma/b$ wherein the system displays the most interesting, novel behavior, specifically a supersolid phase, which we discuss in detail in subsection 4.2. A remarkable feature of this phase diagram is that there is no way of going from the superfluid to the insulating crystalline phase at zero temperature, by varying $\rho$ and/or $\sigma$, without going through a supersolid phase. This is also the case in a 2D system of soft core bosons saccani .

4.2 Supersolid

The supersolid phase is generally characterized by the formation of relatively large, prolate droplets, elongated along the $z$-direction, arranging themselves on a triangular lattice or forming more complex structures, as we discuss below. Quantum mechanical exchanges of identical particles, largely suppressed in the classical filament crystal, become important, initially within a single filament super and progressively across filaments, establishing phase coherence and leading to a finite, three-dimensional superfluid response throughout the system. It is interesting to note that while the superfluid response is anisotropic, and in the $xy$ plane saturates to a value lower than unity as $T\to 0$ (see Fig. 1b), as expected leggett , it is always seen to approach unity, in the same limit, in the direction of elongation of the droplets ($z$).

The range of values of $\sigma$ within which the system displays supersolid behavior depends on the density of the system, as shown in Fig. 1a; in particular, at high density the supersolid occurs for lower values of $\sigma$, as quantum-mechanical exchanges are favored by a lower interparticle distance. The width in $\sigma$ of the supersolid region, however, appears to be roughly independent of $\rho$.

As mentioned above, the occurrence of an ordered arrangement can be established through the calculation of the static structure factor, which develops (Bragg) peaks in correspondence of relevant wave vectors (e.g., the inverse distance between droplets). Additional information comes from the direct visual inspection of the patterns that form in the course of the simulation, which are quantum-mechanical in nature and can be markedly different, depending on both the density and the value of $\sigma$.

Fig. 2b shows the same type of density map as in Fig. 2a, but for a supersolid system. As one can see, the filaments are in this case replaced by larger droplets, which include considerably more particles than the classical filaments important . Quantum-mechanical exchanges among adjacent droplets are frequent at low temperature (i.e., $T\sim 1$ for this particular choice of $\rho,\ \sigma$); consequently, macroscopic exchange cycles (i.e., comprising nearly all the particles in the system) take place. Clearly, the density map shown in Fig. 2b is very reminiscent of those that appear in a purely two-dimensional (2D) system of soft core bosons saccani ; saccani2 , namely the simplest supersolid. Indeed, the physics of the system under study here, projected onto the plane perpendicular to the filaments, could be regarded as equivalent to that of a 2D soft core system, the third dimension serving the purpose of “piling up” particles, thus allowing for the formation of “cluster” unit cells.

It is quite interesting to note that the supersolid phase does not always display the structure shown in Fig. 2b. Rather, at high particle density patterns begin to emerge, such as the inverted droplet structure shown in Fig. 3a, and the striped one cintiboninsegni of Fig. 3b, or others that are evocative of those predicted for 2D dipolar systems, in the context of “microemulsions” vanderbilt ; spivak

4.3 Superfluid

The superfluid phase arises as the range of the repulsive interaction $U_{sr}$ (eq. 3) is progressively increased, at fixed density, as shown in Fig. 1a. As the droplets expand, due to the hard core repulsion, they eventually merge, giving rise to a uniform superfluid phase. It is interesting to study the elementary excitation spectrum of the superfluid phase, especially as crystallization is approached.

Although it is possible to infer the excitation spectrum from the full imaginary time dynamics, computed by QMC, using an inversion method (e.g., MaxEnt ceperley ; maxent ), because we are only interested here in gaining qualitative understanding we make use of a simpler approach, based on the Bijl-Feynman approximation feynman54 . Specifically, we assume that the dynamic structure factor $S({\bf q},\omega)$ is dominated by a single peak, which allows one to obtain the elementary excitation spectrum as

[TABLE]

Eq. 6, known as Bijl-Feynman formula, provides a reasonable qualitative and semi-quantitative account of the elementary excitation spectrum in superfluid 4He, in particular of the presence and position of the roton minimum glyde .

The behavior of the excitation spectrum as a function of $\sigma$ is the same across all values of density investigated in our simulations. Fig 4 shows two examples, for $\rho=1.0$ and $\rho=100$. The phonon-like dispersion that is characteristic of the superfluid phase at low $q$ starts acquiring a negative curvature as the superfluid-supersolid transition is approached from above (with reference to the phase diagram shown in Fig. 1a. At values of $\sigma$ sufficiently close to the critical value, the system experiences stronger density fluctuations at that value of the momentum, giving rise to a roton minimum and signaling incipient crystallization.

These results are at least in qualitative agreement with the experimental data in Ref. ferlainopreprint . As the value of $\sigma$ is lowered, the roton minimum becomes progressively lower, eventually hitting the horizontal axis; that is consistent with the divergence of $S(q)$ at the roton wave vector, i.e., the formation of a dipolar crystal. We come back to a more extensive comparison of our results with experiment in the next section.

It is important to note that the presence of a roton minimum in the elementary excitation spectrum merely signals the proximity of the system to crystallization. As such, it is not a special feature of this particular system, or of the character of the interaction (dipolar). Indeed, a roton minimum can be observed in the large $\sigma$ limit, in which the system is essentially a gas of hard spheres, as the crystal (in this case a conventional one, i.e., with few particles per unit cell) is approached from below. In the latter case, the roton minimum forms at a different, higher value of the wave vector, which reflects the fact that the ensuing crystal structure in this case has a smaller lattice constant.

5 Conclusions

We carried out extensive QMC simulations of a system of dipolar Bose particles of spin zero, in three dimensions, in order to gain insight into the phase diagram of this system. The anisotropic character of the interaction gives rise to novel phases, chiefly a supersolid characterized by various, intriguing density patterns.

The supersolid phase appears to be observable in a relatively broad range of parameter space, obviously making allowance for the difficulty of exploring the high density phase, as three-body recombination must be overcome.

The results presented here are consistent with recent experimental data ferlainopreprint ; tanzi ; pfaupreprint ; newferlaino . In order to assess whether the density range explored here is comparable to that of recent experiments, we consider Ref. newferlaino , reporting measurements carried out on an assembly of 164Dy atoms. Using their numbers, namely $N=3.5\times{10^{4}}$ atoms confined in an anisotropic harmonic potential of characteristic frequencies equal to 300, 16 and 22 Hz, we estimate the density at the center of the trap $\sim 1.3\times 10^{-8}$ Å*-3*. Expressed in units of the dipolar length $a$, which in this case is worth 208 Å factorofthree , this is equal to 0.117, i.e., very close to the lowest density considered in this work. Indeed, the results shown in Fig. 1g of ref. newferlaino indicate that the transition between a uniform BEC and a supersolid occurs in correspondence of a scattering length $a_{s}\sim 0.24\ a$, which is in remarkable quantitative agreement with our prediction of $\sigma\sim 0.4\ a$ for $\rho=0.125$ (Eq. 4 yields $\sigma=0.38$ if $a_{s}=0.24$).

This result provides strong quantitative support for the microscopic model utilized here, as well as for the results of our calculation. It also suggests that the experimental findings of Ref. newferlaino reflect the physics of the bulk, to an appreciable extent, i.e., they are not overly affected by confinement, nor by the relatively small size of the system. This observation is consistent with the conclusions of a recent theoretical (mean-field) study of the system confined in an elongated trap ancilotto , whose results are in qualitative agreement with ours.

Acknowledgments

This work was supported by the Natural Sciences and Engineering Research Council of Canada, as well as ComputeCanada. Useful conversations with F. Cinti, S. Moroni and F. Ferlaino are also gratefully acknowledged.

Conflict of interest

The authors declare that they have no conflict of interest.