Beyond-mean-field corrections for dipolar bosons in an optical lattice

TL;DR

This paper investigates how optical lattices influence beyond-mean-field corrections in dipolar bosons, revealing anisotropic effects and orientation-dependent modifications crucial for understanding ultracold dipolar gases.

Contribution

It provides a detailed analysis of how optical lattices alter beyond-mean-field corrections, including anisotropic effective mass and orientation-dependent effects.

Findings

Deep in the superfluid regime, the equation of state is described by an anisotropic effective mass.

Deep lattices introduce anomalous density dependence terms not present in free space.

The orientation of dipoles relative to the lattice significantly affects the corrections.

Abstract

Recent experiments with ultracold lanthanide atoms which are characterized by a large magnetic moment have revealed the crucial importance of beyond-mean-field corrections in understanding the dynamics of the gas. We study how the presence of an external optical lattice modifies the structure of the corrections. We find that deep in the superfluid regime the equation of state is well described by introducing an anisotropic effective mass. However, for a deep lattice we find terms with anomalous density dependence which do not arise in free space. For a one-dimensional lattice, the relative orientation of the dipole axis with respect to the lattice plays a crucial role and the beyond-mean-field corrections can be either enhanced or suppressed.

Equations

$$H$$

$$+\frac{1}{2}\int d^{3}r\int d^3{r}^{\prime }{\Psi }^{†}(\mathbf{r}){\Psi }^{†}({\mathbf{r}}^{\prime })V(\mathbf{r}-{\mathbf{r}}^{\prime })\Psi ({\mathbf{r}}^{\prime })\Psi (\mathbf{r}).$$

$$V(\mathbf{r})=g_0\delta (\mathbf{r})+\frac{C_{dd}}{4\pi }\frac{1-3{\cos }^2\theta }{r^3},$$

$$\Psi (\mathbf{r})=\sum _{i}w(\mathbf{r}-{\mathbf{r}}_i)a_i$$

$$a_i=\frac{1}{\sqrt{N_L}}\sum _{\mathbf{k}\in K}{a}_{\mathbf{k}}{e}^{-i\mathbf{k}{\mathbf{r}}_i},$$

$$H=\sum _{\mathbf{k}}{ϵ}_{\mathbf{k}}{a}_{\mathbf{k}}^{†}{a}_{\mathbf{k}}+\frac{1}{2}\sum _{i,j,l,m}{V}_{ijlm}{a}_i^{†}{a}_j^{†}{a}_l{a}_m.$$

$${ϵ}_{\mathbf{k}}=-2t\left(\sum _{i=x,y,z}\cos (k_{i}l)-3\right),$$

$$V_{ijlm}$$

$$\times w({\mathbf{r}}^{\prime }-{\mathbf{r}}_l)w(\mathbf{r}-{\mathbf{r}}_m).$$

$$V_{ij}^{\text{contact}}={\delta }_{ij}{g}_0\int d^{3}r|w(\mathbf{r}){|}^4\equiv g{\delta }_{ij}.$$

$$g=g_0\frac{(2\pi {)}^{3/2}}{l^3}{s}^{3/4}.$$

$$w^{\ast }(\mathbf{r}-{\mathbf{r}}_i)w(\mathbf{r}-{\mathbf{r}}_m)\approx {\delta }_{im}\delta (\mathbf{r}-{\mathbf{r}}_i).$$

$$V_{ij}=g{\delta }_{ij}+\frac{C_{dd}}{4\pi }\frac{1-3{\cos }^2{\theta }_{ij}}{|{\mathbf{r}}_i-{\mathbf{r}}_j{|}^3},$$

$$H=\sum _{\mathbf{k}}{ϵ}_{\mathbf{k}}{a}_{\mathbf{k}}^{†}{a}_{\mathbf{k}}+\frac{1}{2}\sum _{i,j}{V}_{ij}{n}_i{n}_j.$$

$$V(\mathbf{k})\equiv V_{\mathbf{k}}=g\left(1+{\epsilon }_{dd}\frac{3}{4\pi }\sum _j{e}^{i\mathbf{k}{\mathbf{r}}_j}\frac{1-3{\cos }^2{\theta }_j}{|\mathbf{j}{|}^3}\right),$$

$${\epsilon }_{dd}=\frac{{\epsilon }_{dd}^0}{l^3\int d^{3}r|w(\mathbf{r}){|}^4}=\frac{{\epsilon }_{dd}^0}{(2\pi {)}^{3/2}{s}^{3/4}}.$$

$$\frac{E^{(2)}}{V}-\frac{1}{2}n{\mu }^{(2)}=\frac{1}{2}{\int }_K\frac{d^{3}k}{(2\pi {)}^3}\frac{({ϵ}_{\mathbf{k}}+n{V}_{\mathbf{k}}-E_{\mathbf{k}})({ϵ}_{\mathbf{k}}-E_{\mathbf{k}})}{2E_{\mathbf{k}}},$$

$$\frac{1}{2}ng\frac{1}{l^3{\alpha }^{3/2}}{\int }_{-\pi \sqrt{\alpha }}^{\pi \sqrt{\alpha }}\frac{d^{3}u}{(2\pi {)}^3}\frac{({ϵ}_{\mathbf{u}}+1-E_{\mathbf{u}})({ϵ}_{\mathbf{u}}-E_{\mathbf{u}})}{2E_{\mathbf{u}}}$$

$$\frac{ng}{l^3{\alpha }^{3/2}}=\frac{(ng{)}^{5/2}(2m^{\ast }{)}^{3/2}}{{\hslash }^3}.$$

$$I(\alpha )=\frac{1}{2}{\int }_{-\pi /\sqrt{\alpha }}^{\pi /\sqrt{\alpha }}\frac{d^{3}u}{(2\pi {)}^3}\frac{({ϵ}_{\mathbf{u}}+1-E_{\mathbf{u}})({ϵ}_{\mathbf{u}}-E_{\mathbf{u}})}{2E_{\mathbf{u}}}.$$

$$\frac{E^{(2)}}{V}-\frac{1}{2}n{\mu }^{(2)}=\frac{(ng{)}^{5/2}(2m^{\ast }{)}^{3/2}}{{\hslash }^3}\left(-\frac{1}{15\sqrt{2}{\pi }^2}\right).$$

$$\frac{E^{(2)}}{V}=\frac{8}{15{\pi }^2}\frac{(m^{\ast }{)}^{3/2}(ng{)}^{5/2}}{{\hslash }^3}$$

$$V_{dd}^{\text{free}}(\mathbf{k})=g{\epsilon }_{dd}\left(3\frac{k_z^2}{k_x^2+k_y^2+k_z^2}-1\right),$$

$$I(\alpha \to \infty )=-\frac{1}{15\sqrt{2}{\pi }^2}F({\epsilon }_{dd})$$

$$e_0^{(2)}+\frac{1}{2}\alpha \frac{d{e}_0^{(2)}}{d\alpha }$$

$$=-\frac{F({\epsilon }_{dd})}{15\sqrt{2}{\pi }^2}{\alpha }^{-5/2}\left[1+f(\alpha )\right]\equiv P(\alpha ).$$

$$e_0^{(2)}=-\frac{1}{{\alpha }^{-2}}{\int }_{\alpha }^{\infty }dy2P(y)y.$$

$$\frac{E_0^{(2)}}{V}=-n^2{\int }_0^{n}dy\frac{2G(y)}{y^3}.$$

$$P(\alpha )=P_0(\alpha )+\Delta P(\alpha )=-\frac{F({\epsilon }_{dd})}{15\sqrt{2}{\pi }^2}{\alpha }^{-5/2}\left[1+f(\alpha )\right].$$

$$\frac{E_0^{(2)}}{V}=\frac{4F({\epsilon }_{dd})}{15\sqrt{2}{\pi }^2}\frac{(2m^{\ast }{)}^{3/2}(ng{)}^{5/2}}{{\hslash }^3}.$$

$$\Delta e_0^{(2)}$$

$$=\frac{4F({\epsilon }_{dd})}{15\sqrt{2}{\pi }^2}\frac{c}{1+2\gamma }{\alpha }^{-5/2-\gamma }.$$

Full text

I Introduction

Dilute gases of weakly interacting ultracold bosons are commonly described in terms of mean-field theory which predicts the formation of a Bose-Einstein condensate (BEC) described by a macroscopic wave function. The presence of long-range interactions can strongly affect the properties of the gas. Dipolar interactions are both long range and anisotropic, which brings in a number of interesting effects and possible applications Lahaye et al. (2009); Baranov et al. (2012). For example, the partly attractive nature of the interaction can lead the gas to collapse Lahaye et al. (2008), which can be experimentally controlled by tuning the strength of the short-range repulsion by means of Feshbach resonances Chin et al. (2010). A confined strongly dipolar gas can exhibit a roton-maxon excitation spectrum Santos et al. (2003); Ronen et al. (2007), which has been experimentally observed Chomaz et al. (2018).

In recent years there has been a tremendous progress in experimental investigations of strongly dipolar systems, enabled by the realization of BEC of lanthanide atoms: erbium and dysprosium Lu et al. (2011); Aikawa et al. (2012). It has then been discovered that when the dipolar gas is close to the stability boundary, it can form a set of stable dense droplets instead of collapsing Kadau et al. (2016). This effect can be explained by effective many-body repulsion induced by the leading beyond-mean-field corrections. These corrections are known as Lee-Huang-Yang (LHY) corrections in the case of contact interactions Lee et al. (1957); Beliaev (1958); Hugenholtz and Pines (1959) and have also been extended to dipolar interactions Lima and Pelster (2012). The physics of these dipolar droplets has been the subject of intense studies Kadau et al. (2016); Ferrier-Barbut et al. (2016); Wächtler and Santos (2016); Baillie et al. (2016); Schmitt et al. (2016); Chomaz et al. (2018). In the case of a single-component contact interacting Bose gas, the beyond-mean-field corrections are rather weak, however, those terms become dominant when the contribution from the mean-field term vanishes as in the case of partially attractive two-species mixtures Petrov (2015); Cabrera et al. (2018) and dipolar gases Ferrier-Barbut et al. (2016). Recently, the behavior of the beyond-mean-field corrections in the case of a dimensional crossover from three to low dimensions has been studied for contact Ilg et al. (2018); Zin et al. (2018) and dipolar Edler et al. (2017) interactions.

Optical lattices are an extremely useful tool for manipulation and control of ultracold gases Jaksch et al. (1998); Lewenstein et al. (2007); Bloch et al. (2008). They allow for creating a perfectly periodic external potential for the atoms with tunable depth and geometry. This enabled the experimental realization of the Bose-Hubbard model and demonstration of the quantum phase transition between the superfluid and Mott insulator phase Greiner et al. (2002). The interplay of dipolar interactions and optical lattice confinement gives rise to a variety of phenomena induced by the long-range nature of the interactions as well as their anisotropy Góral et al. (2002); Büchler et al. (2007); Capogrosso-Sansone et al. (2010); Sowiński et al. (2012); Wall and Carr (2013).

In this paper, we study the influence of an optical lattice on the beyond-mean-field corrections for a dipolar Bose gas. In particular, we investigate whether and how one can control the strength of the quantum fluctuations in a suitable way. We calculate the LHY correction for the case of a deep three-dimensional lattice as well as weak one-dimensional lattice. We find that in general the presence of the lattice enhances the fluctuations, but also introduces a nontrivial density dependence. Moreover, manipulating the relative orientation of the lattice and the dipole axis allows for controlling the strength of the LHY term.

This paper is structured as follows. In Section II, we consider the case of a deep three-dimensional cubic lattice. We compute the LHY correction to the ground state energy of an interacting Bose gas for contact (on-site) interactions, as well as for dipolar density-density interactions neglecting the interaction-induced tunneling effects. We show the emergence of an effective mass in the limit of a large healing length and calculate the corrections arising for larger interaction strength. In Section III, we turn to the case of weak one-dimensional lattice. We show that by manipulating the orientation of the lattice with respect to the dipole orientation axis it is possible to tune the magnitude of the LHY correction.

II Deep optical lattice

II.1 Hamiltonian

We start our considerations with the general many-body Hamiltonian of an interacting Bose gas at zero temperature in second quantized form,

[TABLE]

In this Hamiltonian, $U(\mathbf{r})$ describes an external potential, e.g. a trapping potential or an optical lattice, while $V(\mathbf{r})$ describes the interaction potential (see Fig. 1). For dipolar particles which are polarized along one direction, the interaction potential can be represented in terms of a pseudopotential Yi and You (2001)

[TABLE]

where $\theta$ is the angle between the direction of the polarization of the dipoles and their relative position. The first term of this pseudopotential accounts for the $s$-wave interactions which result from the short-ranged van-der-Waals interactions but contains also the contribution of the dipolar part of the potential. The second part stands for the higher partial waves, which are dominated by the long-range and anisotropic dipolar interaction. The coupling constant $g_{0}$ is related to the $s$-wave scattering length $a_{s}$ via $g_{0}=4\pi\hbar^{2}a_{s}/m$. The dipolar coupling constant $C_{dd}$ is $\mu_{0}\mu^{2}$ for particles having a magnetic dipole moment with $\mu_{0}$ being the permeability of the vacuum, and $d^{2}/\epsilon_{0}$ for particles having an electric dipole moment $d$ with $\epsilon_{0}$ being the permittivity of the vacuum. In free space, the relative strength of the contact and the dipolar interaction is determined by the dimensionless parameter $\varepsilon_{dd}^{0}=a_{dd}/a_{s}$ with the so-called dipolar length $a_{dd}=mC_{dd}/12\pi\hbar^{2}$.

Now, we consider the particles to be trapped in a deep three-dimensional simple cubic optical lattice described by the trapping potential $U(\mathbf{r})=U_{L}\sum_{i=x,y,z}\sin^{2}(q_{L}r_{i})$, where $U_{L}$ is the lattice depth and $q_{L}=2\pi/\lambda$ is the lattice vector with $\lambda$ being the laser wavelength and the lattice period is given by $l=\lambda/2$. The typical energy scale of a particle in a lattice is the recoil energy $E_{R}=\hbar^{2}q_{L}^{2}/2m$ and the strength of the lattice can be characterized by the dimensionless parameter $s=U_{L}/E_{R}$. In the case of a deep lattice, we have $s\gg 1$. In this regime, one can assume that only the lowest Bloch band is occupied (in practice this is a good approximation already for $s\sim 10$ Bloch et al. (2008)) and we can rewrite the field operators $\Psi(\mathbf{r})$ and $\Psi^{\dagger}(\mathbf{r})$ in terms of the Wannier functions $w(\mathbf{r}-\mathbf{r}_{i})$ localized around the lattice site $i$ at position $\mathbf{r}_{i}$:

[TABLE]

with the bosonic annihilation operator at lattice site $i$, $a_{i}$. The annihilation operators in real and quasi-momentum space are connected via a discrete Fourier transform

[TABLE]

where $N_{L}$ denotes the number of lattice sites and the summation over $\mathbf{k}$ is restricted to the first Brillouin zone $K$. Using these transformations, the single-particle term of the Hamiltonian becomes diagonal and the total Hamiltonian can be written as

[TABLE]

The dispersion relation of a particle in a deep lattice is given by

[TABLE]

where the hopping amplitude $t$ is related to the lattice depth $U_{L}$, the recoil energy $E_{R}$, and the lattice spacing $l$ Jaksch et al. (1998). Note that the zero of energy is shifted such that $\epsilon_{0}=0$. In the interaction part of (5), the matrix elements $V_{ijlm}$ are given in terms of the Wannier functions

[TABLE]

For deep lattices, the Wannier functions are well localized and the contribution due to the contact interaction is only significant if $i=j=l=m$, such that we may write

[TABLE]

Approximating the Wannier function at a given lattice site by the ground state wavefunction of an harmonic oscillator, $g$ can be calculated explicitly as Jaksch et al. (1998)

[TABLE]

For the dipolar part, we replace the Wannier functions by $\delta$ functions,

[TABLE]

The matrix elements now only depend on sites $i$ and $j$ and are given by

[TABLE]

with $\theta_{ij}$ being the angle between sites $i$ and $j$. Note that the on-site contribution from the dipolar part is zero for an isotropic confinement at each lattice site. The Hamiltonian (5) then reduces to

[TABLE]

Taking into account the spatial structure of the Wannier states and computing the matrix elements in (7) explicitly gives rise to additional processes such as density-assisted and correlated tunneling Sowiński et al. (2012); Wall and Carr (2013). These processes are strongly suppressed for deep lattices due to the exponential decay of the Wannier functions, but can lead to additional corrections for moderate lattice depths. The role of these terms is discussed in Appendix A.2.

For what follows, we also need to transform the interaction part of (5) into momentum space which requires the discrete Fourier transform of the interaction potential and might be written as

[TABLE]

where the parameter $\varepsilon_{dd}$ is renormalized by the lattice and related to its free space variant by

[TABLE]

The last result was obtained using Eq. (9). The exact form of $V_{\mathbf{k}}$ can be obtained analytically under our approximations (see Appendix A.1) and leads to noticeable differences with respect to the free-space Fourier transform in which $V_{\mathbf{k}}$ only depends on the angle between $\mathbf{k}$ and the direction of the polarization of the dipoles but not on the magnitude of $\mathbf{k}$.

II.2 Beyond-mean-field corrections

In order to calculate the beyond-mean-field energy corrections, we restrict ourselves to the case where the system is in the superfluid phase and the mean-field approach correctly describes its properties. The correction to the mean-field energy is then given by Hugenholtz and Pines (1959)

[TABLE]

where $V=N_{L}l^{3}$ is the volume of the system, $n=N/V$ the density, and $\mu^{(2)}$ denotes the second-order correction to the chemical potential, which is related to the energy by $\mu=dE/dn$. The Bogoliubov dispersion relation $E_{k}$ is given by $E_{\mathbf{k}}=\sqrt{\epsilon_{\mathbf{k}}(\epsilon_{\mathbf{k}}+2nV_{\mathbf{k}})}$ with the non-interacting single-particle dispersion $\epsilon_{\mathbf{k}}$ [6] and the (discrete) Fourier transform of the dipolar potential, $V_{\mathbf{k}}$. We note that as the integration is restricted to the first Brillouin zone, the result is in principle always convergent and no renormalization is needed. We will now study the structure of the correction for different cases, starting for simplicity with the contact interactions.

II.2.1 Contact interaction

For a purely contact interacting Bose gas, we have $V_{\mathbf{k}}=g$ and thus no momentum dependence of the interaction potential. In order to simplify the integral on the right-hand side of (15) and to make it dimensionless, we introduce the effective mass $m^{*}=\hbar^{2}/(2tl^{2})$ and the healing length $\xi^{2}=\hbar^{2}/(2m^{*}ng)$. We further introduce the dimensionless parameter $\alpha=\xi^{2}/l^{2}=t/ng$, which parametrizes the relative strength of the interaction. Staying in the superfluid phase requires $\alpha\gg 1$. Using the substitution $k_{i}l=u_{i}l/\xi=u_{i}/\sqrt{\alpha}$, the integral in (15) reduces to

[TABLE]

with $\epsilon_{\mathbf{u}}=-2\alpha\left(\sum_{i=x,y,z}\cos(u_{i}/\sqrt{\alpha})-3\right)$ and $E_{\mathbf{u}}=\sqrt{\epsilon_{\mathbf{u}}(\epsilon_{\mathbf{u}}+2)}$. The prefactor $ng/(\alpha^{3/2}l^{3})$ in front of the integral can equivalently be written as

[TABLE]

This is, up to a constant numerical factor, exactly the form of the LHY correction in the absence of the optical lattice such that its effects are contained solely in the integral

[TABLE]

In the limit $\alpha\to\infty$, which corresponds to free space, the integral $I(\alpha)$ can be calculated analytically and approaches the constant $I(\alpha\to\infty)=-1/15\sqrt{2}\pi^{2}$ such that

[TABLE]

Solving the differential equation with the initial condition $E^{(2)}(n=0)=0$ ensuring the correct mean-field result, one recovers the well-known free-space LHY term

[TABLE]

with an effective mass accounting for the underlying lattice structure. For finite $\alpha$, one can compute the integral numerically and the result is shown in Fig. 2. As expected, there are corrections to the free-space value at finite $\alpha$ which introduces an additional density dependence to the right-hand side of (15). We discuss this additional density dependence and its influence on the beyond-mean-field energy corrections in Sec. II.2.3.

II.2.2 Dipolar Interaction

We now turn to the case of dipolar interactions. It is interesting to study both the differences which arise with respect to the contact interaction case as well as the effect of the lattice on the dipolar gas. In order to achieve better understanding of the role of the lattice, we first approximate the lattice Fourier transform of the dipolar interaction given in (61) by the free-space Fourier transform

[TABLE]

which only depends on the angle between $\mathbf{k}$ and the direction of polarization of the dipoles (which is assumed to be the $z$ axis). The integral $I(\alpha)$ can only be calculated numerically for finite $\alpha$ and the results for different values of $\varepsilon_{dd}$ are shown in Fig. 2. The asymptotic value for $\alpha\to\infty$ can again be calculated analytically and reads as

[TABLE]

with $F(\varepsilon_{dd})=\frac{1}{2}\int_{-1}^{1}du\,(1+\varepsilon_{dd}(3u^{2}-1))^{5/2}$ accounting for the anisotropic nature of the dipolar interaction (see also Lima and Pelster (2012)). In the absence of dipolar interactions, it reduces to $F(0)=1$. In the case of finite $\alpha$, the integral $I(\alpha)$ leads to corrections qualitatively similar as in the case of contact interaction, while their magnitude increases with increasing $\varepsilon_{dd}$.

The effect of the dipolar interaction is even more enhanced when taking the full lattice Fourier transform of the dipolar interaction as derived in Sec. A.1. For $\alpha\to\infty$, $I(\alpha)$ also approaches the free-space result and is equivalent to taking the free-space Fourier transform (21), whereas the deviations with respect to the case of contact interaction are more prominent for finite $\alpha$. This can be understood by comparing the lattice Fourier transform and the free-space Fourier transform as shown in Fig. 3. In contrast to free space, the lattice Fourier transform of the dipolar potential also depends on the absolute value of the momentum. For large values of $\alpha$, the integral in (15) only gives a contribution near $\mathbf{k}=0$ where there is little difference between both cases, while for decreasing values of $\alpha$, the integral also probes higher momenta where the lattice enhances the effect of the dipolar potential.

II.2.3 General structure of the energy correction

In the previous sections, we have seen that in the presence of a deep lattice, corrections to the usual free-space behavior arise. These corrections lead to a modified scaling in the density which is discussed in this section. In order to do so, we rewrite (15) expressing the energy density $E_{0}^{(2)}/V$ in terms of $t/l^{3}$, that is $E_{0}^{(2)}/V=(t/l^{3})e_{0}^{(2)}$. Substituting $n\to t/\alpha g$ and noting that $d/dn=-(\alpha^{2}t/g)\,d/d\alpha$, we end up with the dimensionless differential equation

[TABLE]

Note that in this expression, all contributions reminiscent of the free-space result are written explicitly and the function $f(\alpha)$ provides all the corrections coming from the lattice and $f(\alpha\to\infty)=0$. The differential equation can formally be solved and using the initial condition $e_{0}^{(2)}(\infty)=0$ such that we obtain the correct mean-field result, the solution reads as

[TABLE]

This is the particular solution of (23), while the homogeneous solution would only affect the mean-field energy in which we are not interested right now. From another point of view, one can also solve (15) with an inhomogeneity $G(n)$ which leads to the solution

[TABLE]

Substituting $\alpha=t/ng$ and expressing the energy density in terms of $t/l^{3}$ leads to the solution (24). For the rest of the discussion, we return to the dimensionless form of the solution, (24). In the presence of the lattice 111The following discussion is essentially more general provided that there exists a dimensionless parameter such that the inhomogeneity has the same scaling as in free space in some limit of this parameter., we can split the inhomogeneous term on the right-hand side of the differential equation into one part containing the free-space result and one part containing the corrections arising from the lattice:

[TABLE]

The first term gives rise to the standard result in free space $e_{0}^{(2),0}=\frac{4}{15\sqrt{2}\pi^{2}}F(\varepsilon_{dd})\alpha^{-5/2}$ which leads to the beyond-mean-field correction

[TABLE]

For the second part, we see from Fig. 2 that we can describe the corrections due to the lattice by a function $f(\alpha)=c\alpha^{-\gamma}$ for $\alpha\gg 1$. The parameters $c$ and $\gamma$ will in general depend on the relative dipolar interaction strength $\varepsilon_{dd}$ and will later be determined by fitting the expression to the numerically calculated $I(\alpha)$. Note that $\gamma>0$ as $f(\alpha\to\infty)=0$. Thus, the energy correction due to the second part reads as

[TABLE]

Finally, the complete beyond-mean-field energy correction reads as

[TABLE]

Since $\gamma>0$, the correction to the beyond-mean-field correction to the ground state energy due to the lattice increases with increasing density. In the limit $t\gg ng$, (29) reduces to the free-space result with a renormalized mass as discussed before. For the present setup of a three-dimensional simple cubic lattice, we determine the coefficients $c$ and $\gamma$ by fitting to the results obtained using numerical integration. For the scaling coefficients, we get $\gamma\approx 1/2$ independent of $\varepsilon_{dd}$ and valid also for contact interactions, whereas the coefficient $c$ ranges from $c\approx 0.3$ for $\varepsilon_{dd}=0$ to $c\approx 0.76$ for $\varepsilon_{dd}=0.7$.

III Weak one-dimensional lattice

Up to now, we have considered a three-dimensional optical lattice. However, an additional intriguing possibility is to restrict the lattice to one dimension and play with the relative orientation of the dipoles and the wave vector of the lattice. In this section, we demonstrate that this generates additional corrections to the usual beyond-mean-field corrections that can be enhanced or diminished depending on the relative orientation between the lattice and the dipoles. The basic assumption of our analysis is that the lattice is weak and can be treated as a perturbation to the free-space case and we do not restrict ourselves to the lowest Bloch band.

III.1 Model

The one-dimensional lattice potential is given as

[TABLE]

where, as before, $q_{L}=\pi/l$ is the lattice vector and $l$ is the lattice period (see also Fig. 4). The orientation of the lattice relative to the dipole axis, which we assume to be parallel to the $z$ axis, can be varied by varying the direction of the lattice vector $\mathbf{q}_{L}$. In the presence of this lattice, the eigenfunctions of the non-interacting problem are given in terms of a product of plane waves (in those directions where no lattice is present) and Bloch functions (in the direction of the lattice),

[TABLE]

In this basis, the single-particle part of the Hamiltonian is diagonal, $H_{0}=\sum_{\mathbf{k}}\epsilon_{\mathbf{k}}a_{\mathbf{k}}^{\dagger}a_{\mathbf{k}}$, where $\epsilon_{\mathbf{k}}=\hbar^{2}\mathbf{k}_{\perp}^{2}/2m+\tilde{\epsilon}_{k}$. The components of $\mathbf{k}$ that are parallel to the lattice are restricted to the first Brillouin zone, $k\in[-\pi/l,\pi/l]$. In general, the lattice dispersion $\tilde{\epsilon}_{k}$ cannot be written down in closed form for arbitrary lattice depth while for weak lattices there are approximations (see Appendix B). As for the deep lattice, we introduce the recoil energy $E_{R}=\hbar^{2}q_{L}^{2}/2m$ and the dimensionless lattice depth $s=U_{L}/E_{R}$ which for a weak lattice is assumed to be much smaller than one, $s\ll 1$.

III.2 Contact interaction

We proceed as in the previous section and start with the contact interaction where the interaction part of the Hamiltonian takes the form

[TABLE]

We now apply the Bogoliubov theory following, e.g., Ref. Abrikosov et al. (2012), which leads to the excitation spectrum

[TABLE]

as well as the beyond-mean-field contributions to the ground-state energy per volume

[TABLE]

Since the lattice is assumed to be weak, the dispersion $\epsilon_{\mathbf{k}}$ will only have small deviations from the dispersion in free space such that we can split the result into a term corresponding to the free space result, the well-known LHY term, and one additional term which includes all the corrections to it. The free space result diverges which can be cured by a proper renormalization of the coupling constant $g$ (see, e.g., Ref. Abrikosov et al. (2012)). We focus here on the remaining corrections to the LHY term due to the lattice which can be expressed as

[TABLE]

where $E_{\mathbf{k}}^{s}$ ($\epsilon_{\mathbf{k}}^{s}$) and $E_{\mathbf{k}}^{0}$ ($\epsilon_{\mathbf{k}}^{0}$) denote the Bogoliubov excitation spectrum (single-particle dispersion) in the presence of the lattice and in free space, respectively.

Similarly to the case of a deep lattice, we introduce a dimensionless quantity

[TABLE]

which compares the healing length $\xi$ with the optical lattice period $l$.

As we are only interested in a weak lattice, we calculate the beyond-mean-field corrections due to the lattice (35) to leading order in $s$, for which analytical results can be obtained (see Appendix B for further details). The general dependence of the beyond-mean-field corrections on the parameter $\beta$ in the case of a contact-interacting Bose gas is shown in Fig. 5 [blue (solid) line] and also includes the leading behavior in the limits $\beta\ll 1$ [green (dash-dotted) line] and $\beta\gg 1$ [orange (dashed) line] which are discussed below.

In the limit of $\beta\to 0$, which corresponds to $l/\xi\gg 1$, the leading behavior of Eq. (97) is given by

[TABLE]

with the usual LHY correction

[TABLE]

The vanishing influence of the lattice in this case is intuitively clear as the limit $l/\xi\gg 1$ indeed should correspond to free space.

In the opposite limit, $\beta\to\infty$ ($l/\xi\ll 1$), the correction instead approaches a constant and reads as

[TABLE]

Combining this with the mean-field energy, the ground state energy per volume is given by

[TABLE]

Both the mean-field term and the beyond-mean-field term are enhanced by the presence of the lattice in the same way in the leading order. The correction to the mean-field term comes from the correction to the $k=0$ mode due to the coupling to modes with $k=\pm q_{L}$ in the presence of the lattice.

Interestingly, this result can also be obtained when assuming that the bosons are in free space but acquire an anisotropic effective mass along the direction of the lattice (say in $z$ direction) with $m_{z}=m/(1-s^{2}/32)=m/\gamma^{2}$ to leading order in $s$. This immediately leads to

[TABLE]

Figure 5 shows that the asymptotic value for $\beta\gg 1$ is already approached for $\beta$ on the order of 10. This suggests that describing the system as a free-space gas with an effective anisotropic mass should be adequate already at values of $\beta$ on the order of 10 which corresponds to $\xi\approx 3l$. Typically, the lattice period is on the order of a few hundred nanometers while the healing length is on the order of micrometers for standard cold-atom experiments Bloch et al. (2008) such that this approximation can be used in typical experiments.

III.3 Dipolar interactions in the $\beta\to\infty$ limit

After having discussed the case of a contact interacting gas, we turn to the the case where the particles interact via dipolar interactions given by the interaction potential

[TABLE]

where $\theta$ is the angle between the direction of the dipole moment, which is assumed to be parallel to the $z$ axis and the relative position of the two particles given by $\mathbf{r}$. As we are only concerned about the leading order corrections for small values of $s$, we can consistently neglect the weak lattice effect on the interaction potential and use the free-space formula

[TABLE]

In contrast to the contact interacting case, the integrals cannot be solved analytically for arbitrary $\beta$. However, the previous analysis suggested to view the limit $\beta\to\infty$, or equivalently $l/\xi\ll 1$, as nothing else but a gas in free space with an anisotropic effective mass. In the following, we make use of this simplification and derive the corrections due to the weak lattice potential along an arbitrary direction. In the end, we discuss two special limits when the lattice is parallel and perpendicular to the polarization axis of the dipoles.

Since the dipole potential is invariant under rotation around the $z$-axis, the final results can only depend on the angle $\eta$ between the direction of the polarization of the dipoles, which is assumed to be the $z$-axis, and the lattice vector $\mathbf{q}_{L}$. For simplicity, we choose the lattice wave vector to be in the $yz$-plane. Introducing the effective mass $m_{\rm eff}=m/\gamma^{2}$ with $\gamma^{2}=1-s^{2}/32$, the dispersion relation reads

[TABLE]

The correction to the ground-state energy can now be calculated according to (15) with the difference that we integrate over the whole momentum space and use (44) and (43) for the single-particle dispersion and interaction potential, respectively. The integrals are most conveniently performed in spherical coordinates and the beyond-mean-field corrections for a dipolar Bose gas in the presence of a weak one-dimensional lattice read as (see Appendix B for details)

[TABLE]

with $F(\varepsilon_{dd})=\frac{1}{2}\int_{-1}^{1}du\,(1+\varepsilon_{dd}(3u^{2}-1))^{5/2}$ and $H(\varepsilon_{dd})=\frac{1}{2}\int_{-1}^{1}du\,(1+\varepsilon_{dd}(3u^{2}-1))^{5/2}(3u^{2}-1)$.

The energy correction thus has the following form: The first term inside the curly brackets is the usual LHY correction in the case of a dipolar gas in free space Lima and Pelster (2012). The second term accounts for the influence of the lattice and has two parts. The first term is isotropic, has the same structure as in free space and can be explained by an effective mass in one direction. The other part proportional to $H(\varepsilon_{dd})$ is clearly anisotropic in the sense that it depends on the orientation of the lattice with respect to the dipoles. In Fig. 6, we plot the function $F(\varepsilon_{dd})+1/2(3\cos^{2}\eta-1)H(\varepsilon_{dd})$ for different values of the tilting angle $\eta$ between lattice wave vector and polarization of the dipoles. When the lattice is oriented parallel to the dipoles [$\eta=0$, blue (solid) line], the corrections are enhanced and monotonically increase with increasing $\varepsilon_{dd}$. In this setup, the dipoles are directed by the trapping potential to arrange in a side-by-side configuration so that the fluctuations have mainly repulsive character. On the other hand, when the lattice is orientated perpendicular to the dipoles [$\eta=\pi/2$, red (dotted) line], the correction first decreases for small $\varepsilon_{dd}$, reaches a minimum, and finally increases for larger $\varepsilon_{dd}$. In contrast to the parallel orientation, the dipoles are now dragged to the head-tail configuration such that the attractive character of the fluctuations is enhanced. In this case, the correction is also much smaller than in the case where the lattice is oriented parallel to the dipoles. At the ”magic angle” $\eta=\arccos(1/\sqrt{3})\approx 54.7\,{}^{\circ}$, the anisotropic correction vanishes and only the isotropic correction contributes [green (dashed) line].

III.4 Mean-field terms

In the above discussion, we have so far neglected the contributions coming from the mean-field terms which also have an isotropic and anisotropic part where the latter comes from the orientation of the lattice relative to the dipoles. The mean-field terms can be written as

[TABLE]

The first term is the usual contribution from the chemical potential in free space in the absence of any lattice. The second part is the leading-order correction to the mean-field energy in the presence of a weak lattice and provides an anisotropic correction. The anisotropic correction has the same functional dependence on the angle between the lattice wave vector and the polarization axis of the dipoles but, apart from the different scaling in the density, a different functional dependence on the relative dipole interaction strength $\varepsilon_{dd}$. While the mean-field term goes linearly with $\varepsilon_{dd}$ for all values, the function $H(\varepsilon_{dd})\approx 2\varepsilon_{dd}+6/7\varepsilon_{dd}^{2}$ is linear only for small values of $\varepsilon_{dd}$.

IV Conclusions

In this paper, we studied the effects of an optical lattice on the beyond-mean-field corrections for a dipolar Bose gas with emphasis on the ability to control and manipulate the strength of these corrections with respect to the depth and the orientation of the lattice. In the case of a deep three-dimensional lattice, the presence of the lattice introduces a nontrivial density dependence of the beyond-mean-field terms whose form is independent of the strength of the dipolar interaction but whose strength can be enhanced by increasing the dipolar interaction strength. For a weak one-dimensional lattice, we find that the strength of the beyond-mean-field corrections can be controlled by manipulating the relative orientation of the lattice and the dipole axis. In view of the current experiments on dipolar quantum droplets, we would like to point out that our results of the study of the deep optical lattice are not directly applicable to these droplets due to the missing cancellation of the mean-field contributions. However, our results for the weak one-dimensional case might be of importance for future experiments with droplets in weak lattices.

V Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft (DFG) within the research unit FOR 2247.

Appendix A Dipolar interactions in a deep lattice

A.1 Lattice Fourier transform

In this section, we present the calculation of the discrete Fourier transform of the dipolar potential as given in (13). The calculation is based on the analogous one performed in two dimensions Peter (2015) extended to three dimensions. Note that in this section, we denote the summation over lattice sites by $\sum_{\mathbf{R}}$ and measure all lengths and momenta in terms of the lattice spacing and the inverse lattice spacing, respectively.

The lattice Fourier transformation of the dipolar interaction can be written as

[TABLE]

which we may rewrite as

[TABLE]

with

[TABLE]

Note that in three dimensions, both $\chi^{3}(\mathbf{k})$ as well as $\frac{\partial^{2}}{\partial k_{z}^{2}}\chi^{5}(\mathbf{k})$ are divergent for $\mathbf{k}\to 0$. However, as will be shown, the sum of both contributions leads to a finite result which is non-analytic for $\mathbf{k}\to 0$.

Calculation of $\chi^{n}$

In order to calculate $\chi^{n}$, we note the following identity

[TABLE]

where $\Gamma(m)$ denotes the Gamma function. In the end, we will be interested in $n=3$ and $5$. Using the above identity and splitting up the integral into regions from [math] to $\eta$ and $\eta$ to $\infty$, we arrive at

[TABLE]

In the second step, we have rescaled $u$ by $\eta/u$ for the first integral and by $\eta u$ in the second one. The parameter $\eta$ is the so-called Ewald parameter and determines the boundary between the summation in real space and the summation in the momentum space. In the end, this parameter should be chosen such that convergence is achieved rapidly for both sums. However, the result is independent of the choice of $\eta$.

Next, we use Poisson’s summation formula to turn the sum of the first part in (51) into a sum in momentum space. In $d$ dimensions, Poisson’s summation formula applied to our case reads as

[TABLE]

With $a=\eta/u$ and $d=3$, we arrive at

[TABLE]

with the exponential integral function

[TABLE]

Setting $n=3$ and choosing $\eta=\pi$, we get

[TABLE]

Note that the function $E_{1}(x)$ diverges logarithmically for $x\to\infty$. This is expected as the sum $\sum_{\mathbf{R}\neq 0}1/R^{3}$ diverges logarithmically in three dimensions.

Calculation of the anisotropic part

Now we turn to the anisotropic part given by $\frac{\partial^{2}}{\partial k_{z}^{2}}\chi^{5}(\mathbf{k})$. For $\eta=\pi$ and $n=5$, we get

[TABLE]

Differentiating the last term twice with respect to $k_{z}$, we obtain

[TABLE]

For the first term, we note the relation

[TABLE]

and thus have

[TABLE]

Finally, we end up with

[TABLE]

In this expression, we again encounter the diverging part $E_{1}(x)$ with exactly the prefactor needed in order to cancel the contribution from the isotropic $1/R^{3}$ part.

Combining everything and using the explicit form of $E_{0}(x)=e^{-x}/x$, the final result reads as

[TABLE]

The first part exactly reproduces the continuous Fourier transform of $V_{dd}(\mathbf{k})$ in the limit $\mathbf{k}\to 0$, while the second part, including the sums, vanishes in this limit as is shown in the next paragraph. Note in addition that this expression can be evaluated numerically very efficiently since the summations converge very quickly. For numerical purposes, it is useful to replace $\mathbf{q}=2\pi\mathbf{R}$ in the second summation such that convergence is achieved using only a few lattice sites.

In Figs. 7 and 8, the angular dependence of (61) is plotted for different values of $k$ and one can see that in contrast to the free space result, the Fourier transform also depends on the absolute value of $\mathbf{k}$. The corrections to the free-space Fourier transform become even more apparent when plotting both functions over the first Brillouin zone of a simple cubic lattice along paths of high symmetry as shown in Fig. 3. Around $\mathbf{k}=0$, there is little difference between both functions, whereas for finite momenta, deviations become clearly visible since the anisotropic structure of the lattice is probed.

Continuous limit $\mathbf{k}\to 0$ and non-analytic behavior

In the long wavelength limit $\mathbf{k}\to 0$, the first term in (61) reduces to the continuous Fourier transform $V_{dd}(\mathbf{k})=g\,\varepsilon_{dd}(3\cos^{2}\theta-1)$, where $\theta$ is the angle between $\mathbf{k}$ and the $z$-axis, which is assumed to be the axis along which the dipoles are polarized.

The other terms in (61) for $\mathbf{k}\to 0$ read as

[TABLE]

where we replaced $\mathbf{q}=2\pi\mathbf{R}$. In order to simplify this expression further, we make use of the recurrence relation

[TABLE]

Choosing $m=3/2$, we get

[TABLE]

due to the symmetry of the lattice. Thus, we end up with

[TABLE]

This is exactly the same result as in free space and is non-analytic for $\mathbf{k}\to 0$.

A.2 Role of the next-order terms

For the analysis of the beyond-mean-field terms in the main text, we only included the density-density interactions. However, dipolar interaction gives rise to terms like density-induced tunneling or pair hopping. Including the interaction-induced nearest-neighbor couplings results in Sowiński et al. (2012); Wall and Carr (2013); Dutta et al. (2015)

[TABLE]

Transforming this to the quasi momentum space, one obtains

[TABLE]

with

[TABLE]

Here, the sums over $\boldsymbol{\delta}$ are perfumed over nearest neighbors only. This Hamiltonian has the same structure as the one without additional terms but with modified effective interaction $f(\mathbf{q},\mathbf{q}^{\prime},\mathbf{k}’)$. We can thus perform the Bogoliubov approximation, assuming the presence of the condensate at zero momentum and replacing $b_{0}\to\sqrt{N_{0}}$, where $N_{0}$ is the number of particles in the condensate with $N=N_{0}+\frac{1}{2}\sum^{\prime}_{\mathbf{p}}{b^{\dagger}_{\mathbf{p}}b_{\mathbf{p}}+b^{\dagger}_{-\mathbf{p}}b_{-\mathbf{p}}}$, where the prime denotes omitting the zero-momentum mode. For the case of 3D cubic lattice this gives

[TABLE]

The relevant values for a cubic lattice in 3D are

[TABLE]

Here, we defined $\sum_{\boldsymbol{\delta}}e^{i\boldsymbol{\delta}\mathbf{p}}\equiv c_{\mathbf{p}}$ and $c_{-\mathbf{p}}=c_{\mathbf{p}}$ for a cubic lattice.

Now, we perform the Bogoliubov transformation and arrive at the formula for the ground state energy

[TABLE]

with

[TABLE]

and

[TABLE]

Defining

[TABLE]

the Bogoliubov spectrum takes the form

[TABLE]

similar to the standard Bose-Hubbard model, where $\tilde{\epsilon}_{\mathbf{p}}=-2t(c_{\mathbf{p}}-3)$ and $\tilde{U}_{\mathbf{p}}=U$. One can also see that the spectrum is gapless and linear for small momenta

[TABLE]

with the renormalized hopping amplitude $\tilde{t}=t+2nT-nP$ which is now density dependent. The sound velocity is given by $c=\sqrt{2\tilde{t}\tilde{U}_{0}}$. The renormalized hopping amplitude also renormalizes the effective mass which is now given by $m_{\rm eff}=\hbar^{2}/2l^{2}\tilde{t}$. We see that the modifications resulting from additional terms do not fundamentally change the properties of the superfluid. In our case the lattice is assumed to be deep and the role of tunneling terms in the interaction is perturbative.

Appendix B Dispersion in the presence of a weak optical lattice

In this appendix, we first state some important results from the physics of a single particle with mass $m$ in the presence of a weak optical lattice. We restrict our discussion to the case of a one-dimensional lattice along the $z$-axis of the form

[TABLE]

where $q_{L}=\pi/l$ is the lattice vector and $l$ is the lattice period. The energy scale of the lattice is given by the recoil energy $E_{R}=\hbar^{2}q_{L}^{2}/2m$. The full spectrum of the resulting single-particle Hamiltonian can be obtained by diagonalization which is in general only possible numerically. However, since we are interested only in weak lattices with $s=U/E_{R}\ll 1$, we can calculate the spectrum analytically using perturbation theory in the parameter $s$. In second order, the correction to the free-space energy is given by

[TABLE]

where $K$ is a vector of the reciprocal lattice, $K=2nq_{L}$ with integer $n$, the free-space dispersion $\epsilon^{0}_{k}=\hbar^{2}k^{2}/2m$, and the Fourier transform of the lattice potential, $U_{K}=U/2(\delta_{K,2q_{L}}+\delta_{K,-2q_{L}})$. Note that we omit the constant shift $U/2$ which leads to a shift in the chemical potential and in the end, we choose the dispersion such that $\epsilon^{s}_{k=0}=~{}0$. Setting $k=zq_{L}$, $z\in[-1,1]$, we can write

[TABLE]

where $n$ denotes the index of the Bloch band.

Note that (86) only holds for non-degenerate energies away from the edges and the center of the Brillouin zone and in a region of energies where $|\epsilon^{0}_{k-K}-\epsilon^{0}_{k-K^{\prime}}|\gg U_{K}$. Close to the edges of the Brillouin zone, $k=\pm q_{L}$, energies are degenerate and non-degenerate perturbation theory cannot be applied. The dispersion relation in this case reads as

[TABLE]

The above results suggest a splitting into a degenerate and non-degenerate region for the lowest two bands where the border of both regions is determined by the condition $\epsilon^{0}_{k-K}-\epsilon^{0}_{k-K^{\prime}}=U_{K}$ leading to $z=1-s/8$. In the following, we will assume that $z$ is positive but not restricted to values smaller than $1$ as indicated above. Thus, the integration has to be split into three regions. The two non-degenerate regions are from $z=0$ to $1-s/8$ and from $z=1+s/8$ to $\infty$ 222Note that to leading order we only consider a splitting up between the lowest and the second lowest Bloch band at $k=\pm q_{L}$ ($z=\pm 1$) while all other splittings are ignored.. The degenerate region extends from $z=1-s/8$ to $1+s/8$.

B.1 Contribution from the degenerate region

Here, we argue that the contribution from the degenerate region will not contribute in leading order in $s$. In order to see that, we expand (89) for $z^{2}-(z-2n)^{2}\ll s^{2}$ which leads to

[TABLE]

for $1-s/8\leq z\leq 1$, and

[TABLE]

for $1\leq z\leq 1+s/8$.

One can now show that the contribution to the integral (96) coming from the degenerate region is of order $s^{3}$. Next, we show that the leading-order corrections are indeed of order $s^{2}$ and come exclusively from the non-degenerate region.

B.2 Contributions from the non-degenerate region

For the non-degenerate region which ranges from $z=0$ to $1-s/8$ and from $z=1+s/8$ to $\infty$, we have to consider an integral of the form

[TABLE]

where $h(z,0)=0$. For small $s$, we can expand the integral

[TABLE]

The constant term $I(0)$ as well as the linear term vanish due to the structure of $h(z,s)$. The quadratic term, using again that $h(z,0)=0$, yields

[TABLE]

where P.V. denotes the Cauchy principal value. When evaluating the integral, divergences will arise in the vicinity of $z=1$ since the non-degenerate approach is no longer valid. However, the divergences for $z>1$ and $z<1$ will cancel each other.

B.3 Expansion for $z\ll\beta$

Since we are also interested in the limit $\beta\to\infty$ and claim that this corresponds to essentially having an effective mass in the direction of the lattice, we look for an expansion of Eq. (96) for $\beta\to\infty$ or, to be more rigorous, for $z\ll\beta$. The function $f$ in the integrand plays the role of our lattice dispersion. For small arguments, we can expand Eq. (87) around $z=0$ for the lowest band, $n=0$:

[TABLE]

The first term is just the constant shift that will be subtracted by an additional chemical potential. The second term gives the leading behavior and corresponds to having an effective mass in the $z$ direction, $m_{z}~{}=m/(1-s^{2}/32)$.

B.4 Beyond-mean-field corrections in a weak lattice

B.4.1 Contact interaction

In a next step, we present the details on how to calculate the beyond-mean-field corrections for a contact-interacting gas in a weak lattice. For this purpose, we consider (35) and, in what follows, we set $\tilde{\epsilon}_{k}=E_{R}f(kl)$, where $f(kl)$ is a dimensionless function of the quasi-momentum along the direction of the lattice. Using this form of the lattice dispersion and changing to dimensionless variables $k_{\perp}\xi=u$ and $kl=z$, we can transform (35) into

[TABLE]

The factor in front of the integral is proportional to the usual LHY term $\sim m^{3/2}(ng)^{5/2}/\hbar^{3}$, while the remaining part accounts for the influence of the lattice. As we are now interested in a weak lattice, we calculate the integrals in leading order in $s$. It turns out that the leading order is proportional to $s^{2}$ and the corrections read

[TABLE]

B.4.2 Dipolar interaction

In this part, we present the details of the calculation of (45). First, we express the dispersion relation for an anisotropic effective mass in spherical coordinates, that is

[TABLE]

The correction to the ground state energy can now be calculated using (15) according to which

[TABLE]

where $E_{k}=\sqrt{\epsilon_{k}(\epsilon_{k}+2nV_{k})}$ is the Bogoliubov spectrum and $\int d\Omega=\int_{0}^{\pi}d\theta\,\sin\theta\int_{0}^{2\pi}d\phi$ denotes the integration over the solid angle .

The integral on the right-hand side can be simplified and made dimensionless by pulling out the factor $ng$ and using the healing length $\xi^{2}=\hbar^{2}/2mng$. Making the substitution $(k\xi)^{2}f=x^{2}$ and performing the integral over $k$ leads to

[TABLE]

The remaining integrals cannot be solved analytically in general. However, we are only interested in the lowest-order corrections which are due to the lattice such that we can expand the function under the integral to the lowest order in $s$, which leads to

[TABLE]

The first term gives rise to the standard LHY correction induced by dipolar interactions Lima and Pelster (2012), while the latter part includes the corrections due to the lattice. After performing the integral over the angle $\phi$, the second-order correction to the ground state energy in the presence of an optical lattice reads as

[TABLE]

with $F(\varepsilon_{dd})=\frac{1}{2}\int_{-1}^{1}du\,(1+\varepsilon_{dd}(3u^{2}-1))^{5/2}$. The latter term, accounting for the influence of the lattice, can be rewritten as

[TABLE]

with $H(\varepsilon_{dd})=\frac{1}{2}\int_{-1}^{1}du\,(1+\varepsilon_{dd}(3u^{2}-1))^{5/2}(3u^{2}-1)$.