Quantum Phases of Dipolar Bosons in Bilayer Geometry

We investigate the quantum phases of hard-core dipolar bosons confined to a square lattice in a bilayer geometry. Using exact theoretical techniques, we discuss the many-body effects resulting from pairing of particles across layers at finite density, including a novel pair supersolid phase, superfluid and solid phases. These results are of direct relevance to experiments with polar molecules and atoms with large magnetic dipole moments trapped in optical lattices.

Recent experimental breakthroughs in the realization of ultracold gases of high-spin aligned atoms with large dipole moments [1][2][3][4], highly excited Rydberg atoms [5,6] and of groundstate polar molecules [7,8] hold considerable promise for investigations of many-body quantum systems where dipolar interactions can become dominant [9][10][11][12]. The anisotropy of dipolar interactions combined with the possibility of confining particles in low-dimensional geometries using optical lattices allows for the study of novel pairing mechanisms and the associated quantum phases in a setup where collisional losses are suppressed. This is particularly intriguing for the case of magnetic atoms, where confinement to lattices with spacings as small as 200 nm is possible [13], which favors inter-site dipolar interactions and pairing.
Pairing of two spin-polarized fermionic dipoles across coupled two-dimensional (2D) layers [14] or bosonic one-dimensional (1D) wires [15] in an optical lattice has already led to the prediction of 2D inter-layer superfluidity [16][17][18][19][20][21], analogous to bi-exciton condensation, and the 1D quantum roughening transition [22] in the case of equal number of particles in each layers. Additional exotic phenomena occur for unequal populations [23], where (spin-rotational) symmetry breaking can induce, e.g., stable liquids and crystals of composite multimers [24] for both bosonic and fermionic species. For bosonic gases in the strongly interacting regime [25], emergent parafermionic behavior has been demonstrated [26,27] in coupled 1D wires. In two dimensions, a recent mean-field study in bilayer geometry [28] has predicted novel quantum phenomena for a model of dipolar bosons on a lattice, including a so-called pair supersolid (PSS) phase. Different from supersolids on a single lattice [29][30][31], the latter implies diagonal solid order coexisting with an off-diagonal superfluid order, both derived from composite pairs of dipoles. The experimental observation of this quantum phase and the associated pair superfluids (PSF) and solids would constitute a breakthrough for condensed matter in the cold atomic and molecular context. Thus, the challenge is now to determine whether these quantum phases can be realized for a realistic Hamiltonian representing the microscopic dynamics of strongly interacting dipolar bosons as realized in experiments.
In this work, we study a system consisting of hardcore dipolar bosons confined to two neighboring 2D layers of a 1D optical lattice (see figure 1(a)). The dipole moment of each particle is polarized perpendicular to the layers, which results in repulsive inplane dipole-dipole interactions. This ensures collisional stability against short-range inelastic collisions in the strongly interacting gas. Out-of-plane dipolar interactions are dominantly attractive, which favors inter-layer pairing. Using exact theoretical techniques based on quantum Monte Carlo methods [32], we demonstrate below that this anisotropy and the long-range nature of interactions can induce crystallization of the dipolar cloud into a charge-density wave for a wide range of trapping parameters and interactions. Exotic quantum phases such as the PSS phase and a PSF are achieved under experimentally realistic trapping conditions. These phases can survive up to temperatures of the order of a few nK for a gas of polar molecules or strongly magnetic atoms.
The system we have in mind is described by the single-band tight-binding Hamiltonian Here α, β = 1, 2 and i, j label the layers and the lattice sites in each layer, respectively, while a i,α (a † i,α ) are the bosonic creation (annihilation) operators, with a † 2 i,α = 0 and n i,α = a † i,α a i,α . The brackets denote summation over nearest neighbors only. The first term in equation (1) describes the kinetic energy with in-plane hopping rate J . The second term is the dipole-dipole Here we show one row in each layer. The separation between the layers is denoted by d z , and a is the lattice constant. We use V dd and V ⊥ dd to denote the intra-layer and inter-layer dipolar interactions, respectively. (b) Phase diagram of Hamiltonian equation (1) as a function of V dd /J and particle density n, computed via QMC simulations, for an inter-layer distance d z /a = 0.36 (see the text). CB, checkerboard solid; PSS, pair supersolid; PSF, paired superfluid; 2SF, independent superfluids. The phase boundaries in the dashed region are not resolved. interaction given by V iα; jβ = C dd (1 − 3 cos 2 θ)/(4π |i α − j β | 3 ), where θ is the angle between particles at positions i α and j β and C dd = d 2 / 0 (C dd = µ 0 d 2 ) for electric (magnetic) dipoles of strength d. We denote the repulsive (attractive) nearest neighbor intra-layer (inter-layer) interaction by , with a the in-plane lattice constant. The inter-layer distance is d z . The relative strength V dd /V ⊥ dd can be tuned over a wide range of values by changing d z /a. The quantity µ α is the chemical potential which sets the number of particles in each layer. Here we fix µ 1 = µ 2 , i.e. N 1 = N 2 .
Hamiltonian equation (1) provides a microscopic description for the dynamics of, e.g., a gas of RbCs molecules (d ≈ 1.25 D) at low density n, such that the initial system has no doubly occupied sites [33]. Collisional stability is ensured for n −1/2 (d 2 /hω ⊥ ) 1/3 130 nm with ω ⊥ 100 kHz the frequency of transverse confinement provided by the in-plane optical lattice [34]. In addition, the choice d 2 /d 3 z < V 0 avoids interaction-induced inter-layer tunneling, with V 0 the depth of the optical potential in the transverse direction. Model equation (1) can also be used to describe the dynamics of a gas of strongly magnetic dipolar atoms, such as Dy (d = 10 µ B ). In this case, the conservative estimate above for collisional stability is satisfied for ω ⊥ 1 kHz.
In the following, we present exact theoretical results based on path integral quantum Monte Carlo simulations using a two-worm algorithm [35] which allows for efficient sampling of paired phases. We have performed simulations of L × L = N sites square lattices with L = 8, 12, 16, 20 and 24. For computational convenience, we have set the dipole-dipole interaction cutoff to the third nearest neighbor and have checked that using a larger cutoff did not change the simulation results within errorbars. Lower cutoff values do not allow for stabilization of, e.g., supersolid phases; see below. In the following we choose d z /a = 0.36 and d z ∼ 200 nm, which is experimentally feasible with, e.g., Cr or Dy atoms [2,3]. We show below that this choice allows one to access a parameter regime where particles on different layers can pair up to form a composite object. Below, we first discuss the phase diagram and then discuss in more detail the various phases.
The phase diagram of equation (1) at temperature T = 0 is shown in figure 1(b) as a function of V dd /J and the density n, in the parameter regime 0.31 > V dd /J > 0.2 and 0.1 < n < 0.9, with d z /a = 0.36. We expect this phase diagram to be representative of situations with d z /a 2, where inter-layer pairing is favored (see figure 2).
At half-filling n = 0.5, an incompressible checkerboard solid of pairs (PCB) is stabilized for sufficiently large values of V dd /J . Similar to the conventional checkerboard phase present in single layers [30], here atoms in each layer occupy every other site of the lattice, due to inplane dipolar repulsion. The checkerboard order is characterized by a finite value of the static structure factor S(k) at the reciprocal lattice vector k = (π, π), with and the system displays zero superfluidity. We find that in the PCB phase, atoms across the layers are strongly paired due to attractive inter-layer interactions. As a result, the position of the two checkerboard solids is strongly correlated, i.e. they sit on top of each other. The system can be thus envisioned as a solid of pairs [28,36], with an effective mass m eff ∼ J 2 /(2V ⊥ dd + zV dd ), where z is the coordination number. The PCB solid is stabilized at (much) lower values of V dd /J compared to the case of checkerboard solids in a single layer [30], in analogy with what is found in [36]. This is due to the higher effective mass of the pairs. A similar robustness of this phase is also found for melting at finite temperature.
Upon doping the PCB solid with extra particles or holes, a so-called PSS phase is immediately stabilized. The latter displays both diagonal long-range order with S(π, π) = 0, off-diagonal long-range order associated with a non-vanishing value of the pair-condensate order parameter = a i,α a i,β = 0 (with α = β) and an associated finite superfluid stiffness for pairs (see below). The single-particle condensate order parameter ψ α i = a i,α = 0 is instead zero. The existence of off-diagonal order is consistent with a picture of delocalized defects [29,37], which here correspond to correlated pairs of holes or extra particles across the layers. The PSS phase forms a lobe structure in the (V dd /J − n)-plane, around the PCB line. Away from the tip of the lobe, we find that by varying n at constant V dd /J the PSS loses its diagonal longrange order by melting into a PSF phase, via an Ising-type transition (red continuous line). The PSF phase, with = 0 and ψ α i = 0, is destroyed in favor of a 2SF (a phase with independent, although correlated, superfluids on each layer) for smaller values of V dd /J . In particular, we note that a tiny PSF region should persist in between the PSS and 2SF phases even close to the tip of the PSS lobe; however, this is within errorbars for V dd /J 0.2. Exactly at filling n = 0.5, our results are consistent with a direct PCB-2SF transition, as discussed below, with no intermediate PSS phase. In particular, we find no evidence of, e.g., possible micro-emulsion phases [31,38], within errorbars.
Finally, we note that a host of other phases are present in the general phase diagram for two layers. In particular, we find that for stronger values of V dd /J 0.3 the system displays a sequence of incompressible phases at various rational fillings of the lattice, similar to the socalled Devil's staircase found in the case of a single layer. We also expect novel PSS phases to appear around lobes at, e.g., filling n = 0.25, in analogy with [30]. In addition, independent solids as well as supersolid phases can be achieved by increasing the layer distance, while mixtures of solid and superfluid phases can be stabilized by modifying the relative particle density in the two layers. The discussion of some of these phases is, however, beyond the scope of the present work. In the remainder of the paper, we discuss in more detail the various phases and their transitions at zero and finite temperature around n = 0.5.
Stability of the PCB phase. As discussed above, the PCB phase at n = 0.5 is characterized by a finite value of the order parameter S(π, π) and no off-diagonal order. The latter is associated with superfluidity in a (2+1)-dimensional interacting system, which can be measured straightforwardly within Monte Carlo (see below). In addition, within the PCB phase inter-layer dipolar attraction strongly correlates the positions of particles in the two layers.
The stability of the PCB phase with respect to intra-plane interactions as well as inter-layer distance d z /a at zero temperature is analyzed in figure 2. There, we numerically determine the minimum dipolar interaction strength V dd /J required to stabilize the PCB phase at a given d z /a. In order to establish whether the solid phase is paired, we have performed several simulations with different initial conditions for each set of parameters and observed whether the equilibrium configuration was dependent on the initial choice or not. The figure shows that a PCB phase is stabilized for d z /a 2 and sufficiently large V dd /J (continuous line). In this parameter regime, the system above (below) the continuous line is a PCB (2SF) phase, respectively; that is, the continuous line visualizes the shift of the PCB-2SF transition point of figure 1(b) as a function of d z /a. Instead, for d z /a > 2 and large enough interactions the insulating phase above the (dotted) line corresponds to two independent checkerboard phases (2CB). This points to the possible presence of a tri-critical point in the phase diagram around d z /a ≈ 2. The shaded area in figure 2 shows the transition region where the PCB phase is replaced by the 2CB phase. The details of this phase transition are beyond the scope of this work. We have confirmed that the computed transition points are independent of the interaction cutoff that we use, within our errorbars, and should thus be quantitatively relevant to experiments. In the following, we focus on d z /a = 0.36 to satisfy V ⊥ dd 10J in the vicinity of the tip of the lobe. We find that this choice ensures pairing at n ∼ 0.5 (in the vicinity of the PCB phase) while keeping V dd relatively low. This corresponds to experimentally optimal conditions to observe PSS phase: a lower effective mass of pairs m eff results in a larger superfluid density, which in turn results in higher critical temperatures (see also below).
PSS phase. Figure 1(b) shows that a PSS lobe is immediately formed by doping the PCB solid with either vacancies (holes) or interstitials (extra particles). The hard-core constraint of equation (1) ensures particle-hole symmetry, and thus reflection symmetry of the lobe, around n = 0.5.
We characterize this PSS phase in figure 3, for a specific choice of interaction strength V dd /J = 0.238. In the figure, the order parameter for the diagonal checkerboard solid order S(π, π) (continuous lines) and the superfluid stiffness of pairs ρ PSS are plotted as a function of n. The quantity ρ PSS = T W 2 /d L d−2 [39] is directly related to a pair condensate, and can be calculated within quantum Monte Carlo, with W = W 1 + W 2 the sum of winding numbers in layers 1 and 2. The figure shows that for an extended range of densities, both the static structure factor and the PSF stiffness are finite and system size independent, showing the existence of a stable supersolid phase in the lobe region. We note that, due to pairing across the layers, in the PSS phase the fluctuation of difference in winding numbers is zero (W 1 − W 2 ) 2 .  figure 1 and the corresponding order parameters: structure factor S(π, π ); single-particle condensate ψ α i = a i,α in each layer α; pair-condensate order parameter = a i,α a i,β , with α = β (see the text).
Superfluid phases. As the system is doped further, the PSS disappears in favor of a PSF phase. The latter displays pair-induced off-diagonal long-range order only (see, e.g., table 1). We find that the PSS-PSF transition is of the Ising-type universality class in (2+1) dimensions, analogous to the case of a single layer [30]. Critical points are determined using finitesize scaling for the static structure factor with scaling coefficients 2β/ν = 1.0366 [40] (see figure 3(b) for the specific choice V dd = 0.238J ). In the figure the scaled quantity S(π, π)L 1.0366 is plotted as a function of n, and the crossing of the curves at n cr = 0.573 ± 0.002 corresponds to the quantum critical point where the finite-size effects disappear (see also panel (a)).
We find, in general, that by lowering the interaction strength V dd /J at constant n from the PSF phase, the system finally develops into two independent superfluids (2SF) with a finite value of the single-component condensate order parameters, ψ α i = ψ β i = 0, via a second-order phase transition in the (2+1) XY universality class. The transition points between the PSF and 2SF phases in figure 1(b) are calculated using finite-size scaling of (W 1 − W 2 ) 2 . The latter quantity is zero inside the PSF phase in the thermodynamic limit due to pairing across the layers, while it has a finite value in the 2SF phase. We note that the pair-order parameter in the 2SF phase is instead trivially non-zero, = 0 (see also table 1).
The phase diagram in figure 1(b) shows that the boundary of the PSF-2SF transition shifts downward approximately linearly in the (V dd /J − n)-plane, as the density becomes sufficiently smaller or larger than n = 0.5. This is easily understood in the limit of very small densities, by noting that inter-plane dipole-dipole interactions always favor the existence of a two-body bound state, even for an arbitrarily small interaction strength. However, we find that many-body effects result in a threshold for the formation of pairs at finite density, where the magnitude of the interaction strength required to stabilize pairing increases with n. This is explained by noting that, in the limit of low density, d z √ n 1, the PSF phase is composed of weakly interacting superfluid dimers. When the density is increased, exchanges between dimers are favored. This destabilizes the dimers, inducing the transition to two independent superfluids in the 2SF phase. Eventually, the presence of diagonal order near n = 0.5 forces the PSF-2SF line to bend down, deviating from the linear dependence on n.
We gain further insight into the structure of correlations in the condensed phases by studying the following four-point correlation function: Here i, j, l refer to sites, 1, 2 refer to layers and denotes a quantum and thermal average as well as site averaging over i. In the presence of pair superfluidity, one expects this correlation function to be short ranged with respect to r jl = |r j − r l |, and simultaneously long ranged with respect to r il = |r i − r l | and r i j = |r i − r j |. In the 2SF phase, instead, f jl is obviously long ranged with respect to r il and r i j , but it is independent of r jl . Figure 4 shows f jl (normalized to unity: ∞ 0 f jl d r jl = 1) as a function of r jl for the PSS (green triangles, n = 0.48, V dd = 0.25J ), PSF (red dots, n = 0.40, V dd = 0.25J ) and 2SF phases (blue squares, n = 0.30, V dd = 0.18J ). As expected, we find that f jl is independent of r jl in the 2SF phase, where pairing is absent, whereas it is peaked at r i j = 0 in both the PSS and PSF phases. The figure shows that an exponential ansatz of the form f 0 e −r jl /ξ 0 fits quite well the large-r i j behavior of f jl in the latter phases, and is essentially exact for all r jl in the PSS phase with ξ 0 = 1.63a. Here, ξ 0 can be interpreted as the spread of the pair wavefunction, and is obtained from figure 4 by fitting the tail of f jl , as obtained numerically. The inset in figure 4 shows ξ 0 as a function of n, as the PSF-PSS phase boundary is crossed. The pair wavefunction is shown to be considerably more tightly bound in the PSS phase than in the PSF phase. The abrupt drop in ξ 0 locates precisely the transition point.
Finite temperature. We have studied the robustness of the quantum phases described above against thermal fluctuations. As expected for 2D systems, we find, in general, that superfluidity in the PSS, PSF and 2SF phases disappears at finite temperature T via a Kosterlitz-Thouless (KT)-type [41] transition. Diagonal long-range order in the PCB and PSS phases is instead lost via a 2D Ising-type transition. We found that, when present, pairing still exists at the transition shown using black squares, red circles, blue triangles and green diamonds, respectively. When the temperature is increased, the 2SF phase undergoes a KT phase transition at the critical temperature T KT, 2SF ≈ 0.255J , indicated by an arrow. The inset shows finite-size scaling [42] where the dashed line is a linear fit of our simulation results (points).
points, suggesting that the temperatures required for breaking pairs are higher than the critical temperatures measured here. Figure 5 shows one example of the SF-normal transition in the 2SF phase. We plot ρ s versus T /J at V dd /J = 0.20 and n = 0.3 for different system sizes. The inset shows the finite-size scaling procedure [42] used to determine the critical temperature. We find that T KT, 2SF = πh 2 ρ s (T KT ) 2 ≈ 0.255J . For the PSF-normal transition we find that T KT, PSF ≈ 0.08J at n = 0.3 and V dd /J = 0.25. The lower KT transition temperature compared to the 2SF-normal transition is due to a larger effective mass of the pairs, i.e. lower effective hopping, which results in a suppression of particle delocalization and consequently smaller ρ s . The disappearance of the PSS phase proceeds in two successive stages. At T KT, PSS the PSS phase melts into a liquid-like phase reminiscent of a liquid crystal, with ρ s = 0 and S(π, π) = 0. Upon further increasing the temperature S(π, π ) becomes zero at a critical temperature T c through an Isingtype transition (2β/ν = 1/4 in 2D). For example, we find that T KT, PSS ≈ 0.06J and T c ≈ 0.3J for n = 0.48, V dd = 0.25J . Similar T c values are found for the critical temperature of the melting of the PCB phase into a featureless normal fluid, e.g. T c ≈ 0.35J for V dd = 0.25J . Clearly, for larger interaction strengths, i.e. away from the tip of the lobe, transition temperatures will increase.
Experimental estimates. Based on our results we estimate under which experimental conditions the phases described can be observed. For example, with a gas of Dy (d = 10µ B ) a choice of lattice parameters a = 500 nm, d z = 200 nm, J = 50 h Hz results in V dd /J ∼ 0.21, which stabilizes the PCB phase. In the case of Er 2 Feshbach molecules [43,44] (d = 14µ B ) with a = 400 nm, d z = 200 nm, J = 100 h Hz, the PCB phase is stabilized at V dd /J ∼ 0.4. In both cases, the PCB phase can be observed at nK temperatures.
Using RbCs (d = 0.3D) and typical trapping parameters a = 500 nm, d z = 300 nm and J = 150 h Hz, we find that V dd /J = 0.7, which is large enough to stabilize the PCB. The latter survives up to T PCB c ∼ 4 nK. By doping away from filling factor n = 0.5 the PSS phase can be reached with a KT transition temperature for PSF-normal transition of the order of nK.
In local density approximation, the presence of weak in-plane harmonic confinement as provided by, e.g., magnetic trapping will essentially result in a scan of the phase diagram over the chemical potential, with coexistence of different phases. Similar to the single-layer case [30], because of their energy gaps, we expect the solid phases CB and PCB to be robust against, e.g., hole doping, and thus observable, in the presence of in-plane harmonic confinement.
All the phases discussed above can be distinguished by monitoring the behavior of the order parameters of table 1. Inter-layer particle correlations will be detected via in situ imaging [45,46] as well as noise correlation measurements. In addition, inter-layer pairing may be measured via spectroscopic techniques.
In conclusion, we have studied the quantum phases of dipolar bosons in a bilayer lattice geometry described by the microscopic Hamiltonian equation (1) for hard-core particles, in a situation where the number of particles in each layer is the same. Relevant to experiments with polar molecules and magnetic atoms, we have established under which conditions pairing for two particles is stabilized across the layers. Our zero-temperature study indicates that the system displays a rich ground state phase diagram including a novel PSS phase for hard-core dipolar bosons, in addition to PSF and checkerboard-like solid phases. Our finite-temperature results indicate that these phases are experimentally observable at temperatures of the order of nK. A four-body correlation function connected with the spread of the pair wavefunction can be used to characterize these phases and their transitions. Future work will include the extension of similar quantum Monte Carlo studies to multilayer geometries as well as to systems with population imbalance in the layers.