Variational ansatz for the superfluid Mott-insulator transition in optical lattices

We develop a variational wave function for the ground state of a one-dimensional bosonic lattice gas. The variational theory is initally developed for the quantum rotor model and later on extended to the Bose-Hubbard model. This theory is compared with quasi-exact numerical results obtained by Density Matrix Renormalization Group (DMRG) studies and with results from other analytical approximations. Our approach accurately gives local properties for strong and weak interactions, and it also describes the crossover from the superfluid phase to the Mott-insulator phase.


Introduction
During the last years a spectacular development in the storage and manipulation of cold atoms in optical lattices [1,2] has taken place. Greiner et al. [1], to name one important example, succeeded in experimentally driving a quantum phase transition between a superfluid and a Mott-insulating phase in bosonic systems. This experimental progress has revived the interest in the Bose-Hubbard model [Eq. (1)] as a generic Hamiltonian for strongly correlated bosons, by which the quantum phase transition can be described [3]. The Bose-Hubbard Hamiltonian has been used previously in condensed matter physics to study the adsorption of noble gases in nanotubes [4], or Cooper pairs in superconducting films with strong charging effects [5,6]. In this context a lot of work has already been done to characterize the quantum phase transition, the statistics, and the low-energy excitations of the Bose-Hubbbard model [7,8]. However, new interesting questions arise now due to the good tunability of the experiments with optical lattices. In particular, it becomes possible to study time-dependent processes such as driven quantum phase transitions [1]. A theoretical understanding of such phenomena is challenging, since the characteristics of the superfluid phase -where the atoms tend to delocalize throughout the lattice and large fluctuations in the local density exist -, and the Mott-insulating phase -where the number fluctuations decrease, and a gap in the excitation spectrum opens -must be covered at the same time. Both regions are separated by a non-analyticity of the spectrum, which implies that a perturbative study [9,10] works best in strong coupling limit, while a Hartree-Fock-Bogoliubov mean field works best in the superfluid regime. In addition it is possible to develop a mean field theory [11,12] based on a Gutzwiller ansatz [3]: this reproduces the mean field theory in the superfluid limit, as well as the limit of infinite interaction, which raises the hope that the theory also interpolates properly between these limits. In this paper we develop a variational description of the ground state of an ensemble of cold atoms in an one-dimensional optical lattice. Our trial wavefunction treats the connections between neighboring sites as entities which decouple in the limit of infinitely large lattices. We apply this technique first to the quantum rotor model -which describes the lattice for large and commensurate occupation per site [5], and has also been used to describe an array of Josephson junctions [13]-, and next to the Bose-Hubbard Hamiltonian. The accuracy of the variational theory in both models is confirmed by several comparisons. In the quantum rotor case, we use a spin wave approximation in the limit of weak interaction, and a first order perturbation theory in the limit of weak tunneling. For the Bose-Hubbard model we compare against results obtained applying the quasi-exact, numerical DMRG method to one-dimensional lattices with up to 128 sites, and also with calculations based on the Gutzwiller ansatz. Our conclusion is that the variational picture of self-regulated connections between sites provides a rather cheap and simple, but very good description of the local properties of the system in the superfluid and insulator regimes, and a fairly good interpolation across the quantum phase transition. It cannot describe, however, the algebraic decay with distance of the off-diagonal elements of the one-particle density matrix. The simplicity of the method suggests a possible generalization to higher dimensionalities and other physical models.
The outline of the paper is as follows: In Sec. 2 we introduce the quantum rotor model as a possible limit of the Bose-Hubbard Hamiltonian. Next, information about the ground state of the quantum rotor model is obtained variationally as the solution of a Mathieu equation. We can estimate energies, correlation functions and length, and the variance of the density as a function of the only free parameter. A comparison with perturbative estimates demonstrates the accuracy of the method when computing local properties. Since the quantum rotor model is only an approximate description of the optical lattice, in Sec. 3 we develop a similar variational theory for the Bose-Hubbard Hamiltonian. After bringing the Bose-Hubbard Hamiltonian to an appropriate form, we can estimate the local properties of its ground state. The variational solutions are compared in Sec. 3.3 with the results of DMRG studies of the Bose-Hubbard model. We confirm that the variational method describes very well the local properties of both the Mott-insulator and the superfluid regime, and provides a fairly good interpolation across the phase transition. Finally, in Sec. 4 we summarize our results and comment on possible extensions.

2
Quantum phase model

Relation to the Bose-Hubbard model
In this section we show the equivalence of the Bose-Hubbard model and the quantum rotor model for large and integer occupationn [5]. In our notation, M is the number of lattice sites and N =nM the number of atoms. Both the Bose-Hubbard model and the quantum rotor model show a phase transition due to the interplay between the kinetic term proportional to J and the interaction term proportional to U . For convenience we have subtracted the ground state energy in the perfect insulator limit U/J → ∞.
If we expand a configuration of the lattice using Fock states |ψ = n c n | n = n c n |n 1 ⊗ · · · ⊗ |n M , and the number of particles per lattice site is large, n k > 1, we may approximate the hopping terms as follows a † l a j |ψ = n(n + 1) Here, A ± |n = |n ± 1 are ladder operators and P projects on states with non-negative occupation numbers, n k ≥ 0. To lowest order the error |∆ lj is and its norm is bound by where σ l = (n l −n) 2 is the variance in the number of particles per lattice site. For the approximation (3) to be valid, the uncertainty in the number of atoms must be small compared to the mean value,n ≫ σ, and the interaction energy must exceed the neglected terms, Un(n − 1) ≫ Jσ. Following the previous procedure the Bose-Hubbard model becomes Here ρ = n(n + 1) is approximately the density, A z j = a † j a j −n is essentially the number operator and we have used that j A z j |ψ = 0 when we work with states that have a fixed, commensurate number of particles. In the following we will define the energy per lattice site as Since the physically interesting states will be concentrated around large occupations, n k =n, the usual step now is to drop the projector, P, and move to the basis of phase states, defined by In doing so, we obtain the identification A ± j → e ±iφ k and A z j → −i∂/∂φ j , which produces the usual representation of the quantum rotor model with the associated state writen as A similar derivation is possible using path integrals [15].
Instead of working directly with the population of each well, n k , we can use other quantum numbers, w k , defined by the relation n k = w k − w k−1 +n, and which behave like a set of chemical potentials acting on the barriers that connect neigboring sites.

Variational ansatz
In this section we estimate the properties of the ground state of H QR variationally. Due to the previous splitting (10), any wavefunction Ψ( φ) can only depend on the phase difference between neighboring wells, ξ j = φ j+1 − φ j . These new quantum variables describe the connections between neigboring sites. In the limit of large lattices it seems reasonable to assume that these connections become independent from each other adopting the product state This representation becomes exact in the Mott-insulating limit, U/J → ∞, where h mott (ξ) = 1, and in the superfluid limit, , as can be verified by direct substitution in Eq. (10). Even though the phase representation is the best one to find a trial wavefunction, it is not the optimal one for performing computations. It is instead more convenient to work with the variables which are conjugate to the phase differences ξ j = φ j+1 − φ j . These are the new quantum numbers, w k , given by In terms of these numbers, the ansatz (11) reveals itself as a simple product wavefunction with coefficients given by the Fourier transform As sketched in Fig. 1, the w k play the roles of chemical potentials which are established between different wells: the difference between the potentials on the extremes of a site gives the fluctuations over the mean and commesurate occupationn. In this picture where Σ ± |w = |w ± 1 are new infinite-dimensional ladder operators and Σ z |w = w|w . By minimizing the energy associated with H QR over all states within a given ansatz we can both obtain an upper bound to the energy of the ground state and approximate its wave function. A simple computation with our product ansatz leads to the result [compare (7)] where the expected values are computed over a single connection, Σ z = w w|h w | 2 , and the wavefunctions are assumed to be normalized, w |h w | 2 = 1. Since the stationary states have a well defined parity,h (−w) = (−1) Ph w , the optimal variational state must satisfy the linear equation which is nothing but the Fourier transform of a Mathieu equation The estimated ground state energy per site is given by the lowest eigenvalue of either equation.
Using the product ansatz and the same approximations required to derive H QR , we can also compute other properties of the ground state. For instance, the variance of the number of atoms per lattice site and the correlation functions which decay exponentially with the distance. This implies that the ansatz (13) only describes properly the decay of the correlations in the Mott-insulating regime, since the correlations in the superfluid regime follow a power law decay. However, as we will see below, local magnitudes (σ, γ 1 , ε . . .) are properly estimated even if long-range ones are not.
We have solved Eq. (17) numerically in a truncated space. The results are summarized in Fig. 2, where we also plot reference estimates arising from two other analytical methods. In the limit U ≫ J we compare with a first order perturbative calculation around the solution |ψ = 1, which is possible thanks to the energy gap of order O(U ) in the excitation spectrum. In the limit J ≫ U we rather use a spin wave or harmonic approximation in which the cosine-term of the Hamiltonian H QR is expanded up to second order in the phase difference between neighboring sites [see Sec. 2.3]. This approximation is valid in the superfluid regime, where the phase does not vary much between neighboring wells. From the graphical comparison we see that, as expected, the variational wavefunction provides a fairly accurate description of the ground state of the quantum rotor model in both the superfluid and insulating limits. As a side note, we must remark that this ground state has a divergent fluctuation of the number of particles per site as J → 0. This disagrees from the expected behavior of the ground state of the original Bose-Hubard Hamiltonian, and it reminds us that H QR can only model the atomic ensemble when the variance, σ, is small compared to the mean occupation number,n.

Harmonic approximations to the quantum phase model
In the limit J ≫ U it is possible to estimate the ground state of the rotor model (6) analytically. Since we are deep in the superfluid regime, the wavefunction will be concentrated around the line φ 1 = φ 2 = . . . = φ M , and we can approximate If we remove the periodic boundary conditions on φ i and change variables, the preceding Hamiltonian may be diagonalized, H = k ω k A † k A k + 1 2 , with frequencies given by The variance of the number of particles is related to the expectation value of the momentum using the same procedure as above 3 The Bose-Hubbard model In this section we apply to the Bose-Hubbard Hamiltonian the techniques that were developed in Sec. 2. We will do it in three steps: First we will develop a phase representation which is valid for all occupation numbers. Next we will prove that this representation is equivalent to a similarity transformation of the Hamiltonian which brings it to a form similar to (6), at the price of losing Hermiticity. Finally we will show how to implement the ansatz of independent connections (13) to produce estimates for the usual set of observables (ε, γ 1 , σ), which are to be validated with DMRG calculations.

Coherent states
The phase coherent states |φ are defined by Unlike the phase states defined in Sec. 2, they are not orthogonal to each other, φ|θ = exp e i(θ−φ) , but they form a complete basis, so that an expansion like (10) is still possible. A nice property of the coherent states is that we can rewrite the operators a, a † , a † a, etc, in terms of the phases very easily. For instance, a|φ = e iφ |φ , a † |φ = −ie −iφ ∂ ∂φ |φ , and a † a|φ = −i ∂ ∂φ |φ . Using this representation, we obtain an effective Hamiltonian for the wavefunction Ψ( φ), i.e.
Here H t coh stands for the tranpose of H coh . This operator was already used in Ref. [28] to study the Bose-Hubbard model with only two sites. On the one hand, it is a non-Hermitian operator 1 and we cannot do a simple variational study. On the other hand the Hamiltonian still depends on the phase differences, and it is reasonable to look for approximate eigenstates which have the form (11). This will be done in the following section. However, since working with phase variables is inconvenient, we will develop a representation similar to that of Eq. (15) in the following.

Variational procedure for non-Hermitian operators
In this section we will find the best variational function which has the product form of Eq. (11). However, as it happened in Sec. 2, instead of working with phase variables it is more convenient develop a representation of the Bose-Hubbard Hamiltonian in terms of connections. This is once more a two-steps process. First we use a similarity transformation suggested by the definition of the coherent states Since The Hamiltonians (26) and (28) are equivalent: while one is defined in terms of phase variables, the other one is defined with occupations numbers, and both are related by a Fourier transform. The second and final step is to rewrite everything in terms of connections, using the quantum numbers from Eq. (12), and the relations Σ x = Σ + +Σ − , Σ y = i(Σ − − Σ + ). The result is a decomposition of the Hamiltonian into terms which are local, H 1 , and terms which involve pairs of connections, H 2 .
For the quantum rotor model we proved that the optimal product wavefunction was an eigenstate of a Hamiltonian which did not couple connections, like H 1 . The difference now is that, since the operator H coh is not Hermitian, we cannot establish a variational principle and that proof is no longer valid. Nevertheless, we will again propose a variational ansatz which is an eigenstate of the local operator H 1 |h ⊗M = M ε est |h ⊗N . Using the following equality and the product ansatz |χ = |h ⊗M , we arrive to an upper bound for the lowest eigenvalue of the Bose-Hubbard Hamiltonian, expressed in terms of the non-Hermitian one The way to use this variational principle is as follows. First, for a given J and U we compute the lowest eigenstate of H 1 and this way obtainh. After the equivalence (14), finding the ground state of the local Hamiltonian H 1 becomes equivalent to solving a modified Mathieu equation which describes exactly the static properties of a pair of sites with open boundary conditions [28]. Once we have ε est we must still compute the correction ∆ε est using a rather straightforward expansion which is shown in Sec. 3.4. Surprisingly, ∆ε est happens to be negative, so that it is actually an improvement over the simple estimate given by ε est [See Fig. 5]. From the optimal variational state, |ψ = O −1 |h ⊗M , and the estimate for the energy, ε var = ε est + ∆ε est , we may compute other observables. For the density fluctuations and nearest neighbor correlations we use the virial theorem whereas for other properties one has to evaluate numerically the matrix products shown in Sec. 3.4. This allows us to prove that for the product states |h ⊗M correlations decay exponentially, opposite to what is expected in the superfluid phase, whose correlations should decay algebraically. Nevertheless, as we will see next, this family of states does estimate accurately the local properties of the optical lattice.

Comparison to DMRG results
We will now compare the results for the ground state energy, the correlation functions, and the variance of the particle number provided by the two variational ansatz, (18) [10]. The mean occupation numbers are denoted with circles (n = 1), diamonds (n = 2) and boxes (n = 3). and (32), and the Gutzwiller ansatz [14] with those obtained by DMRG studies of the Bose-Hubbard model. The DMRG, developed 1992 by White [16,17] in the area of condensed matter theory, is a very powerful numerical tool to investigate static and dynamic properties of strongly correlated quasi-one-dimensional spin, fermionic or bosonic quantum systems. The DMRG is an essentially quasi-exact numerical method. The fundamental ideas stem from real space renormalization methods: the system size is grown iteratively while the (exponentially diverging) size of the Hilbert space is kept constant by decimation. Hereby one tries to retain only that subset of states that is essential to describe the physical quantity under consideration. In DMRG these are expectation values with respect to low-lying states ("target states"), and in particular with respect to the ground state wave function.
DMRG builds up the system linearly: at each growth step, suitable density matrices for the target states are derived that yield information on the relevance of Hilbert space states. Building on this information, the states and operators are projected onto Hilbert subspaces of fixed dimension M containing the most relevant states. M is chosen to be small enough to be handled numerically, but large enough to obtain the desired accuracy; numerical results can be extrapolated in M to the exact limit of infinite M in the thermodynamic limit. However, results presented here have converged for the largest M considered and no further extrapolation was necessary.
Details on the DMRG method may for example be found in [18]. In the case of the Bose-Hubbard model the DMRG has been used to study properties of the system [19,20,21]. The results of DMRG agree very well with exact diagonalization results for small systems, with quantum Monte-Carlo simulations e.g. [22,23,24,25], and with 13th order perturbation theory [26].
We have used the DMRG to study the properties of the ground state of the Bose-Hubbard model on one-dimensional lattices with 128 sites, and commensurate fillings n = 1, 2, and 3. In Figs. 3(a-c) we show the results for the mean energy per site ε, the nearest neighbor correlations, c 1 = a † j+1 a j , and the variance σ of the density, calculated both with the DMRG and with the variational estimates developed above. As expected, there are no indications of the phase transitions in these quantities, neither in the variational results nor in the numerical solutions. Rather, an inflexion of the nearest neighbor correlation points out the location of the superfluid-insulator transition which lies roughly between 3n and 4n (see [19] and ref. therein). The agreement of the two methods is fairly good above the phase transition and below it.
A more detailed comparison is provided in Fig. 4 for the casen = 1. In this figure we plot together results from the DMRG, the variational ansatz derived above, the quantum rotor model and the well-known Gutzwiller ansatz. The Gutzwiller ansatz [27,3] is a variational ansatz which reduces the wave function to a product of single-site wave functions, m |m j and f (j) m are constants. Such a wavefunction cannot be used in the one-dimensional Mott insulator regime, because a perturbative study of the Gutzwiller ansatz shows that the corrections of order O(J/U ) are lost and all correlations become zero. However this ansatz gives good results in the superfluid regime, where the long-range order is well described by |Ψ G , and we can use these results and those of the DMRG to assert the accuracy of our variational estimates. As Fig. 4 shows, in the Mott insulator regime, the DMRG results agrees perfectly with our variational theory for the Bose-Hubbard model and for the quantum rotor model. Close to the phase transition is the point at which the quantum rotor model no longer describes well the atoms in the optical lattice due to the growth of density fluctuations. At this point we also observe a small disagreement between the DMRG and the coherent states, which is due to the growth of long range correlations and vanishes as we get deeper into the superfluid regime.

Numerical evaluation of the upper bound
In this section we will show how to compute the corrections to the local energy, ∆ε est , from Eq. (31). We basically need a method to compute expectation values of the operator O −2 around product states which have the form in which at most two contiguous vectors are affected by two single-connection operators. For instance, this is the case of H 2 |ψ , where A and B are Σ + , Σ − or Σ z , and also of |ψ , where A and B are just the identity. After some manipulations it is possible to write, for the optimal variational state |ψ = O −1 |h ⊗M , with the real matrices and vectors  [10].
We have used this technique to compute numerically the correction ∆ε est using different lattice sizes and found small or no differences for more than 30 sites. Intuitively, this is because in the limit of large powers the matrices (HO) α become projectors on the eigenvector with the largest eigenvalue.
Using the same type of expansion we may compute other correlations For large values of ∆ and large lattices, the numerator will decay exponentially as γ ∆ 1 , where γ 1 is the largest eigenvalue of the matrix H − O.

Conclusions
In this work we have studied analytically and numerically the properties of the ground state of an ensemble of bosonic atoms in an 1D optical lattice. For the study of the The energy of the product ansatz contains a contribution from each connection, εest, plus the interaction between neighbouring connections, ∆εest. In Fig. (a) we show that ∆εest (dash) is actually negative, and improves the estimate εest moving it towards the exact value, ε DM RG (circles). Everything has been computed forn = 1. In Fig. (b) we show that the correction ∆εest does not change much for large lattices.
atomic ensemble we have used both the quantum rotor model and the Bose-Hubbard model. Exploiting the fact that in these models there exists only nearest neighbor hopping and local interactions, we have developed a variational wavefunction that may be used to easily estimate local properties, such as the energy per well, the nearest neigbor correlations and the fluctuations of the density. In the case of the quantum rotor model we have verified our results with perturbative calculations around the strongly interacting regime, and with a spin wave approximation around the superfluid regime. In the case of the Bose-Hubbard model we have compared the variational estimates with numerical results obtained using the DMRG technique for a maximum density of three atoms per well. We have concluded that this procedure leads to fairly good estimates of local ground state properties of both Hamiltonians, in both the superfluid and the insulator regime, the largest disagreement being localized around the phase transition.
In future work we forsee the application of these methods to time dependent problems.