# Dynamics of cold bosons in optical lattices: effects of higher Bloch bands

Mateusz Łącki, Dominique Delande, Jakub Zakrzewski

► To cite this version:

Mateusz Łącki, Dominique Delande, Jakub Zakrzewski. Dynamics of cold bosons in optical lattices: effects of higher Bloch bands. New Journal of Physics, 2013, 15, pp.013062. 10.1088/1367-2630/15/1/013062. [hal-01587054](https://hal-sorbonne-universite.fr/hal-01587054)

HAL Id: hal-01587054

https://hal.sorbonne-universite.fr/hal-01587054

Submitted on 13 Sep 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Dynamics of cold bosons in optical lattices: effects of higher Bloch bands

To cite this article: Mateusz cki et al 2013 *New J. Phys.* **15** 013062

View the [article online](#) for updates and enhancements.

## Related content

- Tutorial

C Trefzger, C Menotti, B Capogrosso-Sansone et al.

- Non-standard Hubbard models in optical lattices: a review

Omjyoti Dutta, Mariusz Gajda, Philipp Hauke et al.

- Multi-orbital and density-induced tunneling of bosons in optical lattices

Dirk-Sören Lühmann, Ole Jürgensen and Klaus Sengstock

## Recent citations

- Loading ultracold gases in topological Floquet bands: the fate of current and center-of-mass responses

Alexandre Dauphin et al.

- Hubbard Model for Atomic Impurities Bound by the Vortex Lattice of a Rotating Bose-Einstein Condensate

T. H. Johnson et al.

- Physics of higher orbital bands in optical lattices: a review

Xiaopeng Li and W Vincent Liu
# Dynamics of cold bosons in optical lattices: effects of higher Bloch bands

Mateusz Łącki¹, Dominique Delande² and Jakub Zakrzewski¹,³,<fnref n="4" />

¹ Instytut Fizyki imienia Mariana Smoluchowskiego, Uniwersytet Jagielloński, ulica Reymonta 4, PL-30-059 Kraków, Poland

² Laboratoire Kastler-Brossel, UPMC-Paris 6, ENS, CNRS; 4 Place Jussieu, F-75005 Paris, France

³ Mark Kac Complex Systems Research Center, Uniwersytet Jagielloński, ulica Reymonta 4, PL-30-059 Kraków, Poland
E-mail: kuba@if.uj.edu.pl

New Journal of Physics 15 (2013) 013062 (21pp)

Received 16 November 2012

Published 29 January 2013

Online at http://www.njp.org/

doi:10.1088/1367-2630/15/1/013062

**Abstract.** The extended effective multiorbital Bose-Hubbard-type Hamiltonian which takes into account higher Bloch bands is discussed for boson systems in optical lattices, with emphasis on dynamical properties, in relation to current experiments. It is shown that the renormalization of Hamiltonian parameters depends on the dimension of the problem studied. Therefore, mean-field phase diagrams do not scale with the coordination number of the lattice. The effect of Hamiltonian parameters renormalization on the dynamics in reduced one-dimensional optical lattice potential is analyzed. We study both the quasi-adiabatic quench through the superfluid-Mott insulator transition and the absorption spectroscopy, that is, the energy absorption rate when the lattice depth is periodically modulated.

<figure><img src="image_2.png" /></figure>

<footnote n="4">⁴ Author to whom any correspondence should be addressed.</footnote>
Contents

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>2. Tight-binding descriptions of an ultra-cold boson gas in an optical lattice</td><td>3</td></tr><tr><td>2.1. Solving the on-site problem</td><td>6</td></tr><tr><td>3. Mean field diagrams for different lattice dimensions</td><td>9</td></tr><tr><td>4. Consequences of coupling constants renormalization for dynamics</td><td>11</td></tr><tr><td>4.1. Modulation of optical lattice—absorption spectroscopy</td><td>11</td></tr><tr><td>4.2. The Florence experiment revisited</td><td>14</td></tr><tr><td>5. Conclusion</td><td>16</td></tr><tr><td>Acknowledgments</td><td>17</td></tr><tr><td>Appendix. Diagonalization of the on-site Hamiltonian</td><td>17</td></tr><tr><td>References</td><td>20</td></tr></table>

1. Introduction

Ultra-cold bosonic atom gases in optical lattices have been an ultra-hot research area in recent years (for recent reviews and an extensive reference list see [1–3]). They provide a means to create and control experimental systems mimicking different condensed matter physics models [4–6]. The interest has been stimulated in part by the fact that there exists [7] an accurate mapping of a continuous Hamiltonian to a Hamiltonian on a lattice—the Bose-Hubbard (BH) Hamiltonian, originally formulated by Gersch and Knollman [8]. The lattice models significantly ease the analytical [9–12] and numerical analyses although promising new ideas were proposed recently that can enable the analysis in continuous variables [13] beyond the mean field level [14, 15].

Bosons in one-dimensional (1D) optical lattices have also been an area of extensive experimental research [16–21]. The corresponding 1D BH Hamiltonian can be effectively addressed numerically by density matrix renormalization group (DMRG) and related techniques [22]. These techniques have broad applications, in particular they enable simulations targeting real-life many-body systems [23], with controllable error and no systematic errors due to unsound assumptions. Their success relies on area laws that control the growth of entanglement [24]. The entanglement is used as a ‘small parameter’ [25, 26] that makes possible the construction of an efficient variational set—the so-called matrix product states (MPS). Several numerical investigations on experimental and ‘close to experimental’ systems have been performed [27–29]. They focused largely on two aspects: a quench through a phase transition and the simulation of absorption spectroscopy [16]. The phase transition from the superfluid (SF) to the Mott insulator (MI) phase occurs when the lattice depth is increased beyond a critical value [4, 30]. The adiabaticity of this process has been addressed in [15, 29, 31]. In the second example, the energy absorption rate is analyzed as a function of the frequency of modulation of the lattice depth [27, 32, 33]. The locations, the number of peaks in the spectrum and their heights are related to the state of the gas: either SF or insulating.

The derivation of the BH Hamiltonian assumes the restriction of physics to the lowest Bloch band of the optical lattice. This assumption is reasonable, as the energy gap between the first and second Bloch bands is in most cases around ten times larger than the energy scale in
the discrete model. It is also much larger than the thermal energy scale $k_B T$ providing additional argument for a zero-temperature analysis. An early analysis [34] suggested that higher bands may be included by an appropriate modification of BH Hamiltonian parameters for large occupation numbers. The effects due to higher Bloch bands have also been studied for lower densities [35–38] in relation to the discovery that higher bands affect the SF-insulator transition in Bose-Fermi mixtures [39, 40]. This gave an explanation for experimental observations of the shift of the SF-MI phase transition [41] which could not be explained by single-band approaches. In [36], the first excited band is included in a two-flavor model; the effects of higher bands could also be built in via an effective three-body interaction in the lowest band [35, 37].

In recent years, another approach for three-dimensional (3D) optical lattices has been proposed [39, 42, 43] (which somehow resembles in spirit [34] while being used for moderate atomic densities). In this approach, the higher bands are included in the on-site Hamiltonian, which is then diagonalized and many-body ground states of this problem for different total numbers of particles are used as a local Hilbert space, replacing the usual Fock basis $|n\rangle$. The multiband lattice Hamiltonian is then expressed in this new basis, yielding an effective single-band model with occupation-dependent parameters—renormalized values of the initial Hamiltonian parameters for the lowest Bloch band. Interestingly, the change in BH Hamiltonian parameters is large in the MI regime, where the energy gap is also large, contrary to naive intuition.

In this paper, we consider the derivation of the BH model’s effective coupling constants, taking a closer look at the underlying, quite challenging, numerical problem—the accurate diagonalization of the on-site Hamiltonian. We study how the dimensionality of the optical lattice affects the renormalization scheme. The dependence on dimension implies that the mean field diagrams of the system no longer depend solely on the coordination number of the lattice. We then take a look at the effects of the renormalization of coupling constants on the dynamics in a 1D optical lattice: both quench through the SF-MI transition and the absorption spectroscopy are analyzed. We find that renormalization of atom-atom interactions significantly shifts and sometimes modifies the absorption peaks. It also appears to have a serious effect on adiabaticity predictions for the SF-MI phase transition.

## 2. Tight-binding descriptions of an ultra-cold boson gas in an optical lattice

The second quantization Hamiltonian for a dilute gas of interacting bosonic atoms in the optical lattice potential $V(\vec{r})$ and external trapping potential $V_e(\vec{r})$ is of the form

$$H = \int d^3\vec{r} \Psi^\dagger(\vec{r}) \left( -\frac{\hbar^2}{2m} \nabla^2 + V(\vec{r}) + V_e(\vec{r}) \right) \Psi(\vec{r}) \\ + \frac{1}{2} \int d^3\vec{r} d^3\vec{r}' \Psi^\dagger(\vec{r}) \Psi^\dagger(\vec{r}') V_{\text{int}}(\vec{r}, \vec{r}') \Psi(\vec{r}) \Psi(\vec{r}'), \quad (1)$$

where $V_{\text{int}}(\vec{r}, \vec{r}')$ is an isotropic short-range pseudopotential modeling s-wave interactions [44]

$$V_{\text{int}}(\vec{r}, \vec{r}') = \frac{4\pi\hbar^2 a_s}{m} \delta(\vec{r} - \vec{r}') \frac{\partial}{\partial|\vec{r} - \vec{r}'|} |\vec{r} - \vec{r}'| \quad (2)$$

with $a_s$ being the scattering length.

The Hamiltonian admits a natural energy scale, set by the recoil energy $E_R = \frac{\hbar^2 k^2}{2m}$, where $k = \frac{2\pi}{\lambda}$ and $\lambda$ denotes the optical lattice wavelength. We assume this energy unit from
now on. The optical lattice potential is typically: $V(x, y, z)/E_R = s_x \cos^2(kx) + s_y \cos^2(ky) + s_z \cos^2(kz)$. If $s := s_x \ll s_y = s_z := s_⊥$, then a classical setup for a 1D optical lattice is obtained [16]. Tunneling in the y- and z-directions is highly suppressed and the system may be considered as a series of independent 1D tubes along the x-direction. Similarly, when $s := s_x = s_y \ll s_z := s_⊥$, a two-dimensional (2D) optical lattice is obtained. When $s_x = s_y = s_z := s$ the potential corresponds to a 3D lattice. Let us emphasize that all calculations of renormalized parameters presented below are truly 3D ones, the labels ‘1D’ or ‘2D’ just referring to situations where $s$ and $s_⊥$ are very different.

Jaksch et al [7] in their seminal paper introduced a mapping of (1) onto the BH Hamiltonian (here $V_i$ is the local energy shift due to $V_e: V_i = V_e(\vec{r}_i)$):

$$H_{\text{BH}} = -J \sum_{\langle ij \rangle} \hat{a}_i \hat{a}_j^\dagger + \text{h.c.} + \frac{U}{2} \sum_i \hat{n}_i (\hat{n}_i - 1) - \sum_i \hat{n}_i (\mu - V_i). \quad (3)$$

The field operator is expanded in the set of the lowest band Wannier functions of the lattice: $\psi(\vec{r}) = \sum w_i^0(\vec{r})a_i$—here $w_i^α$ denotes the (real valued) Wannier function localized on site $i$ of the lattice for the $α$th Bloch band (we shall denote sites by roman subscripts and bands by Greek superscripts). The parameters $U$ and $J$ are expressed by the appropriate integrals of Wannier functions: $J = -\int w_i^0(\vec{r})(-\frac{\hbar^2}{2m}\nabla^2 + V(\vec{r}))w_j^0(\vec{r}) d\vec{r}$ (where $i$ and $j$ are neighboring sites) and $U = \frac{4\pi\hbar^2 a_s}{m} \int w_i^0(\vec{r})^4 d\vec{r}$.

However, it turns out that such a single-band approximation is often insufficient for realistic values of the parameters, and that contributions of higher bands cannot be neglected. This problem has been studied in the literature under various conditions and using various methods. Of special interest is the situation where the scattering length $a_s$ is large (near a Feshbach resonance), where the optical lattice strongly modifies the effective atom–atom interaction, see, e.g., [45–53]. This is not the situation realized for the $^{87}$Rb atom in zero magnetic field used in the Florence experiment [18, 39], where $a_s = 5.2$ nm, much shorter than the lattice spacing $λ/2 = 377$ nm.

In this paper, we will consider only situations where the atomic wavenumber $k$—being evaluated either in the lowest band or in the excited bands included in the calculation—is such that $ka_s \ll 1$ so that only s-wave low-energy interatomic scattering is relevant. In practice, this puts a limit on the number of bands used $B < 20$, beyond which the model is not a good approximation of the real world.

It is well known that, in dimensions higher than 1, contact interactions cannot be modeled by a $δ$-potential, but require a specific regularization to avoid artificial divergences. Mathematically, a self-adjoint extension of the original Hamiltonian is needed [54]. It boils down to a dimension-dependent regularization of the $δ$-potential [55], which, in 3D, is the so-called Fermi pseudo-potential [56] used in equation (2). Even with the correct pseudo-potential, one must be careful when expanding over an infinite set of square integrable smooth basis functions (such as the Wannier functions used below) without renormalization of the interaction strength, because it may lead to incorrect results, such as diverging perturbative expansions, see [57] for the specific example of two interacting particles in an harmonic trap. Numerical diagonalization, used in the following, leads to less severe problems as discussed in [58]. We will use a rather small number of bands (up to $B = 15$), so that the highest atomic wavenumbers effectively included in the calculation are still rather small and the divergence of the Green function at short distance is just a small perturbation. In other words, although the method we use leads, in principle, to divergences when $B \to \infty$, these divergences manifest themselves
only beyond the largest B used in our calculations. Note also that the single-band approximation described above is free of this problem and pure δ interactions can be used in these cases.

Without restriction to the lowest Bloch band, the expansion in the full Wannier functions basis set would give a multiband variant of the BH Hamiltonian

$$H = -\sum_{\alpha, i, j} J_{ij}^{\alpha\alpha}(\hat{a}_i^{\dagger}\hat{a}_j^{\dagger} + \text{h.c.}) + \sum_{\substack{\alpha, \dots, \delta \\ i, \dots, l}} U_{ijkl}^{\alpha\beta\gamma\delta}(\hat{a}_i^{\dagger}\hat{a}_j^{\dagger}\hat{a}_k^{\dagger}\hat{a}_l^{\dagger} + \sum_{\alpha, i}(E_{\alpha} + V_i - \mu)(\hat{a}_i^{\dagger}\hat{a}_i^{\dagger}) \quad (4)$$

with

$$U_{ijkl}^{\alpha\beta\gamma\delta} = \frac{4\pi\hbar^2 a_s}{m} \int d\vec{r} w_i^{\alpha}(\vec{r}) w_j^{\beta}(\vec{r}) w_k^{\gamma}(\vec{r}) w_l^{\delta}(\vec{r}) \quad (5)$$

and

$$J_{ij}^{\alpha\alpha} = - \int d\vec{r} w_i^{\alpha}(\vec{r}) \left( -\frac{\hbar^2}{2m} \nabla^2 + V(\vec{r}) \right) w_j^{\alpha}(\vec{r}).$$

Note that, because Wannier functions are smooth, the potential (2) may be replaced by a contact Fermi potential

$$V_c(\vec{r}, \vec{r}') = \frac{4\pi\hbar^2 a_s}{m} \delta(\vec{r} - \vec{r}') \quad (6)$$

in the integral (5).

By construction $J_{ij}^{\alpha\beta} = 0$ for $\alpha \neq \beta$. For sufficiently deep optical lattices (typically $s > 3$), it is enough to restrict hopping to nearest-neighbor sites, as tunneling amplitudes are exponentially damped with the hopping distance (for shallow lattice next-nearest-neighbors hopping may be necessary—see [21]).

A 3D Wannier function being a product of 1D Wannier functions, the 3D integral in $U_{ijkl}^{\alpha\beta\gamma\delta}$ is a product of three integrals over each coordinate. The interaction parameters differ for each direction, as Wannier functions depend on the lattice depth, which may be different in each direction.

The Hamiltonian (4) is difficult to use in practice, even in 1D systems, because the on-site dimension $d$ of the $n$-particle problem restricted to the lowest $B$ bands is $d = \binom{B^3+n-1}{n}$ (the on-site problem is genuinely 3D even for a quasi-1D model), and moreover, the numerical complexity of the best 1D algorithms scales with $d$ at least as $O(d^3)$. The complexity of more sophisticated approaches (such as multiscale entanglement renormalization ansatz (MERA) and projected entangled pair states (PEPS) [59–62]) is several orders higher.

Thus for computational purposes, one must restrict the local Hilbert space. Assuming that interactions are on-site only, i.e. $U_{ijkl}^{\alpha\beta\gamma\delta} \neq 0$ for $i = j = k = l$ together with considering the lowest Bloch band only ($\alpha = \beta = \gamma = \delta = 0$) [7] leads directly to the BH Hamiltonian, equation (3), provided we chose the zero of the energy axis at $E_0$.

A more sophisticated approach is discussed in [39, 43]. The on-site Hamiltonian, restriction of (4) to a single site

$$\mathcal{H}_{\text{loc}} = H_E + H_U = \sum_{\alpha} E_{\alpha} \hat{n}_{\alpha} + \sum_{\alpha\beta\gamma\delta} U^{\alpha\beta\gamma\delta} \hat{a}_{\alpha}^{\dagger} \hat{a}_{\beta}^{\dagger} \hat{a}_{\gamma} \hat{a}_{\delta} \quad (7)$$

(with $\hat{n}_{\alpha} = \hat{a}_{\alpha}^{\dagger}\hat{a}_{\alpha}$) can be diagonalized to yield a space of $n$ particle ground states. The eigenenergies $\epsilon_0^n$ of the on-site $n$ particle ground states $|\psi_0^n\rangle$ are the starting point in determining
new values of $U$ parameters in the effective Hamiltonian. To define renormalized values of $U$, the energy $\epsilon_0^n$ has to be decomposed into the interaction energy (which, in the case of the BH model, is just $\frac{U}{2}n(n-1)$) and a single-particle energy; (which, in the BH case shifts $\mu$ by the lowest Bloch band energy). The most natural way to define the interaction energy would be to use

$$ \frac{U_n}{2}n(n-1) = \langle \psi_0^n | \sum_{\alpha\beta\gamma\delta} U^{\alpha\beta\gamma\delta} \hat{a}_{\alpha}^{\dagger} \hat{a}_{\beta}^{\dagger} \hat{a}_{\gamma} \hat{a}_{\delta} | \psi_0^n \rangle. \quad (8) $$

Unfortunately, $U_n$ cannot be defined in such a way if we request the Hamiltonian to have a form resembling equation (3). That is because the single-particle energy is no longer a linear function of $n$. This can be circumvented by defining $U_n$ via

$$ \epsilon_0^n = \frac{U_n}{2}n(n-1) + nE_0. \quad (9) $$

This definition makes Hamiltonians (3) and (10) similar. But $U_n$ is no longer the interaction energy, it also contains contributions of higher Bloch bands population to the single-particle energy. From now on, we use definition (9). Note that $U_n$ depends nontrivially on the geometry of the lattice.

The second stage is to reintroduce inter-site couplings. Even if only one band is taken into account (so that $U_n$ is simply $U$), the inter-site interaction $U_{iiij}^{0000}$ induces an effective coupling, which is proportional to the density, that is, a term $U_{iiij}^{0000}a_i a_j^\dagger(n_i + n_j - 1) + \text{h.c.}$ in the effective Hamiltonian (called bond-charge term in [39, 63]). Similarly to the original tunneling term, it is important only between nearest neighbors as soon as $s$ is larger than unity. In the low-density regime (typically $n < 7$), this contribution leads to an increase (since $U_{iiij}^{0000} < 0$ for $s > 1$) of the tunneling amplitude $J \to J - U_{iiij}^{0000}(n_i + n_j - 1)$. The correction is at most of the order of the raw tunneling. When higher bands are taken into account, the modification of $|\psi_0^n\rangle$ induces a renormalization of the standard tunneling term (as well as of the bond-charge term) which becomes also dependent on the occupation numbers of sites between which tunneling occurs.

The effective multiorbital (EMO) Hamiltonian finally becomes

$$ H^{\text{EMO}} = -\sum_{\langle i,j \rangle} \left( a_i a_j^\dagger \sum_{n_i, n_j} J_{n_i, n_j} P_{n_i}^i P_{n_j}^j + \text{h.c.} \right) + \sum_{n,i} \frac{U_n}{2} n(n-1) P_n^i, \quad (10) $$

where $P_n^i = |i, n\rangle\langle i, n|$.

This Hamiltonian will allow us to study the influence of higher bands on the dynamics later. First we discuss the accurate numerical determination of the $U$ and $J$ parameters that, in itself, is a challenging problem, giving an insight into the physics involved.

## 2.1. Solving the on-site problem

The single-site problem is equivalent to finding the $n$ particle ground state of the Hamiltonian (7). If the lowest $B$ Bloch bands are included (typical values of $B$: 4 [42], 9 [39]), the problem quickly becomes too involved computationally to be exactly diagonalized, and truncation of the basis has to be performed.

For small particle numbers, the problem is dominated by the $H_E$ term. The direct transition from the lowest to the first excited band is forbidden by symmetry consideration so the relevant energy scale is $\Delta = E_2 - E_0$. Promoting a particle from the lowest to the second excited
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Comparison of the effective on-site interaction strength $U_n$ obtained using diagonalization of the on-site Hamiltonian on two different basis sets with the same size equal to 40 000. If basis vectors are chosen according to their energy (dashed lines), false saturation effects appear. Estimating the influence by a perturbative-like scheme (solid lines) does not seem to suffer from saturation effects. The 3D case is considered: $s = s_⊥ = 34.8$.</figcaption></figure>

band, via interaction term, is proportional to $n^{3/2}U^{2000}$ where $U^{2000}$ is the corresponding interaction integral. This yields for the validity of the perturbative approach the condition $n^{3/2}U^{2000} \ll \Delta$. For typical parameters corresponding to Rb scattering length and $s$ of the order of 30, $\Delta/U^{2000} \approx 60$, yielding the limiting value of $n \ll 15$.

For small occupation numbers, a perturbative approach seems justified. Perturbation theory enables us to estimate the impact each excited vector $|\psi_p\rangle$ has on the ground-state energy. Let $|\psi_0\rangle$ be an $n$-boson ground (Fock) state of $H_E$. The larger the matrix element $|\langle\psi_0|\mathcal{H}_{\text{loc}}|\psi_p\rangle|^2$ and the smaller the energy $\langle\psi_p|\mathcal{H}_{\text{loc}}|\psi_p\rangle$, the larger is the impact. The perturbative scheme provides a hint at how to choose an ‘optimal’ subset of the basis in which the full problem could be diagonalized. Since the exact diagonalization in the variational basis is the last step, we do not follow perturbation theory exactly, but just use it to choose a close to optimal basis, not to calculate the energy correction. Details of basis generation and variational space sampling are given in the appendix. A more traditional method [39] is to choose a subset according to the least energy principle—with minimal $\langle\psi_p|\mathcal{H}_{\text{loc}}|\psi_p\rangle$.

We have performed a detailed analysis comparing both methods for $s = s_⊥ = 34.8$, a strongly coupled case (3D optical lattice, $2a_s/\lambda = 0.014$, $\lambda = 754$ nm—parameters taken as typical values from [39]). We choose a system with $n = 2-5$ particles and 40 000 basis vectors according to both least energy (as in [39]) and the perturbative method (39 900 vectors are generated within the first-order, and 100 within the second-order perturbative scheme).

The ground-state energy—obtained from numerical diagonalization of the on-site Hamiltonian in a restricted subset—versus the number $B$ of bands included is shown in figure 1, for various $n$. The least energy method clearly gives the false impression of saturation of results if $B \approx 7-9$ bands are included. The false saturation occurs because the least-energy method does not evaluate $|\langle\psi_0|\mathcal{H}_{\text{loc}}|\psi_p\rangle|^2$ and, therefore, fills the variational basis with low-energy
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Renormalized interaction parameters $U_n$ versus the strength of the optical lattice, for various dimensionalities: black, green and red lines corresponding to 1D, 2D and 3D, respectively. Interaction and lattice parameters: $2a_s/\lambda = 0.014$, $\lambda = 754$ nm. The transverse lattice height is $s_\perp = 34.8$. The curves meet at $s = s_\perp$.</figcaption></figure>

irrelevant vectors with vanishingly small matrix elements. The perturbative-like approach does not show similar saturation effects and, moreover, suggests that linear extrapolation of the results may be performed. We find that the best compromise between computational effort and accuracy is to perform extrapolation of results as a function of $1/B$ by means of the ansatz $U_n(B) = U_n^\infty + c_0/B$. The same extrapolation scheme can be used for the $J$ parameters (but leads to less drastic modification of the results). As mentioned above, numerical diagonalization in the set of Wannier functions might lead to unphysical divergences as $B \to \infty$ because of the subtle properties of contact interactions in 3D. One could expect a $1/B$ divergence for large $B$. Figure 1 does not show any indication of such a divergence, which could be visible for larger $B$.

The $U_n$ parameters for the 1D, 2D and 3D lattices are presented in figure 2, while the renormalized tunneling amplitudes are shown in figure 3. We find that, in low dimensions, $U_n$ vary less with $s$ compared to the full 3D lattice. The high transverse lattice causes significant renormalization of $U_n$ even for small $s$, as $s_\perp$ is still large.

Inspection of figures 1 and 2 shows that the difference between consecutive $U_n$ is approximately constant, i.e.

$$U_{n-1} - U_n \approx W \qquad (11)$$

at least for low densities (for typical lattice parameters, $W$ is constant up to 10%). This is easily understood: the alternative effective theory of [35] expresses the correction to the on-site energy term via a three-body interaction term (equation (12) of [35]) and $W$ is simply related to their parameter $\tilde{U}_3$. The deviations from the linear form, equation (11), may be then related to higher order terms in [35], i.e. four-body term, etc. Similarly, the lowest order perturbative term discussed above gives a correction to the interaction energy term $n(n-1)U/2$ of the order $n^3[U^{2000}]^2/\Delta$. That yields a crude estimate for $W \approx U^{2000}/\Delta \approx 1/60$ in good agreement with figure 1.
<figure><img src="image_5.png" /><figcaption>Figure 3. Renormalized tunneling amplitudes $J_{n_1,n_2}$ versus the strength of the optical lattice, for various dimensionalities: black, green and red lines corresponding to 1D, 2D and 3D, respectively. Interaction and lattice parameters: $2a_s/\lambda = 0.014$, $\lambda = 754$ nm, $s_\perp = 34.8$. The curves meet at $s = s_\perp$.</figcaption></figure>

### 3. Mean field diagrams for different lattice dimensions

A homogeneous (without external trap, $V_i = 0$) Hamiltonian may be taken to the thermodynamic limit. Then the particle density is determined by the chemical potential $\mu$ and the Hamiltonian parameters: $J$ and $U$. An MI phase is determined by integer density $\langle \hat{n} \rangle$ and noncompressibility: $\frac{\partial \langle \hat{n} \rangle}{\partial \mu} = 0$. The rest of the phase diagram is the compressible SF phase [9, 11, 64, 65]. We perform now the mean field analysis of the phase diagrams of Hamiltonians (3) and (10). A Gutzwiller analysis of the BH model is a variational minimization of the following functional (i.e. mean ground-state energy):

$$H^{\text{BH}}[\psi] = \langle \psi | H^{\text{BH}} | \psi \rangle = -2zJ \langle \hat{a}_i \rangle \langle \hat{a}_j \rangle + \frac{U}{2} \langle \hat{n}(\hat{n}-1) \rangle - \mu \langle \hat{n} \rangle \quad (12)$$

using the Gutzwiller ansatz: $|\psi\rangle = \bigotimes_{i} |\psi_i\rangle$ with $|\psi_i\rangle = \sum_{n} f_n |n\rangle$ the on-site wavefunction. The influence of the lattice geometry is reduced only to the coordinate number, $z$, in the first term, as $\langle \hat{a}_i \rangle = \langle \hat{a}_j \rangle$ due to the translational invariance in the thermodynamic limit. The phase diagram depends only on the single parameter $\frac{zJ}{U}$.

For the EMO Hamiltonian, the dependence of interaction parameters on the dimensionality of the optical lattice is nontrivial. We shall use the data from figure 2. Thanks to the translational invariance, the Gutzwiller mean-field approach to the effective Hamiltonian (10) is equivalent to the minimization of the following functional:

$$H^{\text{EMO}}[\psi] = -2z \sum_{n_1, n_2} J_{n_1, n_2} \langle \psi_l | a_i^\dagger | n_1 \rangle \langle n_1 | \psi_l \rangle \langle \psi_l | a_i | n_2 \rangle \langle n_2 | \psi_l \rangle \\ + \sum_n \frac{U_n}{2} n(n-1) |\langle \psi_l | n \rangle|^2 - \mu \sum_n n |\langle \psi_l | n \rangle|^2. \quad (13)$$

Clearly, a single parameter is no longer sufficient to describe the mean field problem. Let us denote by $J_{\text{BH}}$ the tunneling amplitude and by $U_{\text{BH}}$ the interaction of the standard BH
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Mean-field phase diagrams for 1D, 2D and 3D lattices. Different curves denote borders between MI and SF phases. Dashed black lines correspond to the standard BH model for any dimension and blue, green and red curves denote 1D, 2D and 3D lattices of the EMO Hamiltonian. Red dashed lines show the result obtained for nine bands as in [39]. The limit $zJ_{BH}(s)/U_{BH}(s)$, $s \to \infty$ is different for each dimension. The $s \to \infty$ limit corresponds to the ill-defined situation in which the transverse lattice is shallower than the main lattice (this formal limit is also dimension-dependent). The perpendicular lattice depth is fixed at $s_{\perp} = 34.8$, $\lambda = 754$ nm, $2a_s/\lambda = 0.014$ as appropriate for $^{87}$Rb [39].</figcaption></figure>

Hamiltonian. As $zJ_{BH}(s)/U_{BH}(s)$ is strictly decreasing with the lattice height $s$, it provides a way to plot results calculated for lattice depth $s$ in the $(zJ_{BH}/U_{BH}, \mu)$ coordinate space used for a traditional phase diagram. This mapping allows also to directly compare the results obtained using equations (12) and (13). Figure 4 shows the Gutzwiller phase diagrams for 1D, 2D and 3D lattices. In contrast with the ordinary BH model, there is a nontrivial dependence on the dimension. It is rather small for the first lobe becoming more significant for higher occupation of sites. Let us stress again that phase diagrams for the multiorbital (MO) parameters $J_{n_i,n_j}/U_n$ do not need to be directly related to $zJ_{BH}/U_{BH}$. The physical observable that is common for the ordinary BH phase diagram and the MO effective theory is the lattice depth $s$. Moreover, for dimensions 1 and 2, there is actually a whole family of different phase diagrams indexed by $s_{\perp}$. We show just a single choice for a generic value of $s_{\perp} = 34.8$.

For 3D, a comparison with the mean field diagram obtained for nine bands [39] is possible. The difference is quite small, the difference in $U_n$ (of the order of a few per cent—see figure 1) manifests itself mostly in shifting the borders between different Mott lobes for higher occupation numbers. Recently, a continuous space quantum Monte Carlo calculation of the SF–MI border in the cubic lattice for unit filling was reported for various scattering lengths [66] and compared with the mean field results of [39]. We consider a single ratio of $a_s/\lambda$ as appropriate for $^{87}$Rb; the small difference between our results and that of [39] for unit filling on the scale of figure 1 of [66] is negligible.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** 1D phase diagram obtained using imaginary time evolution and the TEBD algorithm. Black dashed curves present the standard BH 1D case, and the red solid lines are obtained for the EMO 1D model (10) with $s_{\perp} = 34.8$.</figcaption></figure>

In 1D, the mean field approximation is inaccurate. To obtain a reliable phase diagram, we have used energy minimization through imaginary time evolution using the time evolving block decimation (TEBD) [25, 26] algorithm. We fix the lattice size to be $L = 100$ (we have checked that choosing a larger lattice size $L = 200, 300$ and $400$ does not alter the results significantly, except at MI tips, where an approximate finite-size scaling is performed). This is significantly less computationally demanding than using the infinite, translationally invariant version of TEBD [62, 67]. The transverse lattice height is again $s_{\perp} = 34.8$. Let us denote by $E(N, s)$ the ground-state energy of an $N$-particle system for lattice height $s$. We calculate approximations to the critical values of chemical potentials delimiting an MI region with average filling $n$ by $\mu_{+}(s) \approx E(nL + 1, s) - E(nL, s)$, $\mu_{-}(s) \approx E(nL, s) - E(nL - 1, s)$. We plot the phase diagrams for both the BH and the EMO Hamiltonians in figure 5. We again see similar results: the Mott lobes shrink also in 1D, as predicted by the mean-field approach. As shown below, this is also reflected in the dynamical properties.

# 4. Consequences of coupling constants renormalization for dynamics

## 4.1. Modulation of optical lattice—absorption spectroscopy

By periodically modulating the lattice depth, one transfers energy to the atomic sample in the lattice. Absorption spectroscopy—also incorrectly nicknamed modulation spectroscopy— consists in studying the dependence of the energy absorption rate on the modulation frequency. This absorption is sensitive to the quantum phases present in the system, as shown in early experiments [16]. It has been simulated [27, 32, 33] for atoms in an optical lattice in the presence of a harmonic confinement, using a standard 1D BH model. It seems interesting to see whether excited bands affect the absorption spectra. To this end, we consider the real-time evolution of the ground state of a given system at $s = s_0$, exposed to a time-varying lattice height $s(t) = s_0 + s_m \cos \omega t$, $s_m/s_0 \ll 1$. The simulation is performed in the presence of a
harmonic trap $V_i = \kappa (i - i_0)^2$. The energy of the system is measured after some fixed time. The energy gain (per particle) as a function of the modulation frequency yields the absorption spectrum.

For a deep optical lattice in a predominantly Mott insulating phase, the absorption spectrum for the standard BH Hamiltonian consists of a few peaks located at multiplicities of $U$ [27, 32, 33]. The situation is slightly more complex for the EMO Hamiltonian. The position of peaks can be easily determined in the deep Mott regime ($J \to 0$). States excited during the modulation are mainly those that differ from the ground state by nearest-neighbor transfer of one particle [32, 68]. Up to a small correction due to the difference of the local chemical potential $\mu_i = \mu - V_i$, the excitation energy is determined by the occupation numbers of the source site $i$ and destination site $j$. It is

$$ \Delta E(n_i; n_j) = \frac{1}{2} [-U_{n_j} n_j(n_j - 1) + U_{n_j+1} n_j(n_j + 1) \\ -U_{n_i} n_i(n_i - 1) + U_{n_i-1} (n_i - 1)(n_i - 2) + \Delta \mu_{ij}] \quad (14) $$

with $\Delta\mu_{ij} = \mu_i - \mu_j = V_j - V_i$. Nearest-neighbor excitation means that $|i - j| = 1$, and $|n_i - n_j| \le 1$. If $U_n = U$ (BH Hamiltonian), we have that $\Delta E(n_i; n_j) = (n_j - n_i + 1)U + \Delta\mu_{ij}$. By virtue of equation (11), we may approximate equation (14) by $\Delta E(n_i; n_j) = (n_j - n_i + 1)U_{n_i} + (n_i - n_j - 2)(n_i + n_j - 1)W$. For a trapped gas with maximum occupation number $n = 3$, the relevant values are $\Delta E(1, 1) = U_2$, $\Delta E(1, 2) = 3U_3 - U_2 \approx 2U_2 - 3W$, $\Delta E(2, 2) = 3U_3 - 2U_2 \approx U_2 - 3W$, $\Delta E(2, 3) = 6U_4 - 3U_3 - U_2 \approx 2U_2 - 9W$, $\Delta E(3, 3) = 6U_4 - 6U_3 + U_2 \approx U_2 - 6W$, $\Delta E(n, n - 1) = 0$, with $\Delta\mu$ neglected for clarity.

A qualitative comparison of the expected absorption spectra for the standard BH case and the EMO model is possible. The density profile in the large $s$, low hopping, limit shows the well-known wedding cake structure (see figure 6(a)). Weak, periodic modulation leads mainly to nearest-neighbor excitations between any pair of neighboring sites. For the standard BH Hamiltonian, the excitation spectrum consists of a large peak at energy $U$ (figure 6(d)). The nonzero width of the peak is due to variations of the local chemical potential (the presence of a trap): the shallower the trap, the narrower the peak.

For the EMO Hamiltonian, the particle-hole excitations from different Mott plateaus have different mean excitation energies. The shift with respect to the mean value is determined by the $\Delta\mu_{i,j}$. For shallow traps $\Delta\mu_{i,j} \approx 2(j - i)\kappa$ ($i - i_0$), and the excitation spectrum from Mott plateau $n$ consists of two bands—one corresponding to inward (I) hopping, $\Delta\mu_{i,j} < 0$, the other one to outward (O) hopping, $\Delta\mu_{i,j} > 0$, with respect to the trap center. For the central Mott plateau with density $n_{max}$ a single, broad peak in the excitation spectrum emerges. This is clearly visible in figures 6(b) and (c). Two cases have been studied: a system of $N = 260$ particles in a trap with $\kappa = 0.001$ (b), moderate case) and $N = 700$, $\kappa = 0.0001$ (c), very shallow trap). Both cases were studied for the 1D optical lattice with $s = 15$, $s_\perp = 40$, $\lambda = 830$ nm, $a = 5.1$ nm. The tunneling $J_{n_i, n_j}$ is artificially set to 0 (the deep Mott regime).

A similar but smaller system is analyzed in a subsequent numerical study, for the same parameters taking fully into account the tunneling effects. We choose a much tighter trap with curvature, $\kappa = 0.009$, and use true values of the hopping constants $J_{n_i, n_j}$. The Wannier function calculations give $U_{BH} = 0.662E_R$, and the renormalization procedure gives $U_2 \approx 0.565E_R$, $W \approx 0.0125E_R$.

We fill the trap with $N = 36$ particles. This system is similar to the one studied in [32, 33]. The system states are represented by MPS vectors and evolved using the TEBD algorithm [25, 26]. The ground state of the system, being the initial state for the evolution, is calculated
<figure><img src="image_8.png" /></figure>

**Figure 6.** Effects of higher Bloch bands on absorption spectroscopy in the deep Mott (low *J*) regime, *s* = 15, *s*⊥ = 40. Panel (a) shows the well-known wedding cake structure with *n* = 1, 2 and 3 Mott plateaus. Excitations within each plateau (colored respectively light gray, dark gray and black, for *n* = 1, 2 and 3) have energies depending on the Mott plateau density and the trapping potential. Inward and outward hopping lead to a splitting of the absorption structure, a partial splitting for a moderate harmonic trap ((b), *κ* = 0.001) or a broad well-resolved structure for a shallow trap ((c), *κ* = 0.0001) in contrast to the standard BH case (d).

using an imaginary time evolution with bond dimension $\chi = 50$. The local Hilbert space assumes maximal filling of six bosons per site.

The density profiles of ground states of the BH and EMO models are practically the same with a central plateau of two particles per site. Thus, any change of excitation frequencies can be interpreted as an effect of coupling constants renormalization. We have performed the absorption spectroscopy simulation for time $t = 100\hbar/E_R$. The modulation amplitude of the lattice was $s_m = 1$. The results are presented in figure 7.

The major difference between the spectra obtained for the BH and EMO Hamiltonians is a significant shift of the observed structures. While for the BH case, the main structure is centered at $U_{BH}$, it has a similar shape, but centered around $U_2$ in the EMO case. Because of the steep harmonic trap—thus the large changes in local chemical potentials—the structure of the peaks is rather complex. Note the global broadening for the EMO case, and an additional peak in the main $U_2$ structure, corresponding to the *n* = 1 plateau excitations having an energy larger by roughly $3W$, as discussed above.

A second small peak on the right appears at $E = 1.75U_2$ (the EMO case) and $E = 1.8U_{BH}$ (the BH case). It corresponds to a particle-hole excitation on the edge between the *n* = 1 and 2 Mott plateaus as identified in [33]. The right panel shows that the spectra become quite similar if rescaled by their proper energy scale, $U_{BH}$ or $U_2$.
<figure><img src="image_9.png" /></figure>

**Figure 7.** Absorption spectrum (modulation time $t = 100\hbar/E_R$; modulation amplitude $s_m = 1$). Black dashed lines correspond to the BH model and red solid curves to the EMO theory. The left panel shows spectra on a common energy scale; observe the significant shift of the EMO structure toward smaller energies. Bars above the plot give the mean expected positions of peaks for the $n = 2$ Mott plateau. The right panel shows the same data with the rescaled energy axes ($U_{BH}$ for the black curve and $U_2$ for the red one).

The absorption spectra are quite sensitive to the details of the system. Taking the same parameters for a slightly larger number of particles may create situations where the ground states of the BH and EMO Hamiltonians significantly differ. This is then reflected in the absorption spectra. If the density profile contains an $n = 3$ or higher plateau, then the structure of peaks may become more complicated, as discussed above.

We have also compared the absorption spectra in the SF regime. The lattice height is fixed at $s = 5$, $s_⊥ = 40$ [18]. The system is modulated for $t = 50\hbar/E_R$ with $s_m = 0.2$. The results are presented in figure 8. Unlike in the Mott regime, the positions of the absorption peaks are not determined solely by the interaction. In particular, no global shift of the structure is observed. In both cases, one observes a broad resonance around the recoil energy, with complicated detailed structures. Note that the modulation depth is much smaller than in the MI situation (to avoid significant excitation of the system) and therefore the absorbed energy per particle is much smaller than in figure 7.

## 4.2. *The Florence experiment revisited*

In the Florence experiment [18], the starting point is an ultra-cold gas in a harmonic trap (without optical lattice). The optical lattices are then ramped up (assuming $s(t) = 0.2s_⊥(t)$) with an exponential ramp $s(t) \sim s_0(1 - \exp(t/\tau))$, for $\tau = 30$ ms, and a total ramping duration 100 ms. The system soon becomes quasi-1D producing a set of 1D tubes. If a 'disordered' system is desired, an additional optical lattice, with a different wavelength $\lambda_2$, is superimposed along the tubes. This adds a potential $V_d(x) = s_2 \sin^2(k_2x)$. For $s_2 \ll s$, it acts effectively as the shift of the on-site energy, i.e. an additional pseudo-random disorder $\Delta V_i = s_2 \sin^2(\frac{\lambda}{\lambda_2}\pi i + \phi)$, where $\phi$ represents the offset between the two optical lattices. We take a generic value
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Absorption spectra in the SF case. The absorbed energy per particle is plotted as a function of the frequency of modulation of the lattice depth. Here $s = 5$, $s_m = 0.2$, $s_\perp = 40$, $U_{BH} \approx 0.465 E_R$. The renormalized interaction parameter: $U_2 \approx 0.406 E_R$. The left panel corresponds to 20 particles and the right panel to 36 particles.</figcaption></figure>

$\phi = 0.12345$. If the ratio $\lambda/\lambda_2$ is chosen irrational enough, the bichromatic lattice simulates a disorder well enough for a finite system [69–73]. The dependence of the pseudo-disorder amplitude $s_2$ on time is set by demanding that $s_2(t) \sim s(t)$ (all optical lattices are ramped up simultaneously). We will consider three cases: no disorder ($s_2(t) = 0$), weak disorder ($s_2 = \frac{1}{32}s$) and strong disorder ($s_2 = \frac{5}{32}s$).

Consider first the no disorder case $s_2 = 0$. As the initial state, we take 151 particles on 81 sites in the presence of a harmonic confinement coming both from the trap and the transverse laser profile. The detailed procedure using the TEBD algorithm is described in [29, 33].

After the optical lattice is ramped, absorption spectroscopy is performed for 30 ms (the conversion unit is $20.91\hbar/E_R = 1$ ms). In the recent numerical investigation of this experiment [29], a discrepancy between experimental [18] and numerical results was found. The reported position of the first absorption peak was 1.9 kHz [18], while Wannier function calculations gave 2.3 kHz [29]. The renormalization procedure renormalizes the value of the $U$ parameter to $U_2 = 2$, $U_3 = 1.96$ and $U_4 = 1.91$ kHz. This suggests that the positions of absorption peaks due to the $n = 1, 2$ and $3$ Mott plateaus are: $U_2 = 2$ kHz, $-2U_2 + 3U_3 = 1.85$ kHz, $U_2 - 6U_3 + 6U_4 = 1.74$ kHz. The ‘average’ peak position is 1.87 kHz. Therefore, the EMO Hamiltonian provides an estimate of the peak position in good agreement with the experiment.

A simulation of the absorption spectrum performed for this system confirms this finding as shown in figure 9. Although the initial state when the periodic lattice modulation starts is not the ground state, but a wavepacket dynamically created during the ramping of the lattice, the peak positions are well predicted by the EMO model. The positions of the first and second peaks agree quite well with the experiment (the relative height is different presumably because of the strong modulation used in [18]).

The exponential ramping of the optical lattice in the experiment [18] may not be adiabatic as discussed in [29] using a standard BH description. Due to the discrepancy in the position of the absorption peak, the lattice depth was adjusted in [29]. Instead of ramping the lattice up
<figure><img src="image_11.png" /></figure>

**Figure 9.** Absorption spectrum obtained by applying lattice modulation with amplitude $s_m = 1$ on the wavepacket created by exponential ramp up to $s = 16$. The black dashed line corresponds to the standard BH model and the red line is the result of EMO theory. The position of the absorption peaks in the latter case reproduces well the experimental results [18].

to $s = 16$, the final value $s = 14$ was considered. In some sense, such a simplified approach may be viewed as a renormalization of the BH parameters (without insight into its origin explained in section 4.1). Let us stress that the agreement between the experimental position of absorption peaks and the EMO predictions proves the necessity of using the EMO theory to explain quantitatively the experimental results.

With that modification of the final $s$ value, it was found using the BH model [29] that the overlap of the prepared wavepacket on the ground state at the final $s$ value was about 9% in the absence of disorder. It is most interesting to see how taking into account higher bands within EMO theory affects the adiabaticity of the dynamics. The simulation performed for the similar exponential ramp starting at $s = 5$ up to 16 yields an overlap of the dynamical wavepacket on the ground state at $s = 16$ equal to 17.3%. It may be qualitatively understood: in the EMO model, the effective interactions are weaker and the effective tunneling larger, allowing particles to redistribute more efficiently among sites during the ramp.

We have also tested an optimized $s(t)$ pulse shape as in [29]. By choosing $s(t)$ changing slowly close to the phase transition point for the $n = 3$ Mott lobe, we have been able to enhance adiabaticity up to 33% squared overlap with the ground state.

The presence of disorder has a devastating influence on adiabaticity, similarly to the standard BH case [29]. We have found that for a small disorder, $s_2 = \frac{1}{32}s$ the squared overlap is a fraction of a per cent (0.005), whereas for the strong disorder $s_2 = \frac{5}{32}s$ it becomes vanishingly small (of the order of $10^{-9}$, beyond the accuracy of the calculation).

**5. Conclusion**

The aim of this paper is twofold. In the first part, we have presented an efficient numerical implementation of the approach sketched in [39], which makes it possible to compute the
parameters of the effective Hamiltonian for bosons in optical lattices. The approach goes beyond the standard BH model [7] limited to the lowest Bloch band. The effective Hamiltonian approach which includes contributions from higher lying bands (MO approach) has been shown to lead to new effects and even new phases [37] for bosonic systems (see also [35]). Our scheme of perturbatively generated basis seems clearly superior to the energy-selected basis used in [39] for low and moderate occupation numbers and allows for better estimates of Hamiltonian parameters. These estimates may be extrapolated to an infinite number of bands.

We have applied the method not only to the standard 3D cubic lattice, but also to reduced 1D and 2D problems, where the lattice depth is different in various directions. The effective Hamiltonian obtained depends on the dimensionality of the problem. In effect, mean-field phase diagrams as obtained with the Gutzwiller ansatz differ even if they are rescaled by the lattice coordination number. It turns out that the role of excited bands is even more pronounced for reduced dimensionality problems than for a 3D lattice.

Motivated by this difference, we have investigated whether the dynamics is different in a standard BH model and for the EMO theory. We have considered two cases, the energy absorption created by modulation of the lattice height and the quasi-adiabatic passage from the SF to the MI phase. In both situations, it turns out that taking into account the density-dependent tunneling terms as well as modification of interactions may lead to significant differences between the two approaches. For the same lattice depth, the effective interactions turn out to be significantly weaker than in the standard BH model. This results in profound differences in the absorption spectra such as significant shifts of absorption peaks. Similarly, the full effective theory predicts that the transition from SF to MI is more adiabatic than with the standard BH model [29].

The results presented in this paper should have direct applicability to any experiment using bosons in an optical lattice, with multiple site occupations.

## Acknowledgments

We thank an anonymous referee for pointing out the potential problem of contact interaction potentials in 3D calculations. MŁ acknowledges communications with D-S Lühmann on the details of calculations in [39]. This work was supported by the International PhD Projects Programme of the Foundation for Polish Science within the European Regional Development Fund of the European Union, agreement no. MPD/2009/6. JZ acknowledges partial support from Polish National Center for Science under grant no. DEC-2012/04/A/ST2/00088.

## Appendix. Diagonalization of the on-site Hamiltonian

We describe the approach we use to generate a perturbatively based variational set used in the diagonalization of the on-site $n$ particle problem

$$ \mathcal{H}_{\text{loc}} = H_E + H_U = \sum_{\alpha} E_{\alpha} \hat{n}_{\alpha} + \sum_{\alpha \beta \gamma \delta} U^{\alpha \beta \gamma \delta} \hat{a}_{\alpha}^{\dagger} \hat{a}_{\beta}^{\dagger} \hat{a}_{\gamma} \hat{a}_{\delta}. \quad (A.1) $$

The most elementary possible excitation promotes a boson from the $\alpha = 0$ to 2 band. This defines two limits: $n^{3/2}U^{2000} \ll E_2 - E_0$ and $n^{3/2}U^{2000} \gg E_2 - E_0$. $U^{\alpha\beta\gamma\delta}$ becomes smaller as the $\alpha, \beta, \gamma, \delta$ indices grow.
The first regime corresponds to the ordinary BH model in which $H_U$ can be treated as a small perturbation, with zero-order ground state $|\psi_0\rangle = (\hat{a}^{\dagger}_{(0,0,0)})^n|0\rangle$, just as was assumed for the ordinary BH Hamiltonian. The large density regime is dominated by $H_U$. A straightforward application of the multiband model in that regime is not justified; however, a mean-field-based renormalization scheme leads to the approximate model of the same form with modified parameters [34]. The transition to the nonperturbative regime occurs at approximately 10–15 atoms per site.

Consider the low-energy regime. By vectors reachable in $k$th-order perturbative expansion, we call those Fock states $|\psi_p\rangle$ for which $\langle\psi_0|H_U^k|\psi_p\rangle \neq 0$. In particular, a full basis can be generated with order $\lceil \frac{n}{2} \rceil$. Let us denote by $B_k$ the set of vectors reachable in $k$th order and unreachable in $(k-1)$th order. The full variational basis of Fock states is $B = \cup_k B_k$ with $B_k$ pairwise disjoint, and $B_k = \emptyset$, for $k > \lceil \frac{n}{2} \rceil$.

For numerical diagonalization, a proper, not too large $S \subset B$ has to be chosen. In [39], $S$ consists of vectors $|\psi\rangle$ with the least $E_{|\psi\rangle} = \langle\psi|\mathcal{H}_{\text{loc}}|\psi\rangle$. In the first-order perturbation method, we choose vectors from $B_1$ with the largest values of

$$f_1(\psi) = \ln \frac{|\langle\psi_0|\hat{H}_U|\psi\rangle|}{E_{|\psi\rangle} - E_{|\psi_0\rangle}}. \qquad (A.2)$$

Disregarding vectors from $\cup_{k \ge 2} B_k$ seems to be a crude approach, but still our numerical calculations prove this ‘basic’ perturbative approach to be more efficient for low densities than the least-energy-based selection. The perturbation theory provides a way to evaluate a perturbative contribution of vectors from $B_k$, for arbitrary $k$, which, unfortunately, is computationally involved (summation over intermediate states, degeneracy resolution). As we do not calculate the ground-state energy within the perturbation theory treatment, but only motivate the choice of a variational basis for numerical diagonalization, the following function,

$$f_2(\psi) = \ln \sup_{|\psi_1\rangle \in B_1} \frac{|\langle\psi_0|\hat{H}_U|\psi_1\rangle \langle\psi_1|\hat{H}_U|\psi\rangle|}{(E_{|\psi\rangle} - E_{|\psi_0\rangle})(E_{|\psi_1\rangle} - E_{|\psi_0\rangle})}, \qquad (A.3)$$

is chosen to approximate the relevancy measure for vector $|\psi\rangle$. It minimizes the state’s $|\psi\rangle$ energy and maximizes the overlap.

A detailed comparison of the three methods: the least-energy, first-order perturbative and ‘improved’ first-order perturbative approaches are presented in figure A.1. We have used $K = 10000$ vectors (plus additional second higher order vectors for the ‘improved’ method), set $s = s_\perp = 38$ and compared the three methods as a function of $n$. Two regimes, perturbative and nonperturbative, emerge, as expected. In the small and moderate $n$ regime, the perturbative approach gives a better estimate for the ground-state energy than the least-energy method [39]. The opposite tendency is visible in the nonperturbative regime.

Let us describe in detail how a choice of the basis with the largest $f_1$ and $f_2$ values is performed. We have fixed the maximal number of Bloch bands included at $B = 15$. We use the Markov chain Monte Carlo method, which is quite general and can be applied at any order $k$ of perturbation expansion. $B_k$ is, in general, too large to evaluate function $f_k$ for all elements (its size increases exponentially with the total number of particles). It is usually possible just to scan the whole $B_1$ set, so from this point on we assume that $k \ge 2$. To choose $K$ vectors with the largest $f_k$ values, we construct a random walk based on Metropolis’ Monte Carlo algorithm. A state of the random walk is a finite $k+1$-tuple of $n$-particle Fock states
<figure><img src="image_12.png" /><figcaption>**Figure A.1.** Comparison of $U_n$ obtained from diagonalization of $\mathcal{H}_{\text{loc}}$, equation (A.1), with energy selected basis and with perturbatively chosen sets. Each perturbative set consists of 10 000 vectors of first order and of 0, 100 and 800 vectors of second order (curves: black, red and green). The blue dashed curve denotes the least-energy basis result. For small and moderate $n$, the perturbative based basis is clearly superior to the least-energy set. The failure of the perturbation theory approach for $n$ sufficiently large is apparent, too.</figcaption></figure>

$\mathcal{V} = (|\psi_0\rangle, |\psi_1\rangle, |\psi_2\rangle, \dots, |\psi_{k_1}\rangle, |\psi\rangle)$, for $|\psi_i\rangle \in \mathcal{B}_i$. For any such $\mathcal{V}$, we define the following generalization of (A.3):

$$ 
\begin{align} 
g_k(\mathcal{V}) &= \ln \frac{|\langle\psi_0|\hat{H}_U|\psi_1\rangle\langle\psi_1|\hat{H}_U|\psi_2\rangle|, \dots, \langle\psi_{k-1}|\hat{H}_U|\psi\rangle|}{(E_{|\psi\rangle} - E_{|\psi_0\rangle})(E_{|\psi_{k-1}\rangle} - E_{|\psi_0\rangle}), \dots, (E_{|\psi_1\rangle} - E_{|\psi_0\rangle})}, \tag{A.4} \\ 
f_k(\psi) &= \sup_{|\psi_1\rangle \in \mathcal{B}_1, \dots, |\psi_{k-1}\rangle \in \mathcal{B}_{k-1}} g_k(\mathcal{V}). \nonumber 
\end{align} 
$$

To obtain the random walk, we have to update $\mathcal{V}$. First we choose at random a Fock state $|\psi_l\rangle \in \mathcal{V}$ to be updated. With equal probability, we update one or two particles of $|\psi_l\rangle$ preserving the total parity of the state. One-particle update is done according to $|\psi_l\rangle \rightarrow \hat{a}_{(i_x, i_y, i_z)}\hat{a}^{\dagger}_{(i_x, i'_y, i'_z)}|\psi_l\rangle$, $i_y \equiv i'_y (\text{mod } 2)$, whereas two-particle update is $|\psi_l\rangle \rightarrow \hat{a}_{(i_x, i_y, i_z)}\hat{a}^{\dagger}_{(i_x, i'_y, i'_z)}|\psi_l\rangle$, $i_y + j_y \equiv i'_y + j'_y (\text{mod } 2)$. All vectors are normalized. The direction $y$ is not special in any way: with equal probability any of $x, y, z$ is chosen. After the update a proposition $\mathcal{V}'$ is prepared. We automatically reject updates for which $|\psi_l\rangle \notin \mathcal{B}_l$. If that is not the case, the acceptance probability is determined as in Metropolis algorithm: it is given by $\min\{1, \exp[\beta(g_k(\mathcal{V}) - g_k(\mathcal{V}'))]\}$. The inverse temperature $\beta$ is tuned to optimize the sampling efficiency—we choose it by requiring the acceptance rate to be close to 0.3. After a successful update, the last element of the tuple $\mathcal{V}'$, state $|\psi_k\rangle$ is accepted into the solution set if its perturbative importance $g_k(\mathcal{V})$ is in the $K$ lowest values recorded so far. The accepted vector $|\psi_k\rangle$ is memorized as well as the importance value $g_k(\mathcal{V}')$. If $|\psi_k\rangle$ had been generated before, the memorized value of $g_k$ is updated (only if the new value is larger than the old one). If, in a subsequent few thousand sweeps (empirical value), no vector makes it into the solution set, nor are $g_k$ values updated, then the procedure is restarted. This ensures that all low-energy
excitations are taken into account (the starting point is always the low-energy configuration). Altogether, we make $2 \times 10^9$ Monte Carlo sweeps to generate basis of size 40 000 (as used for the results presented in the main part of the paper).

If all Bloch bands were included, then the set of Fock space would be infinite. On the other hand, only a finite number of them could satisfy the inequality: $f_k > \varepsilon$. The values of $f_i$ for the remaining states are very close to 0, and a singularity in density of states $\frac{\partial f_k}{\partial n}$ arises. Logarithm is used to 'smooth' this singularity for numerical purposes. It does not affect the ordering, as ln is increasing injection.

## References

[1] Lewenstein M, Sanpera A, Ahufinger V, Damski B, Sen(De) A and Sen U 2007 *Adv. Phys.* **56** 243–379

[2] Yukalov V 2009 *Laser Phys.* **19** 1–110

[3] Lewenstein M, Sanpera A and Ahufinger V 2012 *Ultracold Atoms in Optical Lattices: Simulating Quantum Many-Body Systems* (Oxford: Oxford University Press)

[4] Greiner M, Mandel O, Esslinger T, Hansch T W and Bloch I 2002 *Nature* **415** 39–44

[5] Paredes B, Widera A, Murg V, Mandel O, Folling S, Cirac I, Shlyapnikov G V, Hansch T W and Bloch I 2004 *Nature* **429** 277–81

[6] Bloch I 2005 *Nature Phys.* **1** 23–30

[7] Jaksch D, Bruder C, Cirac J I, Gardiner C W and Zoller P 1998 *Phys. Rev. Lett.* **81** 3108–11

[8] Gersch H A and Knollman G C 1963 *Phys. Rev.* **129** 959–67

[9] Fisher M P A, Weichman P B, Grinstein G and Fisher D S 1989 *Phys. Rev. B* **40** 546–70

[10] van Oosten D, van der Straten P and Stoof H T C 2001 *Phys. Rev. A* **63** 053601

[11] Freericks J K and Monien H 1996 *Phys. Rev. B* **53** 2691–700

[12] Damski B and Zakrzewski J 2006 *Phys. Rev. A* **74** 043609

[13] Verstraete F and Cirac J I 2010 *Phys. Rev. Lett.* **104** 190405

[14] Rokhsar D S and Kotliar B G 1991 *Phys. Rev. B* **44** 10328–32

[15] Zakrzewski J 2005 *Phys. Rev. A* **72** 039904

[16] Stöferle T, Moritz H, Schori C, Köhl M and Esslinger T 2004 *Phys. Rev. Lett.* **92** 130403

[17] Kinoshita T, Wenger T and Weiss D S 2004 *Science* **305** 1125–8

[18] Fallani L, Lye J E, Guarrera V, Fort C and Inguscio M 2007 *Phys. Rev. Lett.* **98** 130404

[19] Gemelke N, Zhang X, Hung C L and Chin C 2009 *Nature* **460** 995–8

[20] Bloch I, Dalibard J and Nascimbene S 2012 *Nature Phys.* **8** 267–76

[21] Trotzky S, Chen Y A, Flesch A, McCulloch I P, Schollwock U, Eisert J and Bloch I 2012 *Nature Phys.* **8** 325

[22] Schollwöck U 2011 *Ann. Phys.* **326** 96–192

[23] White S R and Feiguin A E 2004 *Phys. Rev. Lett.* **93** 076401

[24] Eisert J, Cramer M and Plenio M B 2010 *Rev. Mod. Phys.* **82** 277–306

[25] Vidal G 2003 *Phys. Rev. Lett.* **91** 147902

[26] Vidal G 2004 *Phys. Rev. Lett.* **93** 040502

[27] Kollath C, Iucci A, Giamarchi T, Hofstetter W and Schollwöck U 2006 *Phys. Rev. Lett.* **97** 050402

[28] Clark S R and Jaksch D 2006 *New J. Phys.* **8** 160

[29] Zakrzewski J and Delande D 2009 *Phys. Rev. A* **80** 013602

[30] Cucchietti F M, Damski B, Dziarmaga J and Zurek W H 2007 *Phys. Rev. A* **75** 023603

[31] Eurich N, Eckstein M and Werner P 2011 *Phys. Rev. B* **83** 155122

[32] Huo J W, Zhang F C, Chen W, Troyer M and Schollwöck U 2011 *Phys. Rev. A* **84** 043608

[33] Łącki M, Delande D and Zakrzewski J 2012 *Phys. Rev. A* **86** 013602

[34] van Oosten D, van der Straten P and Stoof H T C 2003 *Phys. Rev. A* **67** 033606

[35] Johnson P R, Tiesinga E, Porto J V and Williams C J 2009 *New J. Phys.* **11** 093022
[36] Larson J, Collin A and Martikainen J P 2009 *Phys. Rev. A* **79** 033603

[37] Dutta O, Eckardt A, Hauke P, Malomed B and Lewenstein M 2011 *New J. Phys.* **13** 023019

[38] Sakmann K, Streltsov A I, Alon O E and Cederbaum L S 2011 *New J. Phys.* **13** 043003

[39] Lühmann D S, Jurgensen O and Sengstock K 2012 *New J. Phys.* **14** 033021

[40] Mering A and Fleischhauer M 2011 *Phys. Rev. A* **83** 063630

[41] Ospelkaus S, Ospelkaus C, Wille O, Succo M, Ernst P, Sengstock K and Bongs K 2006 *Phys. Rev. Lett.* **96** 180403

[42] Will S, Best T, Schneider U, Hackermüller L, Lühmann D S and Bloch I 2010 *Nature* **465** 197–201

[43] Bissbort U, Deuretzbacher F and Hofstetter W 2011 arXiv:1108.6047

[44] Bloch I, Dalibard J and Zwerger W 2008 *Rev. Mod. Phys.* **80** 885–964

[45] Fedichev P O, Bijlsma M J and Zoller P 2004 *Phys. Rev. Lett.* **92** 080401

[46] Dickerscheid D B M, Al Khawaja U, van Oosten D and Stoof H T C 2005 *Phys. Rev. A* **71** 043604

[47] Gubbels K B, Dickerscheid D B M and Stoof H T C 2006 *New J. Phys.* **8** 151

[48] Diener R B and Ho T L 2006 *Phys. Rev. Lett.* **96** 010402

[49] Wouters M and Orso G 2006 *Phys. Rev. A* **73** 012707

[50] Büchler H P 2010 *Phys. Rev. Lett.* **104** 090402

[51] Cui X, Wang Y and Zhou F 2010 *Phys. Rev. Lett.* **104** 153201

[52] von Stecher J, Gurarie V, Radzihovsky L and Rey A M 2011 *Phys. Rev. Lett.* **106** 235301

[53] Mark M J, Haller E, Lauber K, Danzl J G, Daley A J and Nägerl H-C 2011 *Phys. Rev. Lett.* **107** 175301

[54] Albeverio S, Gesztesy F and Hoegh-Krohn R 1988 *Solvable Models in Quantum Mechanics* (Berlin: Springer)

[55] Wódkiewicz K 1991 *Phys. Rev. A* **43** 68–76

[56] Fermi E 1936 *Ric. Sci.* **7** 13

[57] Busch T, Englert B G, Rzążewski K and Wilkens M 1998 *Found. Phys.* **28** 549–59

[58] Rontani M, Åberg S and Reimann S M 2008 arXiv:0810.4305

[59] Evenbly G and Vidal G 2009 *Phys. Rev. Lett.* **102** 180406

[60] Verstraete F, Wolf M M, Perez-Garcia D and Cirac J I 2006 *Phys. Rev. Lett.* **96** 220601

[61] Verstraete F and Cirac J I 2004 arXiv:cond-mat/0407066

[62] Vidal G 2007 *Phys. Rev. Lett.* **99** 220405

[63] Sowiński T, Dutta O, Hauke P, Tagliacozzo L and Lewenstein M 2012 *Phys. Rev. Lett.* **108** 115301

[64] Freericks J K and Monien H 1994 *Europhys. Lett.* **26** 545

[65] Kühner T D and Monien H 1998 *Phys. Rev. B* **58** R14741–4

[66] Pilati S and Troyer M 2012 *Phys. Rev. Lett.* **108** 155301

[67] Zakrzewski J and Delande D 2008 *AIP Conf. Proc.* **1076** 292–300

[68] Łącki M, Delande D and Zakrzewski J 2011 *Acta Phys. Pol. A* **12-a** A-178

[69] Diener R B, Georgakis G A, Zhong J, Raizen M and Niu Q 2001 *Phys. Rev. A* **64** 033416

[70] Damski B, Zakrzewski J, Santos L, Zoller P and Lewenstein M 2003 *Phys. Rev. Lett.* **91** 080403

[71] Roth R and Burnett K 2003 *Phys. Rev. A* **68** 023604

[72] Roscilde T 2008 *Phys. Rev. A* **77** 063605

[73] Roux G, Barthel T, McCulloch I P, Kollath C, Schollwöck U and Giamarchi T 2008 *Phys. Rev. A* **78** 023628