An exact reformulation of the Bose-Hubbard model in terms of a stochastic Gutzwiller ansatz

We extend our exact reformulation of the bosonic many-body problem in terms of a stochastic Hartree ansatz to a stochastic Gutzwiller ansatz for the Bose Hubbard model. This makes the corresponding Monte Carlo method more efficient for strongly correlated bosonic phases like the Mott insulator phase or the Tonks phase. We present a first numerical application of this stochastic method to a system of impenetrable bosons on a 1D lattice showing the transition from the discrete Tonks gas to the Mott phase as the chemical potential is increased.


I. INTRODUCTION
The Bose-Hubbard model of lattice bosons with on-site interactions [1] is a widely used theoretical model for the description of a number of different physical systems, ranging from 4 He adsorbed in porous media [2] to granular superconductors [3,4] and Josephson junction arrays [5,6,7,8,9,10]. In the last years, the Bose-Hubbard model has attracted a renewed interest since its experimental realisation with ultracold atoms loaded in the periodic potential of an optical lattice [11,12]. Following the suggestion of [13], this system has led to the observation of the superfluid-Mott insulator quantum phase transition [12] and of quantum revivals of the matter field [14]. In the present paper we propose a new tool for the study of the Bose-Hubbard model in the form of a Quantum Monte Carlo scheme based on the stochastic evolution of a Gutzwiller ansatz. This new scheme generalizes our previously proposed schemes based on a stochastic Hartree ansatz [15,16,17]. Intuitively, it is clear that the Hartree ansatz is not the optimal one for the study of a Mott insulator state: the expansion of a Mott-insulator state on Hartree states involves an average over all the possible relative phases between adjacent sites, so that a large number of realizations is required for an accurate sampling of the state. Furthermore, if one rises the on-site interaction strength to infinity so to realize the so-called impenetrable boson model, the rate of increase of the statistical error with time, being proportional to the strength of the atom-atom interactions, will diverge. On the other hand, approximate analytical mean-field calculations based on a Gutzwiller-type wavefunction have given predictions which are in reasonable agreement with the exact numerical calculations [18]. This suggests that a suitable stochastic generalization of the Gutzwiller mean-field equations could lead to an efficient quantum Monte Carlo scheme. The main difference with respect to the reformulation in terms of a stochastic Hartree ansatz consists in the exchanged roles of the kinetic and the interaction energies: in the Gutzwiller case, interactions can be fully taken into account by the deterministic evolution, while the noise is necessary to keep track of the correlations between different sites which are induced by the kinetic energy term in the Hamiltonian. While the stochastic Hartree reformulation was well suited for the study of weakly interacting Bose gases [19], the present Gutzwiller scheme is expected to be more suited to the study of systems of strongly interacting bosons, such as the Mott phase or the gas of impenetrable bosons. In sec.II we review the Bose-Hubbard Hamiltonian and we introduce the Gutzwiller ansatz. In sec.III we show how it is possible to find a suitable stochastic evolution of the Gutzwiller wave function such that the stochastic evolution of the ansatz is equal, in the average, to the exact evolution given by the Bose-Hubbard Hamiltonian (1). In sec.IV we discuss the evolution of the statistical error of this stochastic reformulation of the many-body problem and, in particular, we show that the error remains finite at any time of the evolution. This guarantees that the stochastic Gutzwiller reformulation can be actually used for numerical Monte Carlo calculations. In sec.V, we extend the reformulation to the case of an imaginary-time evolution so as to obtain the thermal equilibrium state of the Bose-Hubbard system at a given temperature and chemical potential. A first application to the case of impenetrable bosons on a 1D lattice is presented, which shows the transition from a discrete Tonks gas to the Mott phase as the chemical potential is increased: the results are in good agreement with analytical calculations performed by means of fermionization techniques. For the sake of completeness, a brief discussion of these techniques is given in the appendix. In sec.VI we generalize the results to the case of several spin components. Conclusions are drawn in sec.VII.

II. THE MODEL HAMILTONIAN AND THE GUTZWILLER ANSATZ
In its simplest form, the Bose-Hubbard Hamiltonian can be written as: in which the spatial coordinate r runs on a D-dimensional cubic lattice of N = D i=1 N i points with periodic boundary conditions. With l i and L i = N i l i , we shall denote respectively the lattice spacing and the total lattice size along the i-th dimension. The on-site interaction is taken into account by the first term H U , which is proportional to the charging energy U . The hopping from one site to a neighboring one is described by the other term H J . J quantifies the amplitude of the hopping processes. The function δ (1) r,r ′ is equal to one if r and r ′ are neighboring sites, otherwise it vanishes. The annihilation operatorsΨ(r) satisfy the bosonic commutation rules [Ψ(r),Ψ † (r ′ )] = δ r,r ′ . The generic Gutzwiller ansatz can be written in terms of the complex wavefunctions φ(r, n) as: This kind of states has been extensively used within mean-field theory to obtain approximate analytical results for strongly interacting bosons on a lattice: since the wave function has a factorized form, no correlations between different sites are included and the effect of H J can not be treated exactly. This is complementary to what happens in the case of a Hartree ansatz where the mean field theory can deal exactly with the hopping Hamiltonian, but not with the interaction one. The generality of the Gutzwiller ansatz is however large enough to include Glauber coherent states as well as Fock states. A coherent state of macroscopic wavefunction γ(r): where γ 2 = r |γ(r)| 2 , corresponds to a Gutzwiller state with: On the other hand, a Fock state in which each site has a well defined occupation n 0 (r) corresponds to a Gutzwiller wavefunction of the form: The Gutzwiller ansatz is therefore much more flexible than the coherent state ansatz often used in quantum optics. Its normalization is given by: III. REAL TIME EVOLUTION reabsorb the phases of the complex coefficients of this superposition by a redefinition of the Gutzwiller wavefunctions, so that the initial state vector can be written as the average of the Gutzwiller ansatz over some probability distribution for φ: where we represent this average by angular brackets . . . . In the present section, we shall show that a suitable stochastic evolution can be found for the Gutzwiller wavefunctions φ(r, n), such that a relation like (7) holds at any time t > 0, in the sense that the exact state vector at time t is given by the average over the ensemble of Gutzwiller states resulting from all possible stochastic evolutions of the initial ensemble during the time t: In the hopping Hamiltonian H J , we first explicitely isolate the mean-field contribution by rewriting: in terms of the C-number mean-fieldψ(r). defined as: where the ket |φ(r) describes the state of the site r: This corresponds to the general strategy, already implemented for the Hartree-Fock ansatz for bosons [17] and for fermions [20], to have a deterministic evolution identical to the mean-field theory so that the stochastic part includes only deviations from mean-field. Mathematically, as we shall see, this makes the noise term strictly orthogonal to the state |G : φ , which is expected to minimize the growth rate of the statistical error. Note that only the first term of the square bracket (9) can induce correlations between different sites, and it is this term which will be taken into account with the noise term of the evolution. The action of the total Hamiltonian (1) on a Gutzwiller state of the form (2) can then be written as: where the matrices χ(r; n, n ′ ) have vanishing entries for |n − n ′ | > 1, and so to include all the on-site terms coming from both the interaction and the mean-field part of the hopping Hamiltonian (9). The factor 1/N appearing in the last term in (13) comes from the normalization (6) of the Gutzwiller state. We now consider a generic stochastic evolution for φ(r, n): where dφ d (r, n) is a deterministic drift term proportional to dt and dφ s (r, n) is a stochastic noise in the Ito sense. In particular, the mean value of dφ s is zero and its variance proportional to dt. The average evolution of the Gutzwiller ansatz during dt can be obtained by means of an expansion in powers of dt; retaining the terms up to order dt, one is left with: The stochastic evolution (17) can be made to coincide with the exact one given by (12) simply by imposing that the deterministic evolution dφ d is equal to: and that the correlation function of the noise term satisfies: The symmetrized form of (19) with respect to (r, n) and (r ′ , n ′ ) is imposed by the general property of a noise correlation matrix being symmetric. Apart from an additional term giving a mere phase rotation, the deterministic evolution (18) coincides with the usual evolution considered in mean-field theory. It is easy to verify by direct substitution that a noise with the required correlation function (19) can be obtained as: with dB k being independent and zero-mean complex random variables such that: For the lattice here considered, the quasi-momentum k is a D-component vector whose ith component is equal to k i = 2πn i /L i , n i being an integer number 0 ≤ n i < N i . ω(k) is the particle dispersion in the lattice in the absence of interactions: In any practical simulation, a cut-off N max has to be set in the sums over the site occupation numbers, which physically corresponds to projecting the dynamics on the subspace in which all sites have an occupation number smaller or equal to N max . The error induced by imposing this cut-off is negligible provided a sufficiently large value of N max is adopted.

IV. EVOLUTION OF THE STATISTICAL ERROR
As we have discussed in detail in [17], formally imposing the correct value of the averages is not in general sufficient for a practical application: it is essential to prove the finiteness of the statistical error, which guarantees the convergence of the Monte Carlo simulation to the exact result in the limit of an infinite number of independent realizations. As usual, the variance of the statistical error on the state vector of the system is defined in terms of the stochastically evolving Gutzwiller ansatz |G : φ(t) and the exact state vector |ψ(t) as: Using the definition (8) of the stochastic reformulation, one can write the evolution of |G : φ(t) during an interval dt as: where the Ito noise term d|G : φ(t) s depends linearly on the stochastic term dφ s (r, n) of the evolution of the Gutzwiller wave function. Our specific choice of noise terms (20) not only leads to the orthogonality of this Ito term to the Gutzwiller ansatz: but, even more strongly, to a purely orthogonal stochastic variation of the state of each site r: Using the fact that the square modulus of the exact state vector |ψ(t) 2 is constant and equal to 1 during the evolution, the variance of the statistical error can be rewritten as: with: Inserting the explicit form of the stochastic noise term (20) and taking advantage of (27), one obtains an upper bound to the variation of d∆ in the interval dt: and therefore on ∆(t) itself: This shows that the stochastic trajectories for the Gutzwiller ansatz cannot escape to infinity in a finite time. By ensemble averaging, one immediately finds an upper bound to the variance E(t) of the statistical error (28): This inequality guarantees that the statistical error of the stochastic reformulation using Gutzwiller states remains finite at all evolution times. This result is to be compared with the one we have obtained in [15] for the Fock simple scheme: where N is the number of atoms. It is apparent how for N ≈ N and N max ≈ 1, the present scheme is more suited to simulate a system in the Mott-insulator phase U ≫ J, for which the interaction energy dominates over the kinetic energy.

V. IMAGINARY TIME EVOLUTION
The previous sections have been devoted to the discussion of a stochastic reformulation of the real-time evolution of the Bose-Hubbard model. On the other hand, the density matrix ρ = exp(−H/k B T ) describing the thermodynamical equilibrium state of a system at a given temperature T can be obtained by means of the imaginary-time evolution: during an imaginary-time interval τ = 0 → β = 1/k B T starting from the identity ρ(0) = 1 at τ = 0. The expectation value of any operatorÔ can be then obtained from ρ(β) as: As we have done in [16] for the Fock simple scheme, the stochastic reformulation in terms of a Gutzwiller ansatz originally developed for the real-time evolution can be easily extended to the imaginary-time evolution. It follows from the completeness of coherent states that the identity matrix can be expanded as: in terms of coherent states of the form (3). The integration element Dγ is defined as usual as: According to (36), the imaginary-time evolution therefore starts from initially identical Gutzwiller wavefunctions φ 1,2 (r, n; 0) of the coherent state form (4) with a randomly chosen γ(r) and proceeds as described by a pair of stochastic differential equations for φ α=1,2 : Apart from the obvious substitution dt i → − dτ 2 , the imaginary-time stochastic equations have exactly the same form as (18) and (20). As the Gutzwiller ansatz does not have in general a well defined number of particles, it is more convenient to use the Grand-Canonical ensemble with a chemical potential µ, so that H = H − µN ,N being the operator giving the total number of particles. The deterministic part of the evolution then reads: and the stochastic part: with dB k being independent and zero-mean complex random variables such that: In the practical simulation presented below, the dB k have a Gaussian distribution. For the left and right Gutzwiller ansatz α = 1, 2,ψ α (r) is defined as the mean-field corresponding to the wavefunction φ α (r, n): The finiteness of the statistical error at any time of the imaginary-time evolution can be proven by means of the same arguments invoked above for the case of a real-time evolution. The evolution of the state vector during dτ is given by: As the spectrum of the Bose-Hubbard Hamiltonian H is bounded from below in the presence of the cut-off at N max , an upper bound can be set to the growth rate of |G : φ(τ ) 2 , which guarantees that the statistical error remains finite at any time. In both panels, the temperature is related to the (here negative) hopping constant by |J|/kB T = 9/4. In (b) the difference between the chemical potential and the single particle ground state energy J is (9/8)kB T so that the fermionic analog of the Bose system is on the border of being degenerate. The position x is in units of the lattice spacing l.
In the figure, we present a first, simple application of the method, the observation of the finite temperature cross-over from the discrete Tonks gas to the Mott phase in a system of impenetrable spinless bosons. This system for which U = +∞ was totally out of reach for our previous algorithm with Hartree states [16] but can be easily treated by the stochastic Gutzwiller algorithm by simply setting N max = 1. We introduce the name 'discrete Tonks gas' to refer to a phase of impenetrable bosons on a lattice with a filling factor much smaller than unity; as the mean distance between bosons is then much larger than the lattice spacing, we expect this phase to have large scale properties very close to the ones of the continuous model, hence the name. The Mott phase corresponds to the opposite limit of a unit filling factor.

VI. SEVERAL SPIN STATES
All the discussion of the previous sections can be easily extended to the case in which the particles can be in s different spin states. The Bose-Hubbard Hamiltonian (1) is immediately generalized as: Ψ(r, σ) being the destruction operators for a particle at the site r in the spin state σ. The matrix U σ,σ ′ quantifies the on-site interactions between atoms respectively in the σ and σ ′ states. The hopping is assumed not to affect the spin degrees of freedom and is quantified by the parameters J σ . The straightforward generalization of the Gutzwiller ansatz (2) to the case of several spin states has the form: where the Gutzwiller wavefunction φ(r, n) now depends on the spatial coordinate r and on the s-components vector n = (n 1 , . . . , n s ) formed by the occupation numbers of each spin state.
As in the spinless case, both the on-site interactions and the mean-field contribution to the hopping term are taken into account by the deterministic evolution of φ(r, n). The fluctuations around mean-field are described by the stochastic noise term.
For the case of a real-time evolution, the deterministic term reads: All the sums over spin states run over σ = 1 . . . s and e σ is defined as the s-component unit vector whose σ-th component is 1 and all the others are vanishing. The mean-fieldψ(r, σ) is defined as: The noise term can be written as: in terms of the independent and zero-mean complex random variables dB k,σ such that: ω(k, σ) describes the particle dispersion in the lattice in the absence of interactions: As in the previous section, these equations can be extended to the imaginary-time evolution simply by replacing dt i → − dτ 2 . The finiteness of the statistical error can be proven in the same way as previously done for the spinless case.

VII. CONCLUSIONS
In this paper, we have introduced a reformulation of the bosonic many-body problem on a lattice in terms of an stochastic Gutzwiller ansatz and we have proved that the statistical dispersion remains finite at all times in the presence of a cut-off on the number of particles per lattice site. This reformulation allows to perform efficient Monte Carlo simulations for strongly correlated systems like the Mott phase or the Tonks gas. This new method can also be used to study hard sphere bosons in arbitrary spatial dimensions, a situation which was totally out of reach for our previous stochastic Hartree ansatz.
The correlation functions of a gas of impenetrable bosons in one dimension can be computed by mapping the bosonic system onto a fermionic system. The bosonic operatorsΨ(x) on the lattice are written in terms of fermionic operatorŝ Ψ F (x) by means of the Jordan-Wigner transformation [21,22]: In this way, the Hamiltonian (1) for the impenetrable (U = +∞) Bose gas reduces to a free-fermion Hamiltonian of the form: Given the specific ordering given to sites in (A.1), it is in fact necessary to isolate in (A.2) the term corresponding to the tunneling between the last and the first site of the lattice, whose amplitude changes sign depending on the parity of the total number of particles present in the system. This peculiarity of the periodic boundary conditions is well-known in the frame of exactly solvable one-dimensional models, see e.g. [23]. For any value of the total number of particles, the Hamiltonian (A.2) is quadratic in the Fermi field operators and describes a homogeneous system of non-interacting fermions. The single particle states are plane waves of wavevector k and their energy is given by: but the quantization condition on k depends on the parity of N . For odd values of N , the Hamiltonian (A.2) reduces to the usual one for a one-dimensional lattice with tunneling constant J, so that periodic boundary conditions on the single-particle wavefunctions apply and the allowed values for k are given by: in terms of the integer 0 ≤ n < N . On the other hand, for even values of N , anti-periodic boundary conditions are to be applied, so that k is now given by: The calculation of Ô σ = Tr[Ôσ]/Tr[σ] has to be performed using the Jordan-Wigner mapping (A.1) to write the observable in terms of the fermionic operators and then applying Wick's theorem, whose validity can be proven to extend to the present case of a non-hermitian quadratic Hamiltonian. In particular, the first order correlation function g (1) (x) in one of the density operators σ can be shown [22,24] to be given by: Note that we have added a global factor 1/2 which was missing in [22].