Dynamics of cold bosons in optical lattices: Effects of higher Bloch bands

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 with 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 energy absorption rate when the lattice depth is periodically modulated.


Introduction
Ultracold bosonic atom gases in optical lattices have been an ultrahot research area in recent years (for recent reviews and an extensive reference list see [1][2][3]). They provide means to create and control experimental systems mimicking different condensed matter physics models [4][5][6]. The interest has been stimulated in part by the fact that there exists [7] an accurate mapping of 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][10][11][12] and numerical analysis, though promising new ideas were recently proposed 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 the area of extensive experimental research [16][17][18][19][20][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 enable simulations targeting real-life many-body systems [23], with controllable error and no systematic errors due to unsounded 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 it possible the construction of an efficient variational set -the so called Matrix Product States (MPS). Several numerical investigations of experimental and 'close to experimental' systems have been performed [27][28][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 to the Mott insulator 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, number of peaks in the spectrum and their heights are related to the state of the gas: either superfluid 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 10 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 been also studied for lower densities [35][36][37][38] in relation to the discovery that higher bands affect the superfluid-insulator transition in Bose-Fermi mixtures [39,40]. This gave explanation for experimental observations of the shift of the Superfluid (SF) -Mott insulator (MI) phase transition [41] which could not be explained by single-band approaches. In [36], the first excited band is included in a twoflavour model; effects of higher bands could also be built in via an effective three-body interaction in the lowest band [35,37].
Recently, another approach for 3D optical lattices has been proposed [39,42,43] (which somehow resemble 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 number of particles are used as a local Hilbert space, replacing the usual Fock basis |n . 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 of BH Hamiltonian parameters is large in the MI regime, where the energy gap is also large, contrary to a 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 onsite 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.

Tight-binding descriptions of an ultracold boson gas in an optical lattice
The second quantization Hamiltonian for a dilute gas of interacting bosonic atoms in the optical lattice potential V ( r) and external trapping potential V e ( r) is of the form: where V int ( r, r ′ ) is an isotropic short-range pseudopotential modelling s-wave interactions [44] with a s being the scattering length. The Hamiltonian admits a natural energy scale, set by the recoil energy E R =h 2 k 2 2m , where k = 2π λ and λ 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 ≪ s y = s z =: s ⊥ , then a classical setup for a 1D optical lattice is obtained [16]. Tunneling in y, 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 ≪ s z =: s ⊥ a 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 three-dimensional ones, the labels "1D" or "2D" just refering to situations where s and s ⊥ are very different.
Jaksch and Zoller in their seminal paper [7] 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 ( r i )): The field operator is expanded in the set of the lowest band Wannier functions of the lattice: ψ( r) = w 0 i ( 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: where i and j are neighbouring sites) and U = 4πh 2 as m w 0 i ( r) 4 d 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 in 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][46][47][48][49][50][51][52]. 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.2nm, much shorter than the lattice spacing λ/2 = 377nm.
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 ≪ 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 real world.
It is well known that, in dimension higher than 1, contact interactions cannot be modelled by a δ-potential, but require a specific regularization to avoid artificial divergences. Mathematically, a self-adjoint extension of the original Hamiltonian is needed [53]. It boils down to a dimension-dependent regularization of the δpotential [54], which, in 3D, is the so-called Fermi pseudo-potential [55] used in Eq. (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 [56] 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 [57]. 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 → ∞, 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 Bose-Hubbard Hamiltonian: and Note that, because Wannier functions are smooth, the potential (2 may be replaced by a contact Fermi potential in the integral (5). By construction J αβ ij = 0 for α = β. 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 neighbours hopping may be necessary -see [21]).
A 3D Wannier function being a product of 1D Wannier functions, the 3D integral in U αβγδ ijkl is a product of 3 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 onsite dimension d of the n-particle problem restricted to the lowest B bands is d = B 3 +n−1 n , (the onsite problem is genuinely 3D even for a quasi-1D models) 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 MERA, PEPS [58][59][60][61]) 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 = 0 for i = j = k = l together with considering the lowest Bloch band only (α = β = γ = δ = 0) [7] leads directly to the Bose-Hubbard Hamiltonian, Eq. (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: (withn α =â † αâ α ) can be diagonalized to yield a space of n particle ground states. The eigenenergies ǫ n 0 of the on-site n particle ground states |ψ n 0 , are the starting point in determining new values of U parameters in the effective Hamiltonian. To define renormalized values of U, the energy ǫ n 0 has to be decomposed into the interaction energy [which in case of the BH model, is just U 2 n(n − 1)] and a single-particle energy (which in the BH case shifts µ by the lowest Bloch band energy). The most natural way to define the interaction energy would be to use: Unfortunately, U n cannot be defined in such a way if we request the Hamiltonian to have a form resembling Eq. (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: 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 0000 iiij induces an effective coupling, which is proportional to the density, that is a term U 0000 iiij a i a † j (n i + n j − 1) + h.c. in the effective Hamiltonian (called bond-charge term in [39,62]). 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 0000 The correction is at most of the order of the raw tunneling. When higher bands are taken into account, the modification of |ψ n 0 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: where P i n = |i, n 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.

Solving the onsite 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 ∆ = E 2 − E 0 . Promoting a particle from the lowest to the second excited 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 ≪ ∆. For typical parameters corresponding to Rb scattering length and s of the order of 30, ∆/U 2000 ≈ 60 yielding the limiting value of n ≪ 15.
For small occupation numbers, a perturbative approach seems justified. Perturbation theory enables us to estimate the impact each excited vector |ψ p has on the ground state energy. Let |ψ 0 be an n-boson ground (Fock) state of H E . The larger the matrix element | ψ 0 |H loc |ψ p | 2 and the smaller the energy ψ p |H loc |ψ p the larger is the impact. The perturbative scheme provides a hint on 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 least energy principle -with minimal ψ p |H loc |ψ p .
We have performed a detailed analysis comparing both methods for s = s ⊥ = 34.8, a strongly coupled case (3D optical lattice, 2a s /λ = 0.014, λ = 754 nm -parameters taken as typical values from [39]). We choose a system with n = 2 − 5 particles and 40000 basis vectors according to both least-energy (as in [39]) and the perturbative method (39900 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 ≈ 7 − 9 bands are included. The false saturation occurs because the least-energy method does not evaluate | ψ 0 |H loc |ψ p | 2 and therefore fills the variational basis with low-energy 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 +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 → ∞ 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 ⊥ is still large.
Inspection of Fig. 1 and Fig. 2 shows that the difference between consecutive U n is approximately constant, i.e., 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 (Eq. (12) of [35]) and W is simply related to their parameterŨ 3 . The deviations from the linear form, eq. (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

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 µ and the Hamiltonian parameters: J and U. A MI phase is determined by integer density n and noncompressibility: ∂ n ∂µ = 0. The rest of the phase diagram is the compressible SF phase [9,11,63,64]. We perform now the mean field analysis of the phase diagrams of Hamiltonians (3) and (10). A Gutzwiller analysis of the Bose-Hubbard model is a variational minimization of the following functional (i.e. mean ground state energy): using the Gutzwiller ansatz: |ψ = |ψ l with |ψ l = f n |n the on-site wavefunction. The influence of the lattice geometry is reduced only to the coordinate number, z, in the first term, as â i = â j due to the translational invariance in thermodynamic limit.
The phase diagram depends only on the single parameter 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: Clearly a single parameter is no longer sufficient to describe the mean field problem. Let us denote by J BH the tunneling amplitude and by U BH the interaction of the standard BH 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 , µ) coordinate space used for a traditional phase diagram. This mapping allows also to directly compare the results obtained using Eqs. (12) and (13). Figure 4 shows the Gutzwiller phase diagrams for 1D, 2D, 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 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 Bose-Hubbard 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 ⊥ . We show just a single choice for a generic value of s ⊥ = 34.8.
For 3D, a comparison with the mean field diagram obtained for 9 bands [39] is possible. The difference is quite small, the difference in U n (of the order of few %see Fig. 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 has been reported for various scattering lengths [65] and compared with the mean field results of [39]. We consider a single ratio of a s /λ as appropriate for 87 Rb, the small difference between our results and that of [39] for unit filling on the scale of Fig. 1 of [65] is negligible.
In 1D, the mean field approximation is inaccurate. To get 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, 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 [61,66]. The transverse lattice height is again s ⊥ = 34.8. Let us denote by E(N, s) the ground state energy of a N-particle system for lattice height s. We calculate approximations to the critical values of chemical potentials delimiting a Mott insulator region with average filling n by µ + (s) ≈ E(nL + 1, s) − E(nL, s), µ − (s) ≈ 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.

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 dependance of the energy absorption rate with 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 ωt, s m /s 0 ≪ 1. The simulation is performed in the presence of an harmonic trap V i = κ(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  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, blue, green, and red curves denote 1D, 2D and 3D lattices of the EMO Hamiltonian. Dashed red lines show the result obtained for 9 bands as in [39]. The limit zJ BH (s)/U BH (s), s → ∞ is different for each dimension. The s → ∞ 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 ⊥ = 34.8, λ = 754 nm, 2a s /λ = 0.014 as appropriate for 87 Rb [39]. 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 → 0). States excited during the modulation are mainly those that differ from the ground state by nearest-neighbours transfer of one particle [32,67]. Up to a small correction due to the difference of the local chemical potential µ i = µ−V i , the excitation energy is determined by the occupation numbers of the source site i and destination site j. It is: with ∆µ ij = µ i − µ j = V j − V i . Nearest neighbours excitation means that |i − j| = 1, and |n i − n j | ≤ 1. If U n = U (BH Hamiltonian), we have that ∆E(n i ; n j ) = (n j − n i + 1)U + ∆µ ij . By virtue of Eq. (11) we may approximate Eq. (14) by: ∆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:  A qualitative comparison of the expected absorption spectra for 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 6a). Weak, periodic modulation leads mainly to nearest neighbour excitations between any pair or neighbouring sites. For the standard BH Hamiltonian, the excitation spectrum consists of a large peak at energy U (Figure 6d). 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 ∆µ i,j . For shallow traps ∆µ i,j ≈ 2(j − i)κ(i − i 0 ) and the excitation spectrum from Mott plateau n consists in two bands -one corresponding to Inward (I) hopping, ∆µ i,j < 0, the other one to Outward (O) hopping, ∆µ i,j > 0, with respect to the trap centre. same parameters taking fully into account the tunneling effects. We choose a much tighter trap with curvature, κ = 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 ≈ 0.565E R , W ≈ 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 using an imaginary time evolution with bond dimension χ = 50. The local Hilbert space assumes maximal filling of 6 bosons per site.
The density profiles of ground states of the BH and EMO models are practically the same with a central plateau of 2 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 = 100h/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 (EMO case) and E = 1.8U BH (BH case). It corresponds to a particle-hole excitation on the edge between the n = 1 and n = 2 Mott plateaus as identified in [33]. The right panel shows that the spectra becomes quite similar if rescaled by their proper energy scale, U BH or U 2 .
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 a 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 superfluid regime. The lattice height is fixed at s = 5, s ⊥ = 40 [18]. The system is modulated for t = 50h/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 Mott insulator situation (to avoid significant excitation of the system) and therefore the absorbed energy per particle is much smaller than in Fig. 7.

Florence experiment revisited
In the Florence experiment [18], the starting point is a 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) ∼ s 0 (1 − exp(t/τ )), for τ = 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 different wavelength λ 2 , is superimposed along the tubes. This adds a potential V d (x) = s 2 sin 2 (k 2 x). For s 2 ≪ s, it acts effectively as the shift of the on-site energy, i.e. an additional pseudo-random disorder ∆V i = s 2 sin 2 λ λ 2 πi + φ , where φ represents the offset between the two optical lattices. We take a generic value φ = 0.12345. If the ratio λ/λ 2 is chosen irrational enough, the bichromatic lattice simulates a disorder well enough for a finite system [68][69][70][71][72]. The dependence of the pseudo-disorder amplitude s 2 on time is set by demanding that s 2 (t) ∼ s(t) (all optical lattices are ramped up simultaneously). We will consider 3 cases: no disorder (s 2 (t) = 0), weak disorder (s 2 = 1 32 s), and strong disorder (s 2 = 5 32 s). Consider first the no disorder case s 2 = 0. As 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 30ms (the conversion unit is 20.91h/E R = 1ms). 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.3kHz [29]. The renormalization procedure renormalizes the value of the U parameter to U 2 = 2 kHz, U 3 = 1.96 kHz, U 4 = 1.91 kHz. This suggests that the positions of absorption peaks due to the n = 1, 2, 3 Mott plateaus are: U 2 = 2kHz, −2U 2 + 3U 3 = 1.85kHz, U 2 − 6U 3 + 6U 4 = 1.74kHz. 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 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 position 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 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 prove the necessity of using the effective multiorbital 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 effective multiorbital theory affects the adiabaticity of the dynamics. The simulation performed for the similar exponential ramp starting at s = 5 up to s = 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 devastating influence on adiabaticity, similarly to the standard BH case [29]. We have found that for a small disorder, s 2 = 1 32 s the squared overlap is a fraction of percent (0.005) while for the strong disorder s 2 = 5 32 s it becomes vanishingly small (of the order of 10 −9 , beyond the accuracy of the calculation).

Conclusion
The aim of this paper is two-fold. 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 Bose-Hubbard model [7] limited to the lowest Bloch band. The effective Hamiltonian approach which includes contributions from higher lying bands (multiorbital 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 Bose-Hubbard model and for the effective multiorbital theory. We have considered two cases, the energy absorption created by modulation of the lattice height and the quasi-adiabatic passage from the superfluid to the Mott insulator 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 two approaches. For the same lattice depth, the effective interactions turn out to be significantly weaker than in the standard Bose-Hubbard 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 superfluid to Mott insulator is more adiabatic than with the standard Bose-Hubbard model [29].
The results presented in the present paper should have a direct applicability to any experiment using bosons in an optical lattice, with multiple site occupations. For small and moderate n, the perturbative based basis is clearly superior over the least-energy set. The failure of perturbation theory approach for n sufficiently large is apparent too.
Disregarding vectors from k≥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: is chosen to approximate the relevancy measure for vector |ψ . It minimizes the state's |ψ energy and maximizes the overlap. Detailed comparison of the three methods: the least-energy, first order perturbative and "improved" first order perturbative approach is presented in Figure 1. We have used K = 10000 vectors (plus additional second higher order vectors for the "improved" method), set s = s ⊥ = 38 and compared the three methods as a function of n. Two regimes: perturbative and not perturbative emerge, as expected. In the small and moderate n regime, the perturbative approach gives 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 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 k ≥ 2. To choose K vectors with 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 V = (|ψ 0 , |ψ 1 , |ψ 2 , . . . , |ψ k 1 , |ψ ), for |ψ i ∈ B i . For any such V we define the following generalization of (A.3): To get the random walk, we have to update V. First we choose at random a Fock state |ψ l ∈ V to be updated. With equal probability, we update one or two particles of |ψ l preserving the total parity of the state. One particle update is done according to |ψ l →â (ix,iy,iz)â † (ix,i ′ y ,iz) |ψ l , i y ≡ i ′ y (mod 2), while two particle-update is: |ψ l →â (ix,iy,iz)â(jx,jy,jz)â † (ix,i ′ y ,iz)â † (jx,j ′ y ,jz) |ψ l , i y + j y ≡ i ′ y + j ′ y (mod 2). All vectors are normalized. Direction y is not special in any way: with equal probability any of x, y, z is chosen. After the update a proposition V ′ is prepared. We automatically reject updates for which |ψ l ∈ B l . If that is not the case, the acceptance probability is determined as in Metropolis algorithm: it is given by min{1, exp[β(g k (V ′ ) − g k (V))]}. The inverse temperature β is tuned to optimize 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 V ′ , state |ψ k is accepted into the solution set if its perturbative importance g k (V) is in the K lowest values recorded so far. Th! e accepted vector |ψ k is memorized as well as the importance value g k (V ′ ). If |ψ k had been generated before, the memorized value of g k is updated (only if the new value of larger than the old one). If, in a subsequent few thousand sweeps (empirical value), no vector makes it into the solution set, nor g k values are 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 × 10 9 MC sweeps to generate basis of size 40000 (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 > ε. The values of f i for the remaining states are very close to 0, and a singularity in density of states ∂f k ∂ n arises. Logarithm is used to "smoothen" this singularity for numerical purposes. It does not affect the ordering, as ln is increasing injection.