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

###### Abstract

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 phase diagrams of the three-dimensional and two-dimensional Bose-Hubbard model with an accuracy of or better. The superfluid to normal transition at finite temperature is also reproduced well and only slightly less accurately than in B-DMFT.

###### pacs:

71.10.Fd, 02.70.Ss, 05.30.Jp## 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 [Bloch et al., 2008]. Dramatic experimental progress in this field, such as the observation of the Mott insulator to superfluid phase transition in the Bose-Hubbard model [Bloch et al., 2008] or the recent realization of the Hofstadter model [Aidelsburger et al., 2013], 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. In one dimension the density matrix renormalization group (DMRG) [Kühner and Monien, 1998; Rapsch et al., 1999; Kollath et al., 2004, 2007] works very well, while in higher dimensions the numerical complexity represents a problem. Very successful also in higher dimensions have been path integral quantum Monte Carlo (QMC) simulations with worm-type updates [Prokof’ev et al., 1998] leading to identical results for the Bose-Hubbard model as observed in experiment [Trotzky et al., 2010,Pollet, 2012]. Despite all its impressive successes for bosonic systems the worm-algorithm 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 mean-field theory [Anders et al., 2011,Anders et al., 2010] (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 in three dimensions (3D)] with experimental and QMC data [Anders et al., 2011] 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 multi-component systems and systems with long-range interactions and/or dispersions [He et al., 2014,He et al., 2012].

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 [Anders et al., 2011,Anders et al., 2010], where, due to the difference in sign between the normal and the anomalous Green’s function, a sign problem arises in the symmetry broken phase.

In this paper, 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 high-energy community [Akerlund and de Forcrand, 2013,Akerlund et al., 2014] 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 , , and (the condensate order parameter, and the normal and the anomalous self-energy, respectively) 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 (VCA) by also considering non zero values of pair creation and annihilation of depleted particles [Knap et al., 2011]. It is also the simplest accurate extension of the weakly interacting Bose gas theory [Capogrosso-Sansone et al., 2010] 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. [Kleinert et al., 2014] and the two collective fields proposed by Cooper et al. [Dawson et al., 2013; Cooper et al., 2012; Mihaila et al., 2011], and behaves quantitatively much better. It has a similar functional degree of freedom as the projector technique introduced by Trefzger and Sengupta [Trefzger and Sengupta, 2011] for finite lattices. Among other methods which are used to treat the Bose-Hubbard model are also the process chain approach [Eckardt, 2009,Teichmann et al., 2009] and the quantum phase field rotor field [Polak and Kopeć, 2009].

The paper is organized as follows. In Sec. II the B+U formalism is derived for the Bose-Hubbard model in equilibrium, while in Sec. III the details of the variational 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 Sec. IV, while some simple limits of B+U are explained in Sec. V. In Sec. VI we present the results at zero and finite temperature comparing them with QMC and B-DMFT. Finally, in Sec. 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,

(1) |

where is a bosonic single-particle creation operator on lattice-site , is the particle-number operator, denotes the tunneling amplitude, the on-site interaction, the chemical potential, and 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’s function, defined respectively as

(2) | |||||

(3) |

with Nambu notation . The possibility of broken symmetry forces to be expanded around its mean-field value (which we take to be site independent and can always be chosen real) by . If we concentrate on the site at the origin , the Hamiltonian can be rewritten as

(4) | |||||

where contains all terms of the Hamiltonian (1) not containing the origin “”. The notation means that we sum over the nearest neighbors of , and is the coordination number. We wish to separate the full partition function as . Here, is the full (and unknown) partition of the system determined by the terms in the Hamiltonian not involving the site . It is treated as an irrelevant number in the rest of the paper. The partition function contains the full local Hamiltonian as well as the correlations introduced by “ext” on the origin as follows,

(5) |

We approximate the expectation value by the cumulant expansion

(6) | |||||

where and are rewritten in terms of Nambu operators and is an unknown real-valued matrix with two independent components and which describes a correction to the common mean-field impurity Hamiltonian. The anomalous term , containing processes of the type , 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 , containing the terms. By (6) we arrive at the effective impurity Hamiltonian

(7) | |||||

As can be seen in the Appendix, this effective impurity Hamiltonian is equivalent to B-DMFT in the limit

(8) |

Since is independent on imaginary time, the Dyson equation on the impurity has to be evaluated only for a single Matsubara value. The Green’s function which mirrors the symmetry relations assumed for is the one evaluated at . The central characteristic of B+U theory is that we demand that the condensate and the (full) Green’s function of the Bose-Hubbard model evaluated at o for zero (Matsubara) frequency coincide with the one of system (7), i.e.,

(9) | |||||

(10) |

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 and . This can be solved in a unique way if is invertible, or, technically speaking, if the Luttinger-Ward functional is unique. Since the static mean-field (i.e., the decoupling approximation with ) is always a solution, it is easy to convince oneself that multiple (local) minima occur (cf. Ref. [Kozik et al., 2015] 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 and follow from the Dyson equation on the impurity site at zero Matsubara frequency,

(11) |

with an unknown self-energy . The connected Green’s function on the impurity site is calculated through the full Green’s function on the lattice by

(12) |

and the Dyson equation of the full lattice

(13) |

with the bare Green’s function given by , where is the dispersion relation of the lattice and 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 is present but it acts as a shift in chemical potential. For the broken phase, and are combined with different operators. They control the density of the condensed and depleted atoms, whereas 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 self-energy , as is also done in self-energy-functional theory [Potthoff, 2012] and VCA [Knap et al., 2011]. The minimum with respect to the kinetic condensate term , , is already taken care of by the self-consistency condition (9). is thus kept constant during the variational calculation of the self-energy. We therefore have to minimze

(14) |

We are able to find an analytic expression of , since the grand potential is defined as with , giving us

(15) | |||||

(16) |

After integration, this gives us the relation

(17) |

with some irrelevant integration constant and

(18) | |||||

(19) |

In order to avoid unphysical results, we have to introduce upper bounds on . From (6) we see that cannot exceed the kinetic energy of a double hopping process of depleted particles,

(20) |

Furthermore, we require that for all momenta and [where a small is introduced for stability requirements when inverting the x matrix in (13)], giving us additional bounds on the self-energy

(21) | |||||

(22) |

## Iv Full scheme and Observables

By combining Secs. 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 ), we calculate a new value for through Eqs. (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 Eqs. (11)-(13), keeping in mind the bound on (20). Once the optimal value is found, the new value for , , 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 are included in , 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 ,

(23) |

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

(24) | |||||

(25) |

where the total energy per site is given by . It should further be noted that even though we do not need to calculate explicitly in the solver and we do not include any retardation in our formalism, we can still calculate correlation functions of the kind by

(26) |

or directly in energy space through the eigenvalues of . Since the B+U solver consists of a single impurity, in order to compute momentum-dependent quantities, one has to resort to the approximate expression

(27) |

with

(28) |

By a Fourier transformation this enables us to compute such quantities as the momentum-dependent density

(29) |

or the critical quasi particle and quasi hole energies at zero momentum which can be evaluated from the asymptotic behavior of at zero temperature through [Capogrosso-Sansone et al., 2007,Duchon et al., 2013]

(30) |

where .

## V Simple Limits

From the relation (20) it is clear that as . Furthermore, also, as goes to zero, vanishes, since . Therefore in both cases the mean-field limit is recovered. In the case of the mean-field limit is consistent with the weakly interacting Bose gas theory [Capogrosso-Sansone et al., 2010], 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 semicircular density of states given by

(31) |

as was also implemented for B-DMFT [Anders et al., 2011,Anders et al., 2010] , which reduces the self-consistency of B+U to one single equation,

(32) |

As can be seen in Sec. 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 minimization routine described in Sec. III is no longer necessary.

## Vi Results

In Fig. 1 the phase diagram at zero temperature is shown for both a three-dimensional (3D) cubic and a two-dimensional (2D) square lattice. We compare the results with mean-field theory, path integral Monte Carlo simulations with worm-type updates (QMC) from Ref. [Capogrosso-Sansone et al., 2007], and B-DMFT results from Ref. [Anders et al., 2011]. 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 cases. The results for a Bethe lattice with coordination number 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 Fig. 1(c) the temperature-dependent phase diagram for 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 Fig. 2 the temperature dependence of the superfluid order parameter , the kinetic energy , and the total energy is shown and compared to mean-field and QMC for () and . 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 Fig. 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 per site is shown in Fig. 3 for the same parameters, where also the full local Green’s function on the impurity in imaginary time calculated by (26) is shown for the same parameters and temperature . In Fig. 4(a) we plot the density in momentum-space for and different values of at zero temperature, where are the depleted particles calculated from the connected Green’s function and is the condensate fraction. In Fig. 4(b) the quasi-particle and quasi-hole energies in the Mott phase at zero momentum are extrapolated from the imaginary-time dependence of the zero-momentum Green’s function through (30) and compared to QMC results from Ref. [Duchon et al., 2013] and the analytic zeroth-order solution for . We see both good agreement with QMC for finite as with the strongly interacting limit for

## 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 , , and . 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 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.

###### Acknowledgements.

We wish to thank M. Eckstein, H. U. R. Strand, and F. A. Wolf for fruitful discussions. This work was supported by FP7/ERC Starting Grant No. 306897 and FP7/Marie-Curie Grant No. 321918.*

## Appendix A Relation between B+U and B-DMFT

We start from the Hamiltonian (1), but instead of we treat the partition function in the presence of retardation as with the full action given by

(33) |

As in Sec. II we expand around its site- and imaginary-time-independent mean-field value , giving us

(34) | |||||

and expand the full partition function as by

(35) |

As with the Hamiltonian in Sec. II we approximate the expectation value by the cumulant expansion,

(36) | |||||

leading to the effective B-DMFT impurity action [Anders et al., 2011,Anders et al., 2010]

(37) | |||||

By the relation , one can see that the action (37) is equivalent to the effective Hamiltonian (7) for the ansatz

(38) |

In fact, (38) reduces the Dyson equation of B-DMFT on the impurity to

(39) |

If we now take the zero Matsubara frequency of this equation we recover the B+U Dyson equation on the impurity (11).

## References

- Bloch et al. (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. of Mod. Phys. 80, 885 (2008).
- Aidelsburger et al. (2013) M. Aidelsburger, M. Atala, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch, Phys. Rev. Lett. 111, 185301 (2013).
- Kühner and Monien (1998) T. D. Kühner and H. Monien, Phys. Rev. B 58, R14741 (1998).
- Rapsch et al. (1999) S. Rapsch, U. Schollwöck, and W. Zwerger, Europhys. Lett. 46, 559 (1999).
- Kollath et al. (2004) C. Kollath, U. Schollwöck, J. von Delft, and W. Zwerger, Phys. Rev. A 69, 031601(R) (2004).
- Kollath et al. (2007) C. Kollath, A. M. Läuchli, and E. Altman, Phys. Rev. Lett. 98, 180601 (2007).
- Prokof’ev et al. (1998) N. Prokof’ev, B. Svistunov, and I. Tupitsyn, J. Exp. Theor. Phys. 87, 310 (1998).
- Trotzky et al. (2010) S. Trotzky, L. Pollet, F. Gerbier, U. Schnorrberger, I. Bloch, N. Prokof’ev, B. Svistunov, and M. Troyer, Nat. Phys. 6, 998 (2010).
- Pollet (2012) L. Pollet, Rep. Prog. Phys. 75, 094501 (2012).
- Anders et al. (2011) P. Anders, E. Gull, L. Pollet, M. Troyer, and P. Werner, New J. Phys. 13, 075013 (2011).
- Anders et al. (2010) P. Anders, E. Gull, L. Pollet, M. Troyer, and P. Werner, Phys. Rev. Lett. 105, 096402 (2010).
- He et al. (2014) L. He, A. Ji, and W. Hofstetter, (2014).
- He et al. (2012) L. He, Y. Li, E. Altman, and W. Hofstetter, Phys. Rev. A 86, 043620 (2012).
- Akerlund and de Forcrand (2013) O. Akerlund and P. de Forcrand, (2013).
- Akerlund et al. (2014) O. Akerlund, P. de Forcrand, A. Georges, and P. Werner, Phys. Rev. D 90, 065008 (2014).
- Knap et al. (2011) M. Knap, E. Arrigoni, and W. von der Linden, Phys. Rev. B 83, 134507 (2011).
- Capogrosso-Sansone et al. (2010) B. Capogrosso-Sansone, S. Giorgini, S. Pilati, L. Pollet, N. Prokof’ev, B. Svistunov, and M. Troyer, New J. Phys. 12, 043010 (2010).
- Kleinert et al. (2014) H. Kleinert, Z. Narzikulov, and A. Rakhimov, J. Stat. Mech. , P01003 (2014).
- Dawson et al. (2013) J. F. Dawson, F. Cooper, C.-C. Chien, and B. Mihaila, Phys. Rev. A 88, 023607 (2013).
- Cooper et al. (2012) F. Cooper, C.-C. Chien, B. Mihaila, J. F. Dawson, and E. Timmermans, Phys. Rev. A 85, 023631 (2012).
- Mihaila et al. (2011) B. Mihaila, F. Cooper, J. Dawson, C.-C. Chien, and E. Timmermans, Phys. Rev. A 84, 023603 (2011).
- Trefzger and Sengupta (2011) C. Trefzger and K. Sengupta, Phys. Rev. Lett. 106, 095702 (2011).
- Eckardt (2009) A. Eckardt, Phys. Rev. B 79, 195131 (2009).
- Teichmann et al. (2009) N. Teichmann, D. Hinrichs, M. Holthaus, and A. Eckardt, Phys. Rev. B 79, 224515 (2009).
- Polak and Kopeć (2009) T. P. Polak and T. K. Kopeć, J. Phys. B 42, 095302 (2009).
- Kozik et al. (2015) E. Kozik, M. Ferrero, and A. Georges, Phys. Rev. Lett. 114, 156402 (2015).
- Potthoff (2012) M. Potthoff, in Strongly Correlated Systems, Springer Series Solid State Physics, Vol. 171, edited by A. Avella and F. Mancini (Springer, Berlin, Heidelberg, 2012) Chap. 10, pp. 303–339.
- Capogrosso-Sansone et al. (2007) B. Capogrosso-Sansone, N. Prokof’ev, and B. Svistunov, Phys. Rev. B 75, 134302 (2007).
- Duchon et al. (2013) E. Duchon, Y. L. Loh, and N. Trivedi, in Novel Superfluids, Vol. 2, edited by J. Ketterson and K.-H. Benneman (Oxford University Press, Oxford, UK, 2013) p. 193.