Thermodynamics of the Bose-Hubbard model in a Bogoliubov+U theory

We derive the Bogoliubov+U formalism to study the thermodynamical properties of the Bose-Hubbard model. The framework can be viewed as the zero-frequency limit of bosonic dynamical mean-field theory (B-DMFT), but equally well as an extension of the mean-field decoupling approximation in which pair creation and annihilation of depleted particles is taken into account. The self-energy on the impurity site is treated variationally, minimizing the grand potential. The theory containing just three parameters that are determined self-consistently reproduces the T=0 phase diagrams of the three-dimensional and two-dimensional Bose-Hubbard model with an accuracy of 1% or better. The superfluid to normal transition at finite temperature is also reproduced well and only slightly less accurately than in B-DMFT.


I. INTRODUCTION
The properties of cold atomic gases trapped in an optical lattice can be tuned and controlled very precisely, providing a powerful tool for the simulation of the iconical low-energy effective Hamiltonians of condensed-matter models 1 . Dramatic experimental progress in this field such as the observation of the Mott insulator to superfluid phase transition in the Bose-Hubbard model 1 or the recent realization of the Hofstadter model 2 , have galvanized the condensed matter community.
The accuracy of cold atom experiments has laid bare the need for advanced numerical methods in order to tackle these correlated many-body systems quantitatively. Very successful have been path integral quantum Monte Carlo (QMC) simulations with wormtype updates 3 leading to identical results for the Bose-Hubbard model as observed in experiment 4,5 . Despite all its impressive successes for bosonic systems the wormalgorithm suffers from a prohibitive sign problem when cold atoms are subject to (artifical) gauge fields. In such cases no algorithm is known that works and one is forced to resort to approximations. This has been a main driving force for the development of bosonic dynamical meanfield theory 6,7 (B-DMFT).
In B-DMFT, as in standard mean-field theory, the many-body system is projected onto a single impurity, keeping only the local propagators. This provides us with a model in which the self-energy and the local propagators can be determined self-consistently by solving an effective impurity action and a self-consistency condition iteratively. At the moment it has only been formulated for single-site impurities, but the ultimate goal is to formulate it for small clusters that can also deal with larger unit cells. It is known that B-DMFT provides excellent agreement (of the order of 1% in 3d) with experimental and QMC data 6 for the standard Bose-Hubbard model and improves remarkably on static mean-field theory. B-DMFT is hence a promising candidate to deal with more complicated dispersions: the hope is that it deals with the local interactions as accurately as in the standard Bose-Hubbard model whereas a small cluster could treat the full dispersion. Furthermore, one would expect that a cluster method goes beyond real-space DMFT for multicomponent systems and systems with long-range interactions and/or dispersions 8,9 . However, B-DMFT is numerically rather complex due to the imaginary-time dependency of the hybridization term. At finite temperature the impurity problem has to be solved by continuous time Monte Carlo methods 6,7 , where, due to the difference in sign between the normal and the anomalous Green function, a sign problem arises in the symmetry broken phase.
In this work, we filter out the ingredients of B-DMFT that are indispensable for its accuracy and arrive at a simpler formalism. This is the Bogoliubov+U theory (B+U), which makes use of a simplified effective impurity Hamiltonian, similar to the action of extended mean field theory, which was recently developed in the highenergy community 10,11 but differs conceptually from our formalism. B+U has a negligible computational cost and is not prone to numerical instabilities. The premise of our theory is that the Bose-Hubbard model can be fully characterized at zero temperature by the three parameters φ, Σ 11 and Σ 12 if the self-energy is treated as a variational parameter, providing a far better approximation to finite-temperature properties than simple mean-field theory. B+U can be seen as a simplified B-DMFT where only a single Matsubara frequency is kept (and is hence conserving). It is different from the variational cluster approximation by also considering non-zero values of pair creation and annihilation of depleted particles 12 . It is also the simplest accurate extension of the weakly interacting Bose gas theory 13 to lattice systems with a superfluid to Mott insulator transition. It further provides a very natural framework compared to the collective quantum field theory developed by Kleinert et al. 14 and the two collective fields proposed by Cooper et al. [15][16][17] , and behaves quantitatively much better. It has a similar functional degree of freedom as the projector technique introduced by Trefzger and Sengupta 18 for finite lattices.
The paper is organized as follows. In section II the B+U formalism is derived for the Bose-Hubbard model in equilibrium, while in section III the details of the vari-arXiv:1501.07849v1 [cond-mat.quant-gas] 30 Jan 2015 ational calculation of the optimal self-energy are shown. We furthermore summarize the full self-consistent scheme of B+U and show how thermodynamic quantities can be calculated from it in section IV, while some simple limits of B+U are explained in section V. In section VI we present the results at zero and finite temperature comparing them with QMC and B-DMFT. Finally, in section VII we conclude and present a short outlook about future applications of the B+U-formalism.

II. SOLVER AND SELF-CONSISTENCY CONDITION
In this section we derive the B+U-formalism for the Bose-Hubbard model in equilibrium. In order to derive an effective Hamiltonian, we start from the full Bose-Hubbard Hamiltonian, where b † i is a bosonic single-particle creation operator on lattice-site i, n i is the particle-number operator, J denotes the tunneling amplitude, U the on-site interaction, µ the chemical potential and i, j means that we sum over nearest neighbors. In order to determine the thermodynamic properties of the system we have to compute the condensate density and the connected Green function defined respectively as with Nambu . The possibility of broken U (1) symmetry forces b j to be expanded around its mean-field value φ = b j (which we take to be site-independent and can always be chosen real) by b j = φ + δb j . If we concentrate on the site at the origin b o , the Hamiltonian can be rewritten as where H ext contains all terms of the Hamiltonian (1) not containing the origin 'o'. The notation i, o means that we sum over the nearest neighbors of o, and z is the coordination number. We wish to separate the full partition function Z = tr e −βH as Z = Z ext Z o . Here, Z ext is the full (and unknown) partition of the system determined by the terms in the Hamiltonian not involving the site 'o'. It is treated as an irrelevant number in the rest of the paper. The partition function Z o contains the full local hamiltonian H o as well as the correlations introduced by 'ext' on the origin as follows We approximate the expectation value ∆H Hext by the cumulant expansion where ∆H and ∆H∆H are rewritten in terms of Nambu operators and ∆ is an unknown 2x2 real-valued matrix with two independent components ∆ 11 = ∆ 22 and ∆ 12 = ∆ 21 which describes a correction to the common mean-field impurity Hamiltonian. The anomalous term ∆ 12 , containing processes of the type δb 2 , is explicitly taken to be finite in this notation, since it is known from the Bogoliubov theory that deep in the superfluid phase it becomes equally important to the normal (diagonal) term ∆ 11 , containing the δb † δb terms. By (6) we arrive at the effective impurity Hamiltonian As can be seen in the appendix, this effective impurity Hamiltonian is equivalent to B-DMFT in the limit The central characteristic of B+U theory is that we demand that the condensate and the (full) Green function of the Bose-Hubbard model evaluated at 'o' for zero (Matsubara) frequency coincide with the one of system (7), i.e., The paradoxical compatibility of Eqs. (10) and (8) is specific for bosonic systems (see also below Eq. (11)). The equations constitute a self-consistency problem, whose solution also fixes the factors ∆ 11 and ∆ 12 . This can be solved in a unique way if ∆[G c , φ] is invertible, or, technically speaking, if the Luttinger-Ward functional is unique. Since the static mean-field (i.e., the decoupling approximation with ∆ 11 = ∆ 12 = 0) is always a solution, it is easy to convince oneself that multiple (local) minima occur (cf. Ref. 19 for a recent discussion). Nevertheless, we have been able to determine the physically correct solution without problem in all parameter regimes (see below). In practice, one uses an iteration scheme to solve the self-consistency problem. The factors ∆ 11 and ∆ 12 follow from the Dyson equation on the impurity site at zero Matsubara frequency, with an unknown self-energy Σ E . The connected Green's function on the impurity site G c,E is calculated through the full Green's function on the lattice by and the Dyson equation of the full lattice with the bare Green's function given by is the dispersion relation of the lattice and ω n = 2π β n are the Matsubara frequencies.
The approximation (8) shows that the B+U theory has the same functional form as the decoupling approximation in the non-broken phase. In that case only ∆ 11 is present but it acts as a shift in chemical potential. For the broken phase, ∆ 11 and µ are combined with different operators. They control the density of the condensed and depleted atoms, whereas ∆ 12 mainly determines the anomalous density. According to the Bogoliubov theory of the weakly interacting Bose gas, the anomalous propagator is equally important (but opposite in sign) as the normal propagator deep in the superfluid phase. In this way, the deep superfluid regime is taken care of appropriately in our formalism. The Mott localization is enabled by the exact treatment of the density fluctuations on the impurity.

III. VARIATION OF THE SELF-ENERGY
In order to solve the impurity problem, we consider the minimization of its grand potential Ω[Σ, φ] with respect to its selfenergy Σ, as is also done in self-energyfunctional theory 20 and VCA 12 . The minimum with respect to the kinetic condensate term zJφ, δΩ δ(zJφ) = δΩ δ(zJφ * ) = 0, is already taken care of by the selfconsistency condition (9). φ is thus kept constant during the variational calculation of the self-energy. We therefore have to minimze δΩ δΣ = δΩ δ∆ 11 δ∆ 11 δΣ + δΩ δ∆ 12 We are able to find an analytic expression of δΩ δ∆ij , since the grand potential is defined as Ω(∆, φ) = −lnZ(∆, φ) with Z(∆, φ) = Tr[e −βH E (∆,φ) ], giving us After integration, this gives us the relation with some irrelevant integration constant C and In order to avoid unphysical results, we have to introduce upper bounds on |∆ ij |. From (6) we see that ∆ cannot exceed the kinetic energy of a double hopping process of depleted particles, Furthermore, we require that for all momenta k G 11 (ω n = 0, k) −1 ≤ − and det[G(ω n = 0, k) −1 ] ≥ (where a small is introduced for stability requirements when inverting the 2x2 matrix G(ω n = 0, k) −1 in (13)) giving us additional bounds on the self-energy IV. FULL SCHEME AND OBSERVABLES By combining sections II and III we can write down the full self-consistent scheme for the B+U theory. Starting from an initial guess for φ and ∆ (usually the converged mean-field values for ∆ = 0), we calculate a new value for φ through equations (7) and (9). Then we search for the optimal value of the self-energy by minimizing (17) while keeping φ constant. This is done by varying Σ within the bounds (21) and (22) and calculating ∆(Σ) by equations (11)- (13), keeping in mind the bound on ∆ (20). Once the optimal value Σ opt is found, the new value for ∆, ∆(Σ opt ), is plugged into (7), from which a new value for φ is calculated. This procedure is repeated until convergence is reached. In the B+U self-consistency all bonds adjacent to o are included in H E , whereas when computing the quantities per site, all bonds have to be counted only once. In order to calculate the correct thermodynamic quantities per site once convergence is reached, one therefore has to divide ∆ by 2, giving us e.g. for the density per site n = N/V n → n o H E (∆/2,φ) .
We can further divide the Hamiltonian into a kinetic (upper line in (7)) and a potential term (lower line in (7)), giving us expressions for the kinetic and potential energy per site where the total energy per site is given by E tot = E kin + E pot . It should further be noted that even though we do not need to calculate G(τ ) explicitly in the solver and we do not include any retardation in our formalism, we can still calculate correlation functions of the kind A(τ )B(0) by or directly in energy space through the eigenvalues of H E .

V. SIMPLE LIMITS
From the relation (20) it is clear that ∆ → 0 as J → 0. Furthermore also as U goes to zero, ∆ vanishes, since δb † o δb o → 0. Therefore in both cases the mean-field limit is recovered. In the case of U J the mean-field limit is consistent with the weakly interacting Bose gas theory 13 , where the self-energy is frequency-independent as is the case for B+U in our approach. Another simple limit of B+U is the Bethe lattice for a semi-circular density of states given by as was also implemented for B-DMFT 6,7 , which reduces the self-consistency of B+U to one single equation As can be seen in section VI for the 3d case this leads to good agreement with the full self-consistency of a cubic lattice with a much lower numerical cost, since the minimzation routine described in section III is no longer necessary.

VI. RESULTS
In Figure 1 the phase diagram at zero temperature is shown for both a three-dimensional cubic and a twodimensional square lattice. We compare the results with mean-field theory, path integral Monte Carlo simulations with worm-type updates (QMC) from Ref. 21 and B-DMFT results from Ref. 6. The results are identical with the B-DMFT results and agree within a percent with the QMC data both for the 3d and the 2d case. The For simplicity, the B-DMFT results are not shown here, since they overlap with the QMC data within 1% 6 . c) Temperaturedependent phase diagram of a 3d cubic lattice with chemical potential µ/U = 0.4 calculated with B+U (red dots), compared with mean-field (gray line), QMC 21 (white boxes, mostly overlapping with B-DMFT) and B-DMFT 6 (black triangles) results. The B+U results for a semi-circular density of states (Bethe lattice, z=6) are shown as a black dashed line. The systematic error bar is smaller than the size of the dots.
results for a Bethe lattice with coordination number z = 6 are shown as black dashed lines. As can be seen, for the 3d case the simplified self-consistency for the Bethe lattice works very well, showing deviations only near the tip of the Mott lobe. In Figure 1 c) the temperaturedependent phase diagram for µ/U = 0.4 is shown and compared to B-DMFT, QMC and mean-field results for a 3d cubic lattice. In this case the lack of retardation in the B+U formalism leads to a bigger deviation from the B-DMFT results. However, the B+U results are still far more precise than the ones obtained in static mean-field theory. In Figure 2 the temperature dependence of the superfluid order parameter φ, the kinetic energy E kin and the total energy E tot is shown and compared to meanfield and QMC for µ/U = 0.4 ( n ≈ 1) and U/J = 20. It should be noted that, since the optimization in Σ is very sensitive close to the phase transition, we cannot and wish not to make any statements with respect to the order of the phase transition in Figure 2: a local theory such as B+U should be judged for its accuracy on local observables and is by construction unable to capture long wavelength physics. Information on critical phenomena is hence outside its realm of applicability. The kinetic energy is very accurate for low temperatures, but in the  normal phase we find a plateau just as in the decoupling approximation. The corresponding local density n per site is shown in Figure 3 for the same parameters, where also the full local Green's function on the impurity in (26) is shown for the same parameters and temperature T /J = 1.

VII. CONCLUSION
We have presented the B+U-framework for equilibrium studies of the Bose-Hubbard model. It captures the low energy physics of a condensed phase as well as a Mott localization transition when density fluctuations are strongly suppressed on a lattice. The thermodynamics of local quantities of the Bose-Hubbard model can be accurately reproduced everywhere in parameter space by just three parameters φ, Σ 11 and Σ 12 . By treating the self-energy as a variational parameter which minimizes the grand potential, B+U can reproduce both the 3d and 2d phase-transition from the Mott to the superfluid phase at zero temperature with an accuracy of about 1% near the tip of the lobe and better elsewhere. The B+U method can also be applied to finite-temperature systems showing a strong improvement on the simple mean-field limit but it is less accurate than B-DMFT in locating the phase transition line. Just as B-DMFT, being a local theory, it can of course never capture hydrodynamics. Due to its simplicity and low computational cost, B+U is a powerful tool which in the future could be extended to clusters in order to study inhomogeneous or topologically non-trivial bosonic systems, combining exact interactions on clusters, embedded self-consistently in a lattice with any dispersion.

Appendix: Relation between B+U and B-DMFT
We start from the Hamiltonian (1) but instead of Z = tr e −βH we treat the partition function in the presence of retardation as Z = tr e −S with the full action S given by As in section II we expand b j (τ ) = φ + δb j (τ ) around its site-and imaginary-time-independent mean-field value φ, giving us As with the Hamiltonian in II we approximate the expectation value ∆S Sext by the cumulant expansion leading to the effective B-DMFT impurity action 6,7 By the relation If we now take the zero Matsubara frequency of this equation we recover the B+U Dyson equation on the impurity (11).