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

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_BT 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〉. 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.

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

$$\begin{eqnarray} &&\fl H=\int \textrm{d}^3\vec r\ \Psi^\dagger(\vec r)\left(-\frac{\hbar^2}{2\,m } \nabla^2+V(\vec r)+V_e(\vec r)\right)\Psi(\vec r) \nonumber\\ &&+\frac{1}{2}\int\textrm{d}^3\vec r \,\textrm{d}^3\vec r' \Psi^\dagger(\vec r) \Psi^\dagger(\vec r')V_{\mathrm{int}}(\vec r, \vec r') \Psi(\vec r) \Psi(\vec r'), \end{eqnarray} \tag{ 1 }$$

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

$$\begin{equation} V_{\mathrm{int}}(\vec r,\vec r')=\frac{4\pi \hbar^2 a_{\mathrm{s}}}{m}\delta(\vec r-\vec r')\frac{\partial}{\partial |\vec r-\vec r'| }|\vec r-\vec r'| \end{equation} \tag{ 2 }$$

with a_s being the scattering length.

The Hamiltonian admits a natural energy scale, set by the recoil energy $E_{\mathrm {R}} = \frac {\hbar ^2 k^2}{2m},$ where $k=\frac {2\pi }{\lambda }$ 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²(kx) + s_y cos²(ky) + s_z cos²(kz). If s: = s_x ≪ 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 ≪ 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)$):

$$\begin{eqnarray} H_{\mathrm{BH}}&=&-J \sum\limits_{\langle i j \rangle} \hat a_i \hat a_j^\dagger +{\rm h.c.} + \frac{U}{2} \sum\limits_i \hat n_i(\hat n_i-1) - \sum\limits_i\hat n_i (\mu - V_i). \end{eqnarray} \tag{ 3 }$$

The field operator is expanded in the set of the lowest band Wannier functions of the lattice: $\psi (\vec r)=\sum \limits 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)\ \mathrm {d}\vec r$ (where i and j are neighboring sites) and $U= \frac {4\pi \hbar ^2 a_{\mathrm {s}}}{m} \int w_i^0(\vec r)^4 \,\mathrm {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 ⁸⁷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 ≪ 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 → ∞, 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

$$\begin{eqnarray} &&\fl H=-\sum\limits_{\alpha, i,j}J^{\alpha\alpha}_{ij} (\hat a^{\alpha}_{i})^\dagger \hat a^{\alpha}_{j} + \mathrm{h.c.} + \sum\limits_{{\stackrel{\scriptstyle\alpha\ldots\delta}{i\ldots l}}}U^{\alpha\beta\gamma\delta}_{ijkl} (\hat a^{\alpha}_i)^\dagger (\hat a^{\beta}_{j})^\dagger \hat a^{\gamma}_k \hat a^{\delta}_{l} + \sum\limits_{\alpha,i} (E_\alpha +V_i - \mu) (\hat a^{\alpha}_{i})^\dagger \hat a^{\alpha}_{i} \end{eqnarray} \tag{ 4 }$$

with

$$\begin{equation} U^{\alpha\beta\gamma\delta}_{ijkl}= \frac{4\pi \hbar^2 a\mathrm{_s}}{m} \int \mathrm{d}\vec r\ w^{\alpha}_{i} (\vec r) w^{{\beta}}_{j} (\vec r)w^{{\gamma}}_{k} (\vec r)w^{{\delta}}_{l} (\vec r) \end{equation} \tag{ 5 }$$

and

$$\begin{eqnarray*}&&J^{\alpha\alpha}_{ij} =-\int \mathrm{d}\vec r\ w^{\alpha}_{i}(\vec r) \left(-\frac{\hbar^2}{2m} \nabla^2+V(\vec r)\right)w^{\alpha}_{j}(\vec r).\end{eqnarray*}$$

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

$$\begin{equation} V_{\mathrm{c}}(\vec r,\vec r')=\frac{4\pi \hbar^2 a\mathrm{_s}}{m}\delta(\vec r-\vec r') \end{equation} \tag{ 6 }$$

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-neighbors 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 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=\left({}^{B^{3 }+n-1}_{\hphantom{B^{3 }n} n} \right)$ (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³). 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 ≠ 0 for i = j = k = l together with considering the lowest Bloch band only (α = β = γ = δ = 0) [7] leads directly to the BH Hamiltonian, equation (3), provided we chose the zero of the energy axis at E₀.

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

$$\begin{equation} \mathcal{H}_{\mathrm{loc}}= H_E + H_U = \sum\limits_\alpha E_\alpha \hat n_\alpha + \sum\limits_{\alpha\beta\gamma\delta} U^{\alpha\beta\gamma\delta} \hat a^\dagger_\alpha \hat a^\dagger_\beta \hat a_\gamma \hat a_\delta \end{equation} \tag{ 7 }$$

(with $\hat n_\alpha =\hat a^\dagger _\alpha \hat a_\alpha$) can be diagonalized to yield a space of n particle ground states. The eigenenergies ⁿ₀ of the on-site n particle ground states |ψⁿ₀〉 are the starting point in determining new values of U parameters in the effective Hamiltonian. To define renormalized values of U, the energy ⁿ₀ 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 μ by the lowest Bloch band energy). The most natural way to define the interaction energy would be to use

$$\begin{equation} \frac{U_n}{2} n(n-1) = \langle \psi_0^n |\sum\limits_{\alpha\beta\gamma\delta} U^{\alpha\beta\gamma\delta} \hat a^\dagger_\alpha \hat a^\dagger_\beta \hat a_\gamma \hat a_\delta |\psi_0^n\rangle. \end{equation} \tag{ 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

$$\begin{equation} \epsilon_0^n = \frac{U_n}{2} n(n-1) + n E_0 . \end{equation} \tag{ 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 induces an effective coupling, which is proportional to the density, that is, a term U⁰⁰⁰⁰_iiija_ia^†_j(n_i + n_j − 1) + 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 < 0 for s > 1) of the tunneling amplitude J → J − U⁰⁰⁰⁰_iiij(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 |ψⁿ₀〉 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

$$\begin{equation} H^{\mathrm{EMO}}= - \sum_{\langle i,j \rangle} \left(a_i a_j^\dagger \sum_{n_i,n_j}J_{n_i,n_j} P^i_{n_i}P^j_{n_j}+ \mathrm{h.c.}\right) + \sum_{n,i}\frac{U_n}{2}n(n-1) P^i_{n}, \end{equation} \tag{ 10 }$$

where Pⁱ_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.

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 Δ = E₂ − E₀. Promoting a particle from the lowest to the second excited band, via interaction term, is proportional to n^3/2U²⁰⁰⁰ where U²⁰⁰⁰ is the corresponding interaction integral. This yields for the validity of the perturbative approach the condition n^3/2U²⁰⁰⁰ ≪ Δ. For typical parameters corresponding to Rb scattering length and s of the order of 30, Δ/U²⁰⁰⁰ ≈ 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 |ψ₀〉 be an n-boson ground (Fock) state of H_E. The larger the matrix element $|\langle \psi _0 | \mathcal {H}_{\mathrm {loc}}| \psi _p\rangle |^2$ and the smaller the energy $\langle \psi _p | \mathcal {H}_{\mathrm {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}_{\mathrm {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/λ = 0.014,λ = 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 ≈ 7–9 bands are included. The false saturation occurs because the least-energy method does not evaluate $|\langle \psi _0 | \mathcal {H}_{\mathrm {loc}}| \psi _p\rangle |^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₀/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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

$$\begin{equation} U_{n-1}-U_n\approx W \end{equation} \tag{ 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³[U²⁰⁰⁰]²/Δ. That yields a crude estimate for W ≈ U²⁰⁰⁰/Δ ≈ 1/60 in good agreement with figure 1.

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. 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):

$$\begin{equation} H^{\mathrm{BH}}[\psi]= \langle\psi|H^{\mathrm{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 \end{equation} \tag{ 12 }$$

using the Gutzwiller ansatz: $|\psi \rangle = \bigotimes |\psi _{l}\rangle$ with $|\psi _{l}\rangle = \sum 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:

$$\begin{eqnarray} &&\fl H^{\mathrm{EMO}}[\psi]= -2z \sum\limits_{n_1,n_2}J_{n_1,n_2} \langle \psi_l | a_i | n_1 \rangle \langle n_1 | \psi_l \rangle \langle \psi_l | a_i^\dagger | n_2 \rangle \langle n_2| \psi_l \rangle \nonumber\\ &&+\sum\limits_n \frac{U_n}{2} n(n-1)|\langle \psi_l | n \rangle | ^2 - \mu \sum\limits_n n|\langle \psi_l| n \rangle|^2 . \end{eqnarray} \tag{ 13 }$$

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 $\frac {zJ_{\mathrm {BH}}(s)}{U_{\mathrm {BH}}(s)}$ is strictly decreasing with the lattice height s, it provides a way to plot results calculated for lattice depth s in the $(\frac {zJ_{\mathrm {BH}}}{U_{\mathrm {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_ni,nj/U_n do not need to be directly related to $\frac {zJ_{\mathrm {BH}}}{U_{\mathrm {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_⊥. We show just a single choice for a generic value of s_⊥ = 34.8.

Zoom In Zoom Out Reset image size

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/λ as appropriate for ⁸⁷Rb; the small difference between our results and that of [39] for unit filling on the scale of figure 1 of [66] is negligible.

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_⊥ = 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 μ₊(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.

Zoom In Zoom Out Reset image size

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₀, exposed to a time-varying lattice height s(t) = s₀ + s_m cos ωt, s_m/s₀ ≪ 1. The simulation is performed in the presence of a harmonic trap V_i = κ(i − i₀)². 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 → 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 μ_i = μ − V_i, the excitation energy is determined by the occupation numbers of the source site i and destination site j. It is

$$\begin{eqnarray} &&\fl \Delta E(n_i;n_j)= \textstyle\frac{1}{2} \left[ -U_{n_j} n_j(n_j-1) + U_{n_j+1} n_j(n_j+1)\right.\nonumber \\ &&\left.-U_{n_i} n_i(n_i-1) +U_{n_i-1} (n_i-1)(n_i-2) + \Delta \mu_{ij}\right] \end{eqnarray} \tag{ 14 }$$

with Δμ_ij = μ_i − μ_j = V_j − V_i. Nearest-neighbor 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 equation (11), we may approximate equation (14) by ΔE(n_i;n_j) = (n_j − n_i + 1)U_ni + (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 ΔE(1,1) = U₂, ΔE(1,2) = 3U₃ − U₂ ≈ 2U₂ − 3W, ΔE(2,2) = 3U₃ − 2U₂ ≈ U₂ − 3W, ΔE(2,3) = 6U₄ − 3U₃ − U₂ ≈ 2U₂ − 9W, ΔE(3,3) = 6U₄ − 6U₃ + U₂ ≈ U₂ − 6W, ΔE(n,n − 1) = 0, with Δμ 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.

Zoom In Zoom Out Reset image size

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₀), and the excitation spectrum from Mott plateau n consists of 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 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 κ = 0.001 ((b), moderate case) and N = 700,κ = 0.0001 ((c), very shallow trap). Both cases were studied for the 1D optical lattice with s = 15, s_⊥ = 40, λ = 830 nm,a = 5.1 nm. The tunneling J_ni,nj 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, κ = 0.009, and use true values of the hopping constants J_ni,nj. The Wannier function calculations give U_BH = 0.662E_R, and the renormalization procedure gives U₂ ≈ 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 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ℏ/E_R. The modulation amplitude of the lattice was s_m = 1. The results are presented in figure 7.

Zoom In Zoom Out Reset image size

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₂ 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₂ 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₂ (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₂.

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ℏ/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.

Zoom In Zoom Out Reset image size

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 τ = 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 λ₂, is superimposed along the tubes. This adds a potential V_d(x) = s₂ sin²(k₂x). For s₂ ≪ 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 \mathrm {i} + \phi),$ where ϕ represents the offset between the two optical lattices. We take a generic value ϕ = 0.123 45. If the ratio λ/λ₂ 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₂ on time is set by demanding that s₂(t) ∼ s(t) (all optical lattices are ramped up simultaneously). We will consider three cases: no disorder (s₂(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₂ = 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ℏ/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, U₃ = 1.96 and U₄ = 1.91 kHz. This suggests that the positions of absorption peaks due to the n = 1,2 and 3 Mott plateaus are: U₂ = 2 kHz, −2U₂ + 3U₃ = 1.85 kHz, U₂ − 6U₃ + 6U₄ = 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]).

Zoom In Zoom Out Reset image size

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 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⁻⁹, beyond the accuracy of the calculation).

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.

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.

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

$$\begin{equation} \mathcal{H}_{\mathrm{loc}}= H_E + H_U = \sum\limits_\alpha E_\alpha \hat n_\alpha + \sum\limits_{\alpha\beta\gamma\delta} U^{\alpha\beta\gamma\delta} \hat a^\dagger_\alpha \hat a^\dagger_\beta \hat a_\gamma \hat a_\delta. \end{equation} \tag{ A.1 }$$

The most elementary possible excitation promotes a boson from the α = 0 to 2 band. This defines two limits: n^3/2U²⁰⁰⁰ ≪ E₂ − E₀ and n^3/2U²⁰⁰⁰ ≫ E₂ − E₀. U^αβγδ becomes smaller as the α,β,γ,δ 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_{(0,0,0)}^\dagger )^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 kth-order perturbative expansion, we call those Fock states |ψ_p〉 for which 〈ψ₀|H^k_U|ψ_p〉 ≠ 0. In particular, a full basis can be generated with order $\left \lceil \frac {n}{2}\right \rceil$. Let us denote by $\mathcal {B}_k$ the set of vectors reachable in kth order and unreachable in (k − 1)th order. The full variational basis of Fock states is $\mathcal {B}=\bigcup \limits {_k} \mathcal {B}_k$ with $\mathcal {B}_k$ pairwise disjoint, and $\mathcal {B}_k = \emptyset , \textrm { for } k{>}\left \lceil \frac {n}{2}\right \rceil .$

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

$$\begin{equation} f_1(\psi) = \ln \frac{ |\langle \psi_0 | \hat H_U | \psi \rangle |}{E_{|\psi\rangle} -E_{|\psi_0\rangle}}. \end{equation} \tag{ A.2 }$$

Disregarding vectors from $\bigcup \limits {_{k\geqslant 2}} \mathcal {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 $\mathcal {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,

$$\begin{equation} f_2(\psi) = \ln \sup_{|\psi_1\rangle\in\mathcal{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})}, \end{equation} \tag{ A.3 }$$

is chosen to approximate the relevancy measure for vector |ψ〉. It minimizes the state's |ψ〉 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 = 10 000 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 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.

Zoom In Zoom Out Reset image size

Let us describe in detail how a choice of the basis with the largest f₁ and f₂ 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. $\mathcal {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 $\mathcal {B}_1$ set, so from this point on we assume that k ⩾ 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 $\mathcal {V}=(|\psi _0\rangle ,|\psi _1\rangle ,|\psi _2\rangle ,\ldots ,|\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{eqnarray} &&g_k(\mathcal{V})= \ln \frac{ |\langle \psi_0 | \hat H_U | \psi_1 \rangle \langle \psi_1 | \hat H_U | \psi_2 \rangle |,\ldots, \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}),\ldots, (E_{|\psi_1\rangle} -E_{{|\psi_0\rangle}})},\nonumber \\\\ &&f_k(\psi)= \sup\limits_{|\psi_1\rangle \in \mathcal{B}_1,\ldots,|\psi_{k-1}\rangle \in \mathcal{B}_{k-1} } g_k({\mathcal{V}}).\nonumber \end{eqnarray} \tag{ A.4 }$$

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 |ψ_l〉 preserving the total parity of the state. One-particle update is done according to $|\psi _l\rangle {\to } \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' \,(\textrm {mod }2),$ whereas two-particle update is $|\psi _l\rangle {\to } \hat a_{(i_x,i_y,i_z)}\hat a_{(j_x,j_y,j_z)} \hat a^\dagger _{(i_x,i_y',i_z)}\hat a^\dagger _{(j_x,j_y',j_z)} |\psi _l\rangle , i_y+j_y\equiv i_y'+j_y'\, (\textrm {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 \not \in \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 β 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 |ψ_k〉 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 |ψ_k〉 is memorized as well as the importance value $g_k(\mathcal {V}').$ If |ψ_k〉 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 × 10⁹ 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 > ε. 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 \vec n}$ arises. Logarithm is used to 'smooth' this singularity for numerical purposes. It does not affect the ordering, as ln is increasing injection.