Dipole-dipole interactions in optical lattices do not follow an inverse cube power law

We study the effective dipole-dipole interactions in ultracold quantum gases on optical lattices as a function of asymmetry in confinement along the principal axes of the lattice. In particular, we study the matrix elements of the dipole-dipole interaction in the basis of lowest band Wannier functions which serve as a set of low-energy states for many-body physics on the lattice. We demonstrate that the effective interaction between dipoles in an optical lattice is non-algebraic in the inter-particle separation at short to medium distance on the lattice scale and has a long-range power-law tail, in contrast to the pure power-law behavior of the dipole-dipole interaction in free space. The modifications to the free-space interaction can be sizable; we identify differences of up to 36% from the free-space interaction at the nearest-neighbor distance in quasi-1D arrangements. The interaction difference depends essentially on asymmetry in confinement, due to the d-wave anisotropy of the dipole-dipole interaction. Our results do not depend on statistics, applying to both dipolar Bose-Einstein condensates and degenerate Fermi gases. Using matrix product state simulations, we demonstrate that use of the correct lattice dipolar interaction leads to significant deviations from many-body predictions using the free-space interaction. Our results are relevant to up and coming experiments with ultracold heteronuclear molecules, Rydberg atoms, and strongly magnetic atoms in optical lattices.


Introduction
Recent experimental progress in cooling heteronuclear polar molecules with large electric dipole moments [1,2,3,4,5,6,7,8,9,10,11], Rydberg atoms [12], and atoms with large magnetic dipole moments, in particular Chromium [13], Erbium [14], and Dysprosium [15,16], has sparked interest in the properties of ultracold dipolar gases. While in many ultracold atomic systems interactions are short-range and wellmodeled by a contact pseudopotential, the interactions in dipolar gases have a longrange and anisotropic character in free space, decaying as 1/r 3 with the separation r between particles. These features of the dipole-dipole interaction have lead to a variety of intriguing theoretical proposals such as exotic pairing and bound states in ladder geometries [17,18] and the realization of quantum liquid crystal states of matter [19,20,21]. Even for atoms with relatively weak dipole moments, such as Rubidium, dipole-dipole interactions can play a significant role [22,23]. In this article, we show that the effective dipolar interaction in a lattice is not actually 1/r 3 as commonly believed, but has a non-algebraic decay at moderate separations, and only behaves as 1/r 3 for large separations. Corrections of order 36-48% arise for interactions at the nearest-neighbor distance in moderately confined quasi-low-dimensional scenarios.
A key component of our analysis is the presence of a continuous, periodic potential, which for ultracold atomic and molecular gases is provided by an optical lattice. As first discovered by Kohn [24], Wannier functions, the most localized set of orthogonal singleparticle states with the symmetries of the lattice, generally feature an exponential decay. We find that the effective dipolar interaction in an optical lattice depends essentially on the exponential tails of the Wannier functions rather than only on their widths. Hence, approximating the Wannier functions with localized functions which match only the mean width will fail to accurately capture the effective interaction in an optical lattice.
An additional essential ingredient for our findings is an asymmetry in the degree of confinement along the principal axes of the lattice due to the anisotropic character of the dipole-dipole interaction. In this work, we characterize confinement using the curvature of a lattice site minimum. Our work builds on a wealth of confinement-induced phenomena in ultracold quantum gases, such as confinementinduced resonances [25,26,27,28], the fermionization of a 1D Bose gas [29], and the Berezinskii-Kosterlitz-Thouless transition in a quasi-2D Bose gas [30,31]. We stress, however, that our work does not deal with dipolar confinement-induced resonances such as those studied in Ref. [32], but rather on the modification of the effective interaction in an optical lattice due to the localization properties of the single-particle basis. In dipolar gases, the effects of confinement have been studied in harmonic traps [33,34,35], and within the Gross-Pitaevskii approximation in a triple-well potential [36] and a 1D lattice [37]. For the harmonic oscillator, it has been shown that strong confinement along the axis of a field orienting the dipoles leads to purely repulsive interactions in the weakly confined plane. Additionally, because of the anisotropic character of the dipole-dipole interaction, the stability of a dipolar BEC displays a strong dependence on anisotropy in external confinement, and stable solutions can even take surprising forms such as the "red blood cell" dipolar BEC [38]. Similarly, anisotropic lattice confinement has surprising differences from preconceptions based on uniform isotropic systems. Dipolar interactions in confined geometries appear in many other branches of physics, such as in ferromagnetic nanostructures [39], where the dipole-dipole interaction plays a key role in the dispersion relation for spin waves in ferromagnetic films [40] and wires [41].
An important application of our results is in deriving many-body models to describe the low-energy physics of dipolar gases on lattices. Previous derivations, for example Refs. [42,43,44,45,46], assume that the interaction between localized lattice states has the same functional form as in the continuum. Since this amounts to replacing the localized single-particle probability distributions with delta functions, we will call this the delta function approximation (DFA). Performing matrix product state simulations [47] on infinite one-dimensional (1D) lattices, we demonstrate that the DFA can lead to significant errors in the determination of the phase diagram.
This paper is organized as follows. In Sec. 2 we review the theory of Wannier functions and their use in deriving effective many-body lattice models for strongly correlated systems. In particular, we provide a quantitative analysis of the decay properties of Wannier functions as well as properties of their squares interpreted as probability distributions. Sec. 3 provides numerical results for the effective dipole-dipole interactions in the presence of an optical lattice and discusses the effects of confinement. In Sec. 4 we study the phase diagram of hard-core bosons in one dimension with infinitesize matrix product state techniques to exemplify the impact of the confinement-induced modification of dipole-dipole interactions on many-body physics. Finally, in Sec. 5, we conclude. In Appendix A we discuss numerical methods for evaluating matrix elements of nonlocal potentials in a basis of Wannier functions, and in Appendix B, we provide an explicit evaluation of the dipole-dipole interaction in a cylindrically symmetric harmonic trap for a comparison with the results in an optical lattice.

Optical Lattices, Wannier functions, and Hubbard models
The theory of quantum mechanical objects in a continuous periodic potential is well established [48]. Here, we present a review of the basic facts in order to set notation, and also provide some explicit computations for the optical lattice potential Eq. (1) which do not appear elsewhere in the literature, to the best of our knowledge.
We consider that our system is subject to the separable simple cubic potential where a is the lattice spacing. The typical energy scale derived from a is the recoil energy E R = 2 π 2 /2ma 2 . A potential of the form Eq. (1) is produced for ultracold gases by an optical lattice consisting of three sets of counter-propagating laser beams. The interaction of an optical potential which is far detuned from any atomic or molecular resonances may be described by the AC Stark shift [49] where E opt (r) is the optical field andα (ω opt ) is the dynamical polarizability tensor of the object evaluated at the optical frequency ω opt . The quantities V x , V y , and V z in Eq. (1), which we will refer to as the lattice heights along the x, y, and z directions, may be tuned by increasing the intensity of the optical field. Throughout this paper, we will use the notationṼ ≡ V /E R to denote the dimensionless ratio of a lattice height to the recoil energy. The assumption of a separable potential such as Eq. (1) applies when the polarizability tensor is a scalar, as occurs naturally for alkali atoms. When the dynamical polarizability contains non-scalar components, internal states of, e.g., a rotating molecule can be coupled together for certain polarizations of the optical potential [50]. The theory given in this paper can be extended to this case, but the analysis is more complex. For clarity of exposition, we will focus on the separable case, Eq. (1).
For particles subject to a periodic lattice potential, the energy eigenfunctions are Bloch functions ψ nk (r) characterized by a quasimomentum index k in the first Brillouin zone (BZ) and a band index n. Bloch functions are the simultaneous eigenfunctions of the single-particle Hamiltonian and the lattice translation operators, and so represent the analogs of plane waves in free space when the translational symmetry is a discrete, rather than continuous, group. As such, Bloch functions are delocalized objects, and so are often not an appropriate basis for expanding a many-body Hamiltonian with strong local interactions. A more appropriate basis for describing strong interactions in lattices is provided by Wannier functions, which are the quasimomentum Fourier transforms of the Bloch functions, Here L is the total number of unit cells in a lattice with periodic boundary conditions and r i denotes the coordinate of lattice site i. For simplicity, we will restrict our attention to the lowest band in most cases, and drop the band index n. The extension of our results to multi-band situations [51] is straightforward. Additionally, our methods also extend readily to other bases, for example localized bases which take into account strong contact interactions [52,53,54,55].

Properties of Wannier functions
It was shown in a seminal work by Kohn [24] that the phases on the Bloch functions can be chosen such that the Wannier functions are exponentially decaying away from their centers for one-dimensional centro-symmetric lattice potentials in a sense to be discussed in the next paragraph. The Wannier functions with this choice of phases are called maximally localized Wannier functions, and are used throughout this paper. The exponential decay of Wannier functions has been extended to general one-dimensional lattice potentials [56], non-degenerate bands in arbitrary dimensions [57], and general two-and three-dimensional insulators with vanishing Chern number [58]. Hence, the exponential localization of Wannier functions, which will play an important role in our results, is a general property which does not require fine tuning or a specific lattice structure.
Because the potential we consider, Eq. (1), is separable, the Bloch solutions of the 3D single-particle Schrödinger equation are products of the Bloch solutions of the 1D Hamiltonian The linearity of the transformation Eq.
Thus, we can use the results of Kohn [24] to discuss the properties of the 1D Wannier functions w in (x). In particular, the Wannier functions decay exponentially as |w 0n (x) | ∼ exp(−h n x) in the sense that The parameter h n is the distance to the nearest branch point from the real axis in the complex quasimomentum plane. The decay length h n can be determined numerically by locating crossings in the band structure computed with quasimomentum k = k n + iκ, k n = π(1 − (−1) n+1 )/2a and κ real. At the point k where bands n and n + 1 cross, h n = κ. The results of this analysis for the lowest two bands of the 1D potential Eq. (4) are shown in Fig. 1(a). He and Vanderbilt [59] pointed out that Eq. (6) is consistent with an exponential decay multiplied by an algebraic factor. They then demonstrated that the precise asymptotic behavior of maximally localized, orthogonal Wannier functions in 1D is for a lattice consisting of periodically repeating Gaussian wells. By performing numerical fits with the obtained h n from Fig. 1(a), we find that this same algebraic factor appears in the Wannier functions of the potential Eq. (4), and the power of the algebraic correction does not appear to depend on the lattice height or the band index. The resulting fits to the Wannier functions, together with the numerically computed Wannier functions, are shown in Fig. 1(b). Fig. 1 provides a rigorous definition of the decay of the Wannier functions we consider in this paper as a function of the height of the optical lattice potential.

Approximation of Wannier functions by Harmonic oscillator functions
It is difficult to obtain quantitative analytical predictions from Wannier functions due to their complicated form. Hence, approximations to the true Wannier functions are often made for analytical convenience. The most common such approximation is to replace a single site of the optical lattice by a harmonic well with the same local curvature. We shall call this the harmonic oscillator approximation (HOA). The curvature-matching condition amounts to ν = a/(πṼ 1/4 ν ), where ν is the harmonic oscillator length along Cartesian direction ν. The 1D ground state wave function, which is taken to approximate the 1D lowest band Wannier function, is There are several important differences between the HOA Eq. (8) and the true Wannier function. First, Eq. (8) is everywhere positive, while Wannier functions, even for the lowest band, always have nodes in order to maintain orthogonality between lattice sites. Thus, the HOA poorly captures quantities like tunneling, which involve overlaps of derivatives of Wannier functions. Also, Eq. (8) decays much more rapidly than the true Wannier functions, as a Gaussian rather than an exponential. Hence, the HOA will consistently underestimate the overlap of Wannier functions at different sites. Third, the HOA Eq. (8) is more peaked around its center value than the true Wannier function, and so quantities computed on-site are overestimated by the HOA. Finally, the true Wannier functions have all of the lattice symmetries, while the HOA displays the symmetries of an ellipsoid for the simple cubic lattice. The difference in symmetries is especially important when discussing interaction-induced diagonal tunneling [27,28]. We can make the comparison between Wannier functions and the HOA more quantitative by considering the moments of their associated single-particle probability  distributions. We define the p th one-dimensional moments as where ψ = w and g stand for Wannier and Gaussian, respectively. The second moment gives an estimate for the width of the distribution. It is convenient to define the kurtosis which is a measure of the peakedness of a distribution as well as the heaviness of its tails. We show the second and fourth moments of the true Wannier functions and the HOA in Fig. 2. The width of the lowest band Wannier function is well captured by the HOA for deep lattices withṼ 20, although the HOA always has a smaller second moment than the true Wannier function. The harmonic oscillator always has a vanishing kurtosis, κ g = 0, as is known for Gaussian distributions. In contrast, the Wannier functions have an always positive kurtosis which is very sizable for shallow lattices and approaches zero for deeper lattices. For example, we observed κ w = 10 for V = 2 and κ w = 0.12 forṼ = 35. This difference in the kurtosis quantifies the essential difference between the decay of the Wannier functions and the HOA.

Derivation of Hubbard models for dipolar particles
With the identification of the lowest band Wannier functions as the appropriate singleparticle basis for describing strongly interacting particles in a lattice, we derive a manybody lattice model using the well-known procedure [60,61] of expanding the field operatorψ(r) in the basis of Wannier functions and substituting this expansion into the second-quantized expression for the interaction Hamiltonian, Dipole-dipole interactions in optical lattices do not follow an inverse cube power law 8 where V int (r) is the two-particle interaction potential. The expansion of Eq. (10) in Wannier functions yieldŝ whereâ i destroys a particle in a lowest band Wannier state centered at site i and Here, f ii (r) ≡ w i (r)w i (r) is a product of Wannier functions. The matrix elements U i 1 i 2 i 2 i 1 , which we will call Hubbard parameters, describe the interactions between particles localized in lowest band Wannier states. The interaction we consider is the dipole-dipole interaction, commonly given as whered is the dipole operator, r is the relative position of the interacting particles, r = |r|, and e r is a unit vector in the direction of r. It is convenient to instead recast the dipole-dipole potential as the contraction of two spherical tensors as [62] V where C is an unnormalized spherical harmonic in the spherical coordinates of r and d ⊗d is the irreducible tensor product of two dipole operators with j 1 m 2 j 2 m 2 |jm a Clebsch-Gordan coefficient. The operatorŝ are the spherical decompositions of the vector operatord. The interaction Hamiltonian may be written in terms of coupling constants U DD,q and geometrical factors G DD,q aŝ where All of the information about the size of the dipole moment and the lattice spacing are contained in the coupling constant Eq. (18). On the other hand, the geometrical integral Eq. (19) contains the information about the effects of lattice confinement on the effective interaction through the Wannier function products f ii (r). It is these geometrical factors which we are interested in studying in the present article. While our methods can be applied to all components q of the dipole-dipole interaction, for simplicity we focus on the terms in Eq. (17) with q = 0. These terms do not change the projection of the total angular momentum along a spacefixed quantization axis. The q = 0 components are the only terms relevant for ultracold 1 Σ molecules in an optical lattice and oriented in a strong DC electric field, as states with dipole-allowed transitions are separated by energy splittings large compared to the characteristic dipole-dipole interaction energy [63,45]. Also, these are the most relevant processes for magnetic dipoles in an optical lattice with a strong magnetic field to prevent dipolar relaxation and spontaneous demagnetization [64,65,66], or a quantum simulation of magnetic dipoles using symmetric top molecules in a strong electric field [67]. In this case, we have G DD,0 where θ is the polar angle of the relative coordinate (r − r ). The commonly used deltafunction approximation (DFA) [42,43,44,45,46] replaces f ii (r) → δ (r − r i ) δ ii , such that G DD,0 where θ i 1 i 2 is the polar angle of the relative vector between lattice sites i 1 and i 2 with positions r i 1 and r i 2 , respectively. The terms in Eq. (17) may be viewed as scattering processes in which particles in the lowest band at lattice sites i 1 and i 2 move to lattice sites i 1 and i 2 , respectively, in the course of an interaction. The largest magnitude terms in Eq. (17) are those in which i 1 = i 1 and i 2 = i 2 , which we will call direct interactions using the scattering process analogy. The direct terms arise as density-density interactions in the interaction Hamiltonian, proportional ton i 1n i 2 . Another class of terms which are proportional tô n i 1n i 2 are the exchange interactions, in which i 1 = i 2 and i 2 = i 1 with i 1 = i 2 . All other processes involve changing the position of a Wannier state from its initial position. Examples of such processes are assisted tunneling [46] and pair hopping [27,68], both of which can introduce new quantum phases in shallow lattices [68]. However, assisted tunneling, pair hopping, and exchange terms are all suppressed by factors related to the exponential decay of the Wannier functions, see Fig. 1, and the behavior of these Hubbard parameters with lattice confinement is qualitatively similar. Because the focus of this work is on the qualitative behavior of the Hubbard parameters with the lattice confinement, we will restrict our attention to the direct and exchange terms with the understanding that the behavior of assisted tunneling and pair hopping is similar to that of the exchange term. It should be noted that our methods apply to any Hubbard parameter. Hence, keeping only the direct and exchange interactions and taking into account the statistics of the particles, the expansion Eq. (11) may be written for the dipole-dipole interaction aŝ where we have defined U ≡ U DD,0 and The plus (minus) sign on the exchange term in Eq. (25) refers to bosons (fermions). Results for the effective interactions, Eq. (25), computed using the methods of Appendix A, are presented in the next section.

Effective dipole-dipole interactions
With respect to the length scale of the lattice constant a, the short range physics is given by the interactions I 0 which occur within a unit cell. A point of comparison for the on-site dipolar interaction I 0 is provided by the dimensionless interaction energy I ho 0 = Ĥ DD /U of two bosonic particles in the ground state of a harmonic trap, where the harmonic trap is chosen to match the local curvature of a lattice site as detailed in Sec. 2.2. In the case where the oscillator lengths in the xy plane are equal, x = y = ⊥ , we have that where¯ ≡ ( z 2 ⊥ ) 1/3 /a is the geometric mean oscillator length in units of a, α ≡ z / ⊥ measures the confinement asymmetry, and β ≡ α(1 − α 2 ) −1/2 . Notably, the interaction energy I ho 0 vanishes for isotropic confinement, α = 1. For α < 1, corresponding to stronger confinement along the quantization axis, contributions from θ > arccos 1 3 , where the dipole-dipole potential is negative, are suppressed. Hence, I ho 0 is positive for α < 1. In contrast, for α > 1 where confinement is weakest along the quantization axis, I ho 0 is negative. These qualitative features are shared by the true effective interaction using Wannier functions, as is shown in Appendix B. Of particular interest is that the on-site dipole-dipole interaction vanishes when the lattice heights are equal, V ≡Ṽ x =Ṽ y =Ṽ z . Also, as noted in Sec. 2.2, the HOA energy I ho 0 is always greater than the on-site energy computed using Wannier functions.
For the moderate-to long-range interactions I i ,i , let us define I j ≡ I i+j,i with j along the x direction. The parameters I j correspond to effective dipole-dipole interactions separated by a distance of j lattice spacings along a principal axis. The DFA predicts that the exchange contribution vanishes and the direct term contributes a factor of j −3 such that I j = j −3 . Hence, there are two possible sources of deviation from the DFA. The first source is a non-vanishing contribution from the exchange term. The second source is a deviation of the direct term from j −3 . The exchange term is nonvanishing due to overlap between Wannier functions on different sites, and so decreases with increasing lattice height as the exponential localization factor h 1 increases, see Fig. 1. We can expect that the exchange processes will be negligible in magnitude even at the nearest-neighbor distance when 2h 1 a > 1, as then the Wannier functions are welllocalized within a unit cell. From Fig. 1, this gives an estimation ofṼ ≈ 5. In contrast to the modification of the effective interaction by the exchange term, modification of the effective interaction due to the direct term relies essentially on asymmetry in threedimensional confinement. This can be understood in analogy with the effective shortrange interaction discussed in the last paragraph, and is a consequence of the fact that the dipole-dipole interaction has relative d-wave symmetry.
To demonstrate the effects of asymmetric lattice confinement, we evaluate the parameters I j for two quasi-low dimensional scenarios. In the quasi-2D scenario, we take the z direction to be tightly confined with a lattice heightṼ z = 45 and the lattice heights along the x and y directions to be equal,Ṽ ⊥ ≡Ṽ x =Ṽ y . In the quasi-1D scenario we take both the z and y lattices to be tightly confining withṼ z =Ṽ y = 45, and call the x lattice heightṼ ⊥ . These quasi-dimensional reductions are quite moderate and are commonplace in current experiments [69,70]. The confinement-induced modifications of effective dipole-dipole interactions with respect to the DFA in the quasi-1D and quasi-2D scenarios are shown in Fig. 3. Here, the modifications to the direct term are shown as (G 0110 − 1), which is the deviation from the DFA prediction of 1 at the nearest-neighbor distance. The exchange term is identically zero in the DFA, and so any non-vanishing value for the exchange term represents a deviation from the DFA. The red curves are the numerical computation with Wannier functions, and the blue curves use the HOA. We note that the exchange contributions, shown with dashed lines, are drastically underestimated by the HOA, as expected. The modifications of the direct term, shown with solid lines, display qualitatively similar features between the true solution and the HOA. In particular, both show vanishing modifications of the direct term when there is no asymmetry in confinement,Ṽ ⊥ = 45. However, the HOA underestimates the modifications of the direct term, quite significantly for shallow quasi-low-dimensional confinement.
In Fig. 3, we showed the behavior of the confinement-induced modifications of the effective dipole-dipole interactions with confinement asymmetry at the nearest-neighbor distance. These modifications also depend on the distance between sites in the lattice, and this dependence has a length scale set by the confinement. In order to get a complete picture of the dependence of confinement effects on both distance between sites and confinement asymmetry we fit the numerically obtained data for I j to the form Quasi-2D confinementṼ ⊥  27) is exact. We also stress that the exponential constant b e has no a priori relation to the exponential decay constant h n of the Wannier functions. Across a wide range of quasi-low dimensional confinement, a c ∼ 0.2a, and so the moderate range over which confinement modifies interactions is a few lattice sites. It should also be noted that the confinementinduced modifications of the effective interactions will also have a nontrivial angular dependence due to the fact that Wannier functions are not spherically symmetric, but rather have the symmetries of the lattice. We leave investigations of the nontrivial angular dependence for future work.

Many-body physics
To illustrate the implications of our findings for many-body physics, we study a model of quasi-1D hard-core bosons with long-range dipolar interactions [63] Here, the nearest-neighbor tunneling amplitude is t, i, j denotes nearest-neighbor pairs i and j, and µ is the chemical potential. We compute the phase diagram of Eq. (28) for two realizations of I j using the infinite size variational ground state search algorithm for matrix product states (iMPS) [71]. In one realization, we use the DFA result, I j = j −3 . In the second realization, I j is the effective interaction in an optical lattice computed with Wannier functions, Eq. (25). We take the lattice heights to bẽ V z =Ṽ y = 45,Ṽ x = 6, which fixes the tunneling energy t and ensures that assisted tunneling terms are small. The large confinement asymmetry enforces the hard-core constraint through a large on-site interaction, see Appendix B. The coefficient U can be tuned for ultracold polar molecules by an applied DC electric field. One can determine the strength of the interaction U = U DD,0 knowing only the expected dipole moment and the lattice constant from single-particle physics, and so scaling the phase diagram to U is appropriate. The iMPS method assumes that the many-body ground state has translational invariance under shifts by q sites, and represents the wavefunction of the q-site unit cell as a matrix product state (MPS) [47] with entanglement cutoff χ. One minimizes the energy functional of the unit cell variationally in the parameters of the matrix product, using a sweeping procedure across sites in the unit cell reminiscent of the densitymatrix renormalization group (DMRG) procedure for finite lattices [72]. Alternating with unit cell minimization is an updating of the effective environment of the unit cell using the most current parameters of the unit cell. The procedure is halted when the 2-norm distance between the unit cell in successive optimizations drops below a desired tolerance. Long-range interactions are facilitated within iMPS by using matrix product operators and fitting the interaction I j to a sum of exponentials [73,74,75]. The advantage of the iMPS method for the problem at hand is its ability to handle strong long-range interactions without the significant boundary effects that plague DMRG-type procedures on finite lattices. All results given here have been checked for convergence in the entanglement cutoff χ and the unit cell length q.
The phase diagrams of Eq. (28) for the DFA and lattice-modified dipolar interactions are shown in Fig. 5. The bordered areas represent predicted regions of crystalline phase (CP) with density ρ = 1/2 and a non-vanishing single-particle gap. The remainder of the plots are a gapless superfluid (SF) phase. In the CP, the density correlation function N (r) ≡ (n 0 − ρ)(n r − ρ) behaves as N (r) → cst.(−1) r as r → ∞ and the single-particle density matrix A(r) = â † 0â r is exponentially decaying with r. In panel (a), the area bounded with dashed lines represents the region of CP computed using the DFA in the parameters µ/U and µ/U . The green area bordered with solid lines represents the CP boundaries computed with the actual lattice dipolar interaction. The region predicted by the DFA is shifted both in chemical potential and tunneling by approximately I 1 − 1 ≈ 36%. However, the difference between the DFA and the true solution is not a simple rescaling of the axes. This is shown in panel (b), which rescales the chemical potential and tunneling to the nearest-neighbor interaction energy, µ/(U I 1 ) and t/(U I 1 ). Similar changes in the phase diagram will occur for the 2D case, with shifts of about 48% for the parameters of Ref. [43]. For soft-core particles which also possess local, isotropic interactions, the modified dipolar interaction will be relevant to the convexity of the interaction potential and hence to the formation of supersolid phases [76].

Conclusions
In conclusion, we have shown that the dipole-dipole interaction is strongly modified by imperfect localization of particles in an unequally confined lattice, as found in many experiments in ultracold quantum gases. In particular, we demonstrated that the effective dipole-dipole interactions in optical lattices, given by the matrix elements of the dipole-dipole interaction in a basis of lowest band Wannier functions, are significantly modified from the free-space form of the dipole-dipole interaction. The modifications of the dipole-dipole interaction arise both from interactions in the exchange channel, controlled by the overlap between Wannier functions at different sites, and interactions in the direct channel. The latter interactions depend crucially on the asymmetry in lattice confinement due to the d-wave symmetry of the dipole-dipole interaction. We compared our results to results obtained by approximating Wannier functions as localized Gaussians. Interactions in the direct channel are qualitatively reproduced by this approximation, though estimates of the magnitude can be off by significant factors. Interactions in the exchange channel are quantitatively very poor in the exchange channel, underestimating interactions by orders of magnitude. Based on numerical simulations, we put forward a simple characterization of the modified interaction as being exponential at moderate separations of a few lattice sites and power-law for large separations. Using iMPS simulations, we showed that the modified interaction can significantly alter the predictions of many-body systems, including the determination of phase diagrams.
We acknowledge useful discussions with Mahadevan Ganesh, Kenji Maeda, and Zhigang Wu. This research was supported in part by the National Science Foundation under Grants PHY-1207881 and NSF PHY11-25915, by AFOSR grant number FA9550-11-1-0224, the Heidelberg center for Quantum Dynamics, and by the Alexander von Humboldt Foundation. We also acknowledge the Golden Energy Computing Organization at the Colorado School of Mines for the use of resources acquired with financial assistance from the National Science Foundation and the National Renewable Energy Laboratories. We would like to thank the KITP for hospitality.

Appendix A. Numerical procedures to compute Hubbard parameters
In this section, we discuss numerical techniques for computing geometrical integrals such as Eq. (19) which determine the effective interaction between particles localized in lowest band Wannier states. Such integrals may be written in the general form where I (r − r ) is the dimensionless two-particle potential. The first method employs the convolution theorem to express the integration over the primed coordinates in Eq. (A.1) as a function of the unprimed coordinates with high accuracy. The remaining three-dimensional integral is then performed with standard numerical quadrature. We will call this method the real-space method. In the second method, we compute the interaction matrix elements of the dipole-dipole interaction potential in the basis of lowest band Bloch functions. The matrix elements in the Wannier basis are then obtained by quasimomentum Fourier transform. We will refer to this latter method as the Bloch expansion method. Both methods exhibit steep scaling with the linear domain size L which prohibits studies of very large systems. In order to ensure wellconverged results, we restricted our analysis to a separation of at most 7 lattice sites, see Sec. 3.

Appendix A.1. Real-space method
In the real-space method, we begin with Eq. (A.1) and apply the convolution theorem to find where F k [g(r)] denotes the Fourier transform of the function g(r) as a function of k and likewise for the inverse transform F −1 r {•}. For example, the Fourier transform of the spatial part of the dipole-dipole potential is F k [C (2) q (r)/r 3 ] = −4πC (2) q (k)/3, and the Fourier transform of the delta-function potential is a constant. Hence, the evaluation of the Hubbard parameter may be computed by three-dimensional Fourier transforms followed by a three-dimensional integration in real space rather than by a six-dimensional real space integral. To perform these procedures numerically, we consider each Cartesian dimension to be a symmetric finite interval S = [−L/2, L/2] with periodic boundary conditions, and discretize each interval with n g grid points. The grid spacing in the discrete Fourier conjugate domain is 2π/L and the extent of the domain in Fourier space is controlled by πn g /L, the inverse real space step size. The transformation from a function to its discrete Fourier conjugate is performed by the fast Fourier transform (FFT) algorithm in O(n 3 g log n g ) time [77]. Because the Wannier functions on a finite domain are periodic and band-limited, their discrete and continuous Fourier transforms The Fourier transform of the dipole-dipole potential is ill-defined at k = 0, but our results do not depend on this value. are related by a scaling constant provided we sample the entire domain at a frequency of at least twice their largest frequency component [77]. As is known for spectral methods, convergence of the Fourier space calculation in Eq. (A.2) is exponential in L provided that n g /L is large enough to capture the full support of the function in Fourier space. Defining g to be the support of the lowest band Wannier function in the discrete Fourier space, the choice n g = 2gL + 1 ensures that the function is fully captured in Fourier space. The support of a Wannier function in Fourier space can be determined by using Parseval's theorem on finite Fourier subintervals to determine that the norm is unity to a desired tolerance. For typical g ∼ 5 − 7, the real space integration in Eq. (A.2) is of acceptable precision using a high-order Simpson integrator [78].
There are two dominant sources of error in the real space method. The first error is due to the discretization of the real space domain and the associated discretization error of the numerical quadrature. This error may be controlled by increasing g. The second source of error is spurious interactions due to periodic boundary conditions. While these interactions vanish as the domain becomes infinite, convergence may be slow due to, e.g. ,the power-law decay of the dipole-dipole interaction at long range. Hence, we have instead used a finite-size scaling analysis to extrapolate our results to the limit of an infinite lattice. The main limitation on the system sizes that we can reach using the real space method is a memory requirement which scales as O (g 3 L 3 ) due to the non-separability of the dipole-dipole potential.

Appendix A.2. Bloch expansion method
In the Bloch expansion method, we study the matrix elements of the interaction potential in the basis of lowest band Bloch functions ψ q (r), The geometric integral Eq. (A.1) is then given by . Using the fact that the interaction potential depends only on the relative coordinate, we find where R is any Bravais lattice vector and Q ≡ q 1 + q 2 − q 1 − q 2 modulo 2π. This implies that there are only N 3 independent non-vanishing matrix elements in Eq. (A.3) for a lattice with N unit cells, as opposed to the N 4 possible configurations of the four quasimomenta. As in the main text, we will study the case of the simple cubic lattice in which the Bloch functions separate along principal axes, dr ν [ψ q 1ν (r ν ) ψ q 2ν (r ν )] I (r − r ) ψ q 1ν (r ν ) ψ q 2ν (r ν ) . (A.7) Changing integration variables to 2ξ ν ≡ r ν + r ν , 2η ν ≡ r ν − r ν with Jacobian 2 along each Cartesian direction and expanding the one-dimensional Bloch functions as ψ q (x) = lim Here, is a finite Fourier cutoff used in numerics and we have defined The integrals over ξ ν are L−fν ην fν ην dξ ν e 2πitν ξν = Lδ tν ,0 − sin (2πη ν t ν ) πt ν . (A.13) We will keep only the term proportional to the Kronecker delta, as this is the dominant contribution for large lattices. This approximation becomes exact in the limit of an infinite lattice, L → ∞, as has been shown for the delta-function potential in Ref. [28]. For simplicity, we now also require that the interaction potential is invariant under inversion by any Cartesian coordinate. This is true for the q = 0 component of the dipole-dipole interaction, as well as for the delta function potential. However, the method directly extends to more general interactions. With these two conditions, we find in the limit as L → ∞, V q 2 q 2 q 1 q 1 = 1 L 3 ν=x,y,z p 1ν p 2ν p 1ν p 2ν δ p 1ν − p 1ν + p 2ν − p 2ν + q 1ν − q 1ν + q 2ν − q 2ν 2π × c p 1 q 1 c p 2 q 2 c where F k [I] is again the Fourier transform of I as a function of k and ∆ = (∆ x , ∆ y , ∆ z ). The Bloch expansion method has the advantage of not introducing any discretization error, and also requires significantly less memory when the lattice is separable. However, the computational scaling of this method is O (L 9 ) as opposed to O (L 3 ) for the real space method. The Bloch expansion method suffers from spurious interactions due to periodic boundary conditions, just as in the real space method. We find that both methods agree when we extrapolate to the limit L → ∞, and allow us to put a conservative bound of 1% on our estimated error. Additionally, we have benchmarked both methods with the case of a delta-function potential, where both analytical and numerically exact results are available. Effective on-site interaction Planar confinementṼ ⊥ Figure B1. The on-site dipolar interactions in a quasi-2D lattice, I 0 (red solid line), and the dipolar interaction energy of a cylindrically symmetric harmonic oscillator with the same local curvature as the lattice site, I ho 0 (green dashed line), have similar qualitative behavior with respect to confinement asymmetry. In particular, both vanish as the confinement becomes isotropic,Ṽ ⊥ =Ṽ z = 45, and are significant for large confinement asymmetry.
oscillator is always larger than the corresponding energy in an optical lattice, often by 10%. Additionally, the on-site interaction is large for small to moderate V , enforcing the hard-core constraint used in the many-body study in the main text.