Exact solution of the Bose-Hubbard model on the Bethe lattice

# Exact solution of the Bose-Hubbard model on the Bethe lattice

Guilhem Semerjian LPTENS, Unité Mixte de Recherche (UMR 8549) du CNRS et de l’ENS, associée à l’UPMC Univ Paris 06, 24 Rue Lhomond, 75231 Paris Cedex 05, France.    Marco Tarzia Laboratoire de Physique Théorique de la Matière Condensée, Université Pierre et Marie Curie-Paris 6, UMR CNRS 7600, 4 place Jussieu, 75252 Paris Cedex 05, France.    Francesco Zamponi LPTENS, Unité Mixte de Recherche (UMR 8549) du CNRS et de l’ENS, associée à l’UPMC Univ Paris 06, 24 Rue Lhomond, 75231 Paris Cedex 05, France.
###### Abstract

The exact solution of a quantum Bethe lattice model in the thermodynamic limit amounts to solve a functional self-consistent equation. In this paper we obtain this equation for the Bose-Hubbard model on the Bethe lattice, under two equivalent forms. The first one, based on a coherent state path integral, leads in the large connectivity limit to the mean field treatment of Fisher et al. [Phys. Rev. B 40, 546 (1989)] at the leading order, and to the bosonic Dynamical Mean Field Theory as a first correction, as recently derived by Byczuk and Vollhardt [Phys. Rev. B 77, 235106 (2008)]. We obtain an alternative form of the equation using the occupation number representation, which can be easily solved with an arbitrary numerical precision, for any finite connectivity. We thus compute the transition line between the superfluid and Mott insulator phases of the model, along with thermodynamic observables and the space and imaginary time dependence of correlation functions. The finite connectivity of the Bethe lattice induces a richer physical content with respect to its infinitely connected counterpart: a notion of distance between sites of the lattice is preserved, and the bosons are still weakly mobile in the Mott insulator phase. The Bethe lattice construction can be viewed as an approximation to the finite dimensional version of the model. We show indeed a quantitatively reasonable agreement between our predictions and the results of Quantum Monte Carlo simulations in two and three dimensions.

## I Introduction

The bosonic version of the Hubbard model was introduced for the first time in the seminal paper FiFi (). It describes a system of bosons hopping between neighboring sites of a lattice, and subjected to a local repulsive interaction disfavouring the multiple occupancy of a site. The competition between these two effects leads to a quantum phase transition Sachdev () at zero temperature: when the hopping term dominates the groundstate is a Bose-Einstein Condensate (BEC) delocalized over the lattice, also known as a superfluid (SF) phase. On the contrary for large local repulsion it becomes energetically favorable to form a commensurate state where the average number of bosons per site is fixed to an integer value. This incompressible phase is called a Mott insulator (MI).

The introduction of this model was motivated by experimental results on Helium adsorbed on disordered substrates. The recent progresses of the experimental techniques for the manipulation of cold atoms, and in particular the possibility of devising optical lattices loaded with cold gases, revivified the interest for the Bose-Hubbard model. Following the proposal made in bh_proposal_exp (), the experimental observation of this Mott transition was first achieved in bh_exp (). We refer the reader to review_exp () for a review of such experiments bridging a gap between atomic and condensed-matter physics.

From a more theoretical point of view the Bose-Hubbard model has been studied with various techniques, notably mean-field approximations FiFi (); Gutzwiller (), perturbative expansions in the hopping strength bh_perturb (); bh_perturb2 (); freericks1 (), Random Phase Approximations (RPA) Kampf93 (); Sheshadri93 (); vanOosten01 (); Sengupta05 (); Menotti08 (), and Quantum Monte Carlo simulations bh_num (); bs_num (); bh_qmc_2d (); bh_qmc_3d (); bh_qmc_number (); hohenadler (). The mean-field approach of FiFi () yields a qualitatively correct prediction of the phase diagram of the model. However its description of the Mott insulator is oversimplified, the hopping being completely neglected in this phase at the mean-field level. As often the assumptions underlying the mean-field approximation become valid when the connectivity of each site (the number of its neighbors) is very large, either of the order of the system size itself (in which case a rigorous treatment of the model is possible bh_rig_1 (); bh_rig_2 ()), or in the limit of large dimensionality. The latter case was recently investigated in BDMFT_1 (); BDMFT_2 () where a bosonic version of the Dynamic Mean-Field Theory DMFT () was developed, including first order corrections in the inverse of the spatial dimension.

In this paper we follow a related but slightly different road to improve on the mean-field treatment of FiFi (), by treating the Bose-Hubbard model at the level of the Bethe approximation, i.e. by solving it on a Bethe lattice. We show that, in the thermodynamic limit, the computation of all observables amounts to solve a single self-consistent functional equation. By writing this equation in the coherent state basis one recovers the B-DMFT in the large connectivity limit BDMFT_1 (); BDMFT_2 (). Unfortunately, this equation cannot be analytically solved for finite connectivity, except in some special limits. Hence, we show how to rewrite it in a more convenient way using the occupation number basis; this opens the way to its resolution through an efficient numerical algorithm.

The Bethe approximation is the next step in a hierarchy of approximations for lattice systems known as the Cluster Variation Method CVM (); it is in fact exact when the model considered is defined on a Bethe lattice, that is a graph which is locally a tree (no short loops can be closed on it). Bethe lattices appear in the DMFT derivations DMFT (); BDMFT_1 (); BDMFT_2 () as an intermediate step, before taking the limit of large connectivity. On the contrary in the present paper the connectivity of the Bethe lattice will be kept finite. From the point of view of the universality classes of critical exponents Bethe lattice models fall in the mean-field category, yet their finite connectivity makes them richer than the (fully-connected) mean-field version of FiFi (). In the latter any site is adjacent to all others, whereas on a Bethe lattice one has a well defined notion of distance between two sites, though it does not correspond to the usual Euclidean distance. Moreover the MI phase is non-trivial for the Bethe lattice version of the model, at variance with the mean-field picture. Though we only consider here the ordered version of the model a motivation for studying the Bethe lattice is the possibility that it exhibits a Bose Glass phase FiFi (); inguscio1 () in presence of disorder, which cannot be described at the simplest mean-field level. This hope is backed up by the relevance of the Bethe lattices for the localization problem localization_Bethe ().

The methodology we use has its roots in the extensive research effort which took place over the last decade in the community studying statistical mechanics of disordered systems. Under the name of cavity method a set of techniques has been developed to solve classical spin models on Bethe lattices (in the sense of sparse random graphs), with applications to spin glasses MePa_Bethe () and to random optimization problems arising from theoretical computer science ksat_science (). An extensive presentation of this field can be found in the recent book cavity_book (). The cavity method has been recently extended from classical to quantum spin models in presence of a transverse field qcav_scardi (); qcav_knysh (); qcav_our (). We widen here the range of applicability of the method by including lattice boson models.

The rest of the paper is organized as follows. In Sec. II we define precisely the model studied and the physical observables of interest, and we recall the usual mean-field treatment and its qualitative physical predictions. Sec. III is an overview conveying the main ideas of the cavity method, starting for pedagogical reasons from the ferromagnetic Ising model before explaining its extension to quantum problems, and in particular its connection with the B-DMFT of BDMFT_1 (); BDMFT_2 (). Sec. IV contains the core of our contribution; we first present in Sec. IV.1 the main equations describing the Bose-Hubbard model on the Bethe lattice, and the principles of their numerical resolution. In Sec. IV.2 we present our predictions for several physical observables and compare the phase diagram we obtain with the Quantum Monte Carlo results in two bh_qmc_2d () and three bh_qmc_3d () dimensions, as well as with the B-DMFT prediction of BDMFT_2 () in three dimensions. For the convenience of the reader not interested in technical aspects we postpone the detailed derivation of the equations to Sec. V. Finally we present our conclusions and perspectives for future work in Sec. VI.

## Ii Definitions and reminder on the traditional mean-field approach

### ii.1 Definition of the model and of the observables

We shall consider the Bose-Hubbard model for spinless bosons on a graph of vertices (sites) defined by the following Hamiltonian,

 H=−J∑⟨i,j⟩(a†iaj+a†jai)+U2N∑i=1a†ia†iaiai−μN∑i=1a†iai=HK+HL , (1)

where the first sum runs over a subset of the possible edges (links) between pairs of sites, is the hopping amplitude between neighboring sites, (resp. ) is the boson annihilation (resp. creation) operator on site , is the chemical potential fixing the density of particles, and controls the strength of the on-site interaction between particles. It is convenient to separate the kinetic term (proportional to ) in the Hamiltonian, from the local term including the Hubbard interaction and the chemical potential term. Note that the kinetic energy defined in such a way will be negative: indeed, the discrete lattice version of the kinetic energy would be given by where is the connectivity of site .

Although the method we will discuss in this paper allows in principle to compute very general observables (such as multi-point imaginary time correlations), in the following we will be mainly interested in standard observables such as the mean density, kinetic and on-site energy, condensate fraction, and Green functions. For the sake of clarity we recall now their definitions. The partition function at temperature is defined by

 Z=Tr[e−βH] , (2)

where (we set ), the corresponding free energy is ; the free-energy per site is . The thermodynamic average of an operator is defined as

 ⟨O⟩=1ZTr[Oe−βH] . (3)

In particular the average number of particles on site reads and the density of particles is ; the kinetic (resp. local) energy per site is (resp. ). To define the order parameter for BEC, one possibility that is convenient for our purposes is to introduce a small perturbation to the Hamiltonian in the form and define

 ⟨a⟩=limε→0⟨a⟩ε , (4)

where denotes the thermodynamic average in presence of the perturbation.

The imaginary-time evolution of an operator is given by ; we then define the two-times correlation functions as Negele ()

 Gi>(τ)=⟨ai(τ)a†i(0)⟩=1ZTr[e−(β−τ)Haie−τHa†i] ,Gi<(τ)=⟨a†i(0)ai(τ)⟩=1ZTr[e−(β+τ)Ha†ieτHai] . (5)

The Green function is defined, for , by

 (6)

where is the usual time-ordering operation, which should not be confused with the temperature. Note that the cyclic property of the trace implies that , hence

 Gi(τ)=θ(τ)Gi>(τ)+θ(−τ)Gi>(β+τ) , (7)

and the knowledge of for is enough to determine the Green function.

The cavity method allows also the computation of spatial correlation functions, we shall in particular determine the one-particle density matrix .

### ii.2 Review of mean-field results

The inexistence of an analytical solution of the model for finite dimensional lattices has triggered a large amount of research effort on numerical simulations bh_num (); bs_num (); bh_qmc_2d (); bh_qmc_3d (); bh_qmc_number (), pertubative expansions bh_perturb (); bh_perturb2 (), or mean-field treatments FiFi (); Gutzwiller () of the problem. Let us briefly recall the various points of view and methodology that the mean-field approach usually encompasses. It can first be seen as an approximation to finite-dimensional models. Maybe the most satisfactory way to perform this approximation is by means of a variational method. One introduces a simpler trial (Gutzwiller) Hamiltonian where the sites are decoupled,

 H0=N∑i=1[−ψia†i−ψ∗iai+U2a†ia†iaiai−μa†iai] , (8)

the being here free parameters. Note that breaks the particle conservation symmetry. The free-energy at inverse temperature of the original model, , can be upper-bounded as , where is the free-energy of the trial Hamiltonian and denotes a thermal average with respect to the trial Hamiltonian. The best variational description is thus obtained by minimizing the upper bound with respect to the parameters . Assuming that all sites of the graph have the same number of neighbors ( for a -dimensional hypercubic lattice), the bound is minimized by taking the same value on all sites, which can be chosen real without loss of generality. The free-energy per site can then be upper-bounded as , with

 fvar(β)=infψ[1cJψ2−1βlnTr[eβμa†a−βU2a†a†aa+βψ(a+a†)]] , (9)

where the trace is over the Hilbert space of a single site. The physical properties of this variational free-energy are well-known FiFi (): at all temperatures and chemical potentials there is a transition encountered upon increasing the hopping intensity from an “insulating” phase, characterized by , to a “superfluid” phase with . At zero temperature this transition line draws a series of lobes in the plane, inside each lobe the particle density being constrained to a given integer. These Mott Insulator phases are thus incompressible. There exists alternative ways to obtain the estimation (9) for the free-energy per site of the Bose-Hubbard Hamiltonian. The usual mean-field approximation consists indeed in replacing with in the hopping term of the original Hamiltonian, neglecting correlations between neighboring sites. This transforms the Hamiltonian into the site-decoupled form (8), with given by a sum of over the neighbors of . These expectation values are then computed with respect to the decoupled Hamiltonian, the self-consistency equations leading finally to the same expression (9) of the free-energy per site as the variational approach; the latter has however the advantage of being better controlled, in the sense that it provides a rigorous bound on the true free-energy of the system. Finally, another reasoning yielding this mean-field result consists in devising a model which has exactly (9) as its free-energy per site in the thermodynamic limit, instead of taking it as an approximation for the finite-dimensional case. As could be expected this model corresponds to the fully-connected version of the Hamiltonian (1), with the sum in the hopping term running over all possible pairs of sites, with a coupling constant inversely proportional to the size of the system in order to have an extensive thermodynamic limit. It has been shown rigorously in bh_rig_1 (); bh_rig_2 () that in the thermodynamic limit the free-energy of this fully-connected model converges indeed to (9).

The above described mean-field treatment has limitations both of a quantitative nature (the approximation cannot be expected to be very precise for small dimensions) and of a qualitative one. In particular the MI phase is rather trivial; as the infimum in (9) is reached in , the hopping of particles is completely suppressed in this phase. This drawback is particularly severe in the case of the disordered Bose-Hubbard model, for which it forbids the description of the Bose Glass phase predicted in FiFi (). In order to cure this defect of the mean-field treatment one has to account in some way for the correlations between neighboring sites. Thinking in terms of classical spin models, the mean-field approximation is the lowest level of a hierarchy of descriptions (known as Cluster Variation Method CVM (), or Kikuchi approximations Kikuchi ()) which treats exactly larger and larger regions of the interaction graph. In this work we shall deal with the Bose-Hubbard model at the next level of the hierarchy, known as the Bethe approximation, in which correlations between nearest neighbors are explicitly taken into account. In the same way as the simplest mean-field approximation was exact for the fully-connected model, the Bethe approximation is exact for a family of graphs, known as Bethe lattices. In these graphs each site interacts with a finite number of neighbors, say in order to mimic a -dimensional hypercubic lattice, but the short loops present in the latter are discarded: Bethe lattices have a local tree structure, see Fig. 1 for a picture of the local appearance of a Bethe lattice of connectivity . For at least two reasons it is however better not to picture a Bethe lattice as a finite regular tree (usually called Cayley tree in this context). Indeed a regular tree is strongly inhomogenous, a finite fraction of its “volume” being very close to its “surface”, and only the center of the Cayley tree has the bulk properties one is interested in. Moreover in the case of frustrated models the boundary conditions imposed on the leaves of the Cayley tree have to be chosen with particular care. For these two reasons it has become customary in the community of statistical mechanics of disordered systems, following MePa_Bethe (), to define a Bethe lattice as a random -regular graph rgraphs (), that is a graph chosen uniformly at random from the set of the graphs on vertices where all sites have precisely neighbors. These Bethe lattices are locally tree-like, their loops have typically a length diverging with the size of the system, yet they do not have any boundary, all sites playing the same role (in the same way as periodic boundary conditions impose translation invariance on a finite cubic lattice). The absence of an underlying finite-dimensional structure gives them a mean-field character, but their finite connectivity leads to a richer content than the fully-connected models. The free-energy of Bethe lattice models is self-averaging with respect to their random character in the thermodynamic limit. In other words for large enough a single sample is a good representative of the ensemble average. The so-called cavity method has been developed to treat classical spin models defined on such random graphs, and found important applications for optimization problems of theoretical computer science ksat_science ().

## Iii Overview of the quantum cavity method

### iii.1 The cavity method for ferromagnetic Ising models

For the sake of pedagogy we shall begin our presentation by the cavity method treatment of the ferromagnetic Ising model on the Bethe lattice. We consider the following energy function,

 H(σ1,…,σN)=−J∑⟨i,j⟩σiσj−h∑iσi , (10)

where the classical degrees of freedom can take the values , the first term describing ferromagnetic () interactions between neighbors on a graph of vertices, the second one corresponding to an uniform magnetic field. We shall repeatedly use in the following the notation for the set of vertices adjacent to a given vertex , i.e. for the sites which interact with , and for those vertices around distinct from . The Gibbs-Boltzmann probability measure is

 η(σ1,…,σN)=1Ze−βH(σ1,…,σN) , (11)

with the partition function ensuring its normalization. The goal is to compute the free-energy per site and the local magnetizations , the brackets denoting an average with respect to the law (11).

Let us first consider the case where the interaction graph is a finite tree. In this case the computation can be organized in a very simple way, taking benefit of the natural recursive structure of a tree. We define the quantity , for two adjacent sites and , as the partial partition function for the subtree rooted at , excluding the branch directed towards , with a fixed value of the spin variable on the site . We also introduce , the partition function of the whole tree with a fixed value of . These quantities can be computed according to the following recursion rules, see Fig. 2 for an example,

 Zi→j(σi)=eβhσi∏k∈∂i∖j(∑σkZk→i(σk)eβJσiσk) ,Zi(σi)=eβhσi∏j∈∂i(∑σjZj→i(σj)eβJσiσj) . (12)

It will be useful for the following discussion to rewrite these equations in terms of normalized quantities which can be interpreted as probability laws for the random variable , namely and . The recursion equations read in these notations

 ηi→j(σi)=1zi→jeβhσi∏k∈∂i∖j(∑σkηk→i(σk)eβJσiσk) ,ηi(σi)=1zieβhσi∏j∈∂i(∑σjηj→i(σj)eβJσiσj) , (13)

where and are normalization constants. On a given tree the recursion equations on the for all directed edges of the graph have a single solution, easily found by propagating the recursion from the leaves of the graph. Moreover the quantity is exactly the marginal probability law of the Gibbs-Boltzmann distribution (11), hence the local magnetizations can be computed as .

The reasoning above was made under the assumption that the interaction graph was a tree. One can however look for a solution of the recursion equations (13) on any graph, even in the presence of loops. This approach is known as Belief Propagation in inference problems fgraphs (), and corresponds to the Bethe approximation of statistical mechanics Yedidia (). The cavity method is a set of prescriptions to use these local recursion equations to predict the behavior of models defined on sparse random graphs, that are locally tree-like (the typical length of the loops diverging in the thermodynamic limit). In its simplest version (known as the replica symmetric one) one makes the assumption of the existence of a single pure state which implies that the effect of the loops does not spoil the local recursions (13). Their presence simply provides self-consistent boundary conditions and avoids the ill-defined behaviour due to the dominant surface effects of trees. In the case of unfrustrated ferromagnetic models it has been shown rigorously ferro_rigorous () that this assumption is correct, the predictions of the cavity method being exact in the thermodynamic limit, both for the local magnetizations and for the free-energy per site. Let us be more explicit for the case of the Bethe lattice, where all vertices have the same connectivity . The probability laws are then the same on all edges; denoting their common value, one finds the self-consistent equation

 ηcav(σ)=1zcaveβhσ∑σ1,…,σc−1ηcav(σ1)…ηcav(σc−1)eβJσ(σ1+⋯+σc−1) , (14)

with a normalization constant. A pictorial representation of this equation can be found in Fig. 3. The local magnetization is then computed as

 ⟨σ⟩=∑ση(σ)σ ,η(σ)=1zeβhσ∑σ1,…,σcηcav(σ1)…ηcav(σc)eβJσ(σ1+⋯+σc) , (15)

including the neighbors of a central site as represented in Fig. 4. The term cavity comes from the fact that in Eq. (14) one site has been removed from the neighborhood of the considered vertex. As the Ising variable can only take two values, each of the probability distributions and can be parametrized by a single real number, a cavity (resp. effective) magnetic field,

 ηcav(σ)=eβhcavσ2cosh(βhcav) ,η(σ)=eβheffσ2cosh(βheff) , (16)

solutions of

 hcav=h+c−1βarctanh[tanh(βJ)tanh(βhcav)] ,heff=h+cβarctanh[tanh(βJ)tanh(βhcav)] . (17)

thus making the resolution of (14),(15) extremely simple. In particular at zero external field , one finds that a phase transition occurs at , separating a high temperature paramagnetic phase () from a low temperature ferromagnetic phase (). The critical temperature is easily obtained linearizing the equation for around and is the solution of .

It is worth noting that one can compute explicitly the spin-spin correlation function, which is given in the paramagnetic phase by , where is the distance between sites and , defined as the length of the shortest path connecting sites and . The associated correlation length is which is finite at the transition point . Nevertheless, the associated magnetic susceptibility diverges: one can show using the fluctuation-dissipation theorem that , where is the number of points at distance from a given reference point and scales as for large . Therefore, if decays slower than , the corresponding susceptibility is divergent; this is indeed consistent with the equation for obtained above. Hence phase transitions on Bethe lattices are always associated to diverging susceptibilities and finite correlations lengths (see Farhi () for a discussion of this fact in the context of quantum spin models). For finite dimensional lattices, grows as a power of , and one needs a diverging correlation length to obtain a diverging susceptibility.

### iii.2 Suzuki-Trotter formalism and quantum cavity method

We have just seen how the cavity method deals with classical spin models on locally tree-like graphs. The extension of the method to quantum models is based on the Suzuki-Trotter formula. Generically speaking the resolution of a quantum model amounts to compute its partition function , where the Hamiltonian is now an operator, which can usually be split as with two operators which do not commute. For instance in the case of quantum spin 1/2 models one can have an interaction term involving only the component of the spins, while is a transverse field acting in a perpendicular direction. In the case of the Bose-Hubbard model this decomposition splits the Hamiltonian between the hopping term and the on-site potential energy . The non-commutativity of can be cured with the application of the Suzuki-Trotter formula suzuki (),

 Z=limNs→∞Tr[(e−βNsH1e−βNsH2)Ns] . (18)

The computation then proceeds with the insertion of representations of the identity between the elements of this product. It is convenient to express the identity operator in a basis where one of the two operators is easily evaluated. In the case of quantum spin 1/2 models one can for instance use the eigenvectors of the Pauli matrix in the direction; the quantum partition function is then turned into the partition function of a classical Ising spin model, where each quantum spin is replaced by a set of classical spins, indexed by their position on the original graph and a supplementary “discrete imaginary time” coordinate (which becomes a continuous parameter in the limit). The important point is that the spatial structure of the graph of interactions is preserved in this construction, at the price of the introduction of imaginary time-dependent classical degrees of freedom. In particular, if the interactions of the quantum model are defined on a tree-like graph, the cavity method still applies to this extended classical model. This line of thought was first followed for quantum spin models in qcav_scardi () (see also qcav_knysh ()) for a finite number of Suzuki-Trotter slices, the continuous imaginary time limit was then taken in qcav_our (). In this paper we shall extend this method to deal with lattice bosons models.

In this case the decomposition of the identity operator in the Suzuki-Trotter can be expressed using either coherent states or occupation numbers. The latter has the advantage of being discrete, and we shall use it in the rest of the paper. For the sake of clarity and in order to make contact with the recently proposed B-DMFT BDMFT_1 (); BDMFT_2 () we discuss first the application of the cavity method within the coherent states representation in the rest of this subsection. Inserting such a decomposition of the identity for each of the sites at each of the Suzuki-Trotter slices leads, in the continuous time limit , to the coherent state path integral expression of the partition function of the Bose-Hubbard model Negele ():

 Z = ∫N∏i=1DaiD¯aie−S , (19) S = ∫β0dτ⎡⎣N∑i=1(¯ai(τ)(∂τ−μ)ai(τ)+U2(¯ai(τ)ai(τ))2)−J∑⟨i,j⟩(¯ai(τ)aj(τ)+¯aj(τ)ai(τ))⎤⎦ . (20)

Here and in the following we use bold symbols to denote imaginary-time dependent quantities. and are two (formally conjugate) complex numbers indexing the coherent state at site and imaginary time , with the path-integral measure of this site. Following the same steps as in the study of the Ising model, the cavity method for a Bethe lattice of connectivity leads to the following self-consistency equation:

 ηcav(a,¯a)=1zcavw(a,¯a)∫c−1∏i=1DaiD¯aiηcav(ai,¯ai) exp{J∫β0dτ[¯a(τ)c−1∑i=1ai(τ)+a(τ)c−1∑i=1¯ai(τ)]} , (21)

with the on-site weight of the path given by

 w(a,¯a)=exp[−∫β0dτ(¯a(τ)(∂τ−μ)a(τ)+U2(¯a(τ)a(τ))2)] . (22)

This self-consistent equation on is formally similar to the corresponding Eq. (14) of the Ising model. It is however much more complicated: the Ising degree of freedom could only take two distinct values, whereas belongs to a space of functions of the imaginary time . Therefore is a functional measure whose representation is much more difficult; a complete parametrization requires the knowledge of all -points correlations . On the other hand, this is one of the most interesting features of the quantum cavity method: on-site quantum fluctuations are fully kept into account, without any truncation of higher order correlations.

Unfortunately, an exact solution of (21) can be easily obtained only in the case of free bosons (). acquires in this case a Gaussian form, with averages and two-point functions which can be computed exactly and reproduce the results obtained by direct diagonalization of the adjacency matrix of the Bethe lattice spectrum_Bethe ().

### iii.3 Large connectivity limit and the connection with B-DMFT

In the interacting case () a solution of (21) can be looked for in the limit of large connectivity, and this is precisely the road followed by the B-DMFT studies BDMFT_1 (); BDMFT_2 (). For completeness we shall explain in this subsection how to recover the B-DMFT formalism from Eq. (21), before turning in the next section to the occupation number basis which will allow to solve the model for any connectivity.

To state the large expansion let us rewrite (21) in a more convenient way, using the lighter Nambu notation , and consequently . Then we can rewrite (21) as

 (23)

is the generating functional of the connected correlation functions of and . It can be expanded as

 Γ(Φ)=∫β0dτ⟨Ψ†⟩Φ(τ)+∫β0dτdτ′Φ†(τ)ˆGc(τ−τ′)Φ(τ′)+O(Φ3) , (24)

where the averages are with respect to , we have used the cyclic invariance in imaginary time, and is the connected part of the two point correlator of .

In the large connectivity limit, the superfluid-insulator transition happens for a critical value of , with a finite . This can be argued by looking at Eq. (9), where it is clear that the dependence on hopping and connectivity is only through , which is the real control parameter. For large and , the cumulant expansion of thus becomes a systematic expansion in powers of . By keeping only a finite number of terms in the cumulant expansion, we obtain an expression of that is not Gaussian (because of the term in ); still, we can obtain closed equations for the cumulants by computing them self-consistently as averages over the non-Gaussian .

The leading order in gives, assuming without loss of generality that the average value of the order parameter is real, ,

 (25)

This last equation can be rewritten in the operator representation, which gives back the equation for corresponding to the minimization of the variational free-energy (9), up to a multiplicative constant in the definition of . Note that a generalization of the discussion above and of Eq. (25) to the disordered case leads to the stochastic mean field theory devised in SMFT ().

We will now show that the next-to-leading order in the cumulant expansion of (24) gives the B-DMFT equations recently derived in BDMFT_1 (); BDMFT_2 (). Note that the truncation at this two-point level was also used in the context of spin models in qcav_scardi (). Plugging the expansion (24) in (23), we obtain to order :

 ηcav(Ψ)=1zcavexp[−Sloc] ,Sloc=∫β0dτdτ′Ψ†(τ)ˆG−1(τ−τ′)Ψ(τ′)+∫β0dτ[U8(Ψ†(τ)Ψ(τ))2−Jc−1c⟨Ψ†⟩Ψ(τ)] ,ˆG−1(τ−τ′)=12(∂τ−μ00−∂τ−μ)δ(τ−τ′)−J2cˆGc(τ−τ′) . (26)

Then, and have to be computed self-consistently as averages with the local action . This set of equations correspond exactly to the B-DMFT of BDMFT_1 (); BDMFT_2 () for the special case of a Bethe lattice.

Away from these two limits ( and ) it seems difficult to obtain a solution of the cavity equation (21) as written in the coherent state basis. As a consequence we shall turn in the following to the representation number basis to apply the Suzuki-Trotter formula and thus obtain a more tractable equation for all values of and .

## Iv The quantum cavity method in the occupation number basis

### iv.1 The equations and the procedure for their numerical resolution

The insertion of a decomposition of the identity expressed in the occupation number basis in the Suzuki-Trotter formula (18) leads to an expression of the partition function of the Bose-Hubbard model as a sum over occupation number trajectories in imaginary time, . These trajectories are defined on an imaginary time interval of length , with the periodicity condition . The weight (action) of these trajectories has two origins: the local part of the Hamiltonian (1) yields a contribution of the form for each of the sites, where is the local energy term in the Hamiltonian. In addition the hopping term of the Hamiltonian imposes constraints between the occupation number trajectories: each time is raised (resp. decreased) by 1, the occupation number of one of the neighbors must decrease (resp. increase) of 1, meaning one particle has jumped from to (resp. from to ). Moreover each hopping event multiplies the weight of the trajectory by and by a coefficient depending on the instantaneous occupation numbers of the sites involved in the hopping. An explicit derivation of this representation shall be given in Sec. V. Here we wish to present the results of the computation in a lighter way for the ease of the reader not interested in the technical details (we follow the methodology developed for quantum spin models in qcav_our ()).

The recursion equations of the Bose-Hubbard model defined on a tree (or on a locally tree-like graph with the assumptions of the cavity method) can be expressed in terms of the probability distribution of the hopping trajectories defined on the edges of the graph. This quantity encodes the imaginary times at which a particle has crossed the edge, and the directions of the jumps; examples are given in Fig. 5. In the case of a regular Bethe lattice of connectivity one obtains the following self-consistent equation on ,

For simplicity we denote here with a sum symbol what is actually an integral over continuous degrees of freedom . The explicit expressions for the weights shall be given in Sec. V, see Eqs. (46), (47). Let us emphasize the formal similarity with the Ising equivalent Eq. (14), and of course the greater complexity of the basic degree of freedom in the Bose-Hubbard case, the hopping trajectory assuming values in a much larger space than the Ising variable . We already mentioned this problem while discussing the related equation (21) in the coherent state basis. Compared to this latter case we are however facing now an easier problem: the number of hopping events on a given edge is a (random) number which remains finite as long as the temperature is positive (it is actually of order ). A single hopping trajectory can thus be encoded as an integer , i.e. the number of particle jumps occuring on this edge during the imaginary time interval , reals precising the imaginary times where these jumps occur, and binary variables giving the direction of the jumps, see Fig. 5. In contrast the coherent state trajectories were those of two continuous functions with a priori no compact representation. This remark opens the way to an efficient numerical method for the resolution of (27111the code of our implementation is available upon request.. We can indeed follow the population dynamics strategy localization_Bethe (); MePa_Bethe (). The idea of this method is to represent numerically the probability distribution as a (weighted) sample of a large number of hopping trajectories, namely

 ηcav(y)=Ntraj∑i=1giδ(y−yi) , (28)

where the weights of the trajectories are normalized according to

 Ntraj∑i=1gi=1 . (29)

Sampling an element from the probability distribution corresponds in this representation to extract an integer with probability , and setting . This representation of is an approximation, which yields better and better numerical accuracy when grows. The determination of a sample of weights and trajectories which turns the representation (28) into a good approximation of the solution of (27) can be performed iteratively. To explain this point let us first rewrite the self-consistent equation (27) as

 ηcav(y)=∑y1,…,yc−1ηcav(y1)…ηcav(yc−1)P(y|y1,…,yc−1)Z(y1,…,yc−1)zcav , (30)

where we have defined

Constructed in this way is a conditional probability distribution over . Given the values of the hopping trajectories it is actually possible to perform an exact sampling from and to compute the normalization constant . In fact this reduces to a relatively simple single site problem: one has to construct the occupation trajectory of a site, and the associated hopping trajectory towards one of its neighbors, given the hopping trajectories on the other adjacent edges (see Fig. 5 and recall the pictorial representation given for the Ising case in Fig. 3). The hopping trajectories on the upper edges impose that changes by at the times particles arrive or depart from the considered central site. Between these times we shall show in Sec. V that the single site problem is described by an effective Hamiltonian ; each change in the value of provoked by the creation/annihilation operators in this effective Hamiltonian is associated to an hopping event in the new (downwards) trajectory . An example of this construction is given in Fig. 5. The computation of can be performed recognizing it as the partition function of the single-site effective Hamiltonian.

Suppose now a representation of is available under the form (28); one can plug it into the right-hand-side of Eq. (30) and construct a new set of trajectories and weights corresponding to its left-hand-side, let us call them . Independently and identically for each , this can be done by:

• drawing integers , each of them independently with probability

• drawing from

• setting ,

The new weights are then normalized to fulfill Eq. (29). This process can be iterated, plugging in the r.h.s. of Eq. (30) the newly obtained representation , and so on and so forth. After a certain number of steps, starting from an arbitrary initial condition222The initial condition should however allow for the possible spontaneous symmetry breaking ; we shall discuss this point further in Sec. V., a stationary solution is reached.

The physical observables can be computed from this representation of . Let us first explain the determination of the average occupation number of one site. Following our convention we denote the occupation number imaginary-time trajectory . Its probability distribution is expressed by considering the complete environment of a vertex with its neighbors, as was done for the local magnetization of the Ising model in Eq. (15) and schematized in Fig. 4,

 η(n)=1zsite∑y1,…,ycηcav(y1)…ηcav(yc)wsite(n,y1,…,yc) . (32)

The explicit expression of the weight shall be given in Sec. V. It is however intuitive that it will contain a factor arising from the local part of the Hamiltonian, and requires a consistency between and . In fact all the discontinuities in the occupation trajectory are fixed by the hopping events on the neighboring edges. If the number of hoppings towards the considered site equals the number of jumps outside, then the trajectory is fixed up to a global shift with independent of time. Whenever there is an unbalanced number of hopping towards/outside the vertex no periodic trajectory can be constructed, for such a configuration of . The average occupation number of a site is easily obtained from this probability ,

 ⟨a†a⟩=∑nη(n)n(0)=∑y1,…,ycηcav(y1)…ηcav(yc)∑nwsite(n,y1,…,yc)n(0)∑y1,…,ycηcav(y1)…ηcav(yc)∑nwsite(n,y1,…,yc) , (33)

where in the last expression we have explicited the normalization constant . Note that the time where is evaluated is arbitrary because of the cyclic invariance along the imaginary time axis, hence can be equivalently replaced by the temporal average . Given a representation of as a weighted sample (28) the numerical determination of is straightforward. Both the numerator and the denominator of (33) have the form of the average of a function with the ’s independently drawn from their distribution . As already explained drawing from a trajectory amounts to draw an index with probability and picking the element of the population. In formula this method of sampling leads to

 ∑y1,…,ycηcav(y1)…ηcav(yc)u(y1,…,yc)=1NtriesNtries∑j=1u(yij,1,…,yij,c) , (34)

where are integers drawn independently with the probability , with and . The numerical accuracy of the sampling is increased by taking a large value of . Moreover one can interleave average steps with iteration steps, determining a new set of weights and trajectories and recomputing independent averages on this new representation of .

We shall show in Sec. V that all the other observables defined in Sec. II (order parameter , free-energy, kinetic/potential energy and Green function) can also be expressed as averages of the form (34), and thus can be efficiently computed by the method developed here.

### iv.2 Results

In the following we present the results obtained for the Bose-Hubbard model on the Bethe lattice solved with the cavity method. We have shown in the previous sections that in the thermodynamic limit the exact resolution of the model amounts to solve Eq. (27) and we explained how this can be done with arbitrary numerical precision using the population dynamics algorithm described in subsection IV.1, encoding as a population of a large number of trajectories. We typically used and checked that the results were unchanged when using larger . We recall that in the thermodynamic limit the local observables (say, ) are independent on the site . In the following we focus on Bethe lattices of connectivity and , which mimic the connectivity of the two and three-dimensional square and cubic lattices, respectively.

Note that in order to implement the summations over involved in the computation of the various weights appearing in Eqs. (27) and (32), we have to introduce a cutoff on the possible values of , . The technical details will be discussed in section V. Since all the results presented below corresponds to regions of the parameter space such that , we work with a cut-off on the occupation number in the interval . We have checked for various values of the parameters that all the results are stable with respect to changes of the cutoff (in some cases also up to ).

In this section we also compare our results with other analytic approaches, such as the variational mean-field treatment described in Sec. II.2 and the bosonic Dynamical Mean Field Theory (B-DMFT)333 Note that we compare our results for a Bethe lattice with finite connectivity with those obtained in BDMFT_2 () using the semi-circular density of states. The latter is strictly valid in the limit of large connectivity, but this effect should be of subleading order in the large connectivity expansion. discussed in Sec. III.3, that correspond to the large limit of our method BDMFT_1 (); BDMFT_2 (). We also compare them to Quantum Monte Carlo simulations on square lattice bh_qmc_2d () and cubic lattice bh_qmc_3d ().

#### iv.2.1 Thermodynamic observables

Once the numerical iterative procedure has converged to a fixed point solution, one can compute the thermodynamic observables. Fig. 6 shows the behavior of and , as a function of the chemical potential, for , and for a Bethe lattice of connectivity . We have also computed and at lower values of the temperature (, , and ), but we find that the results are unchanged within our numerical accuracy, except very close to the lowest- superfluid-insulator transitions, where a little temperature effect is detected444 In particular the cusp in the kinetic energy at transition from the SF to the MI becomes sharper and becomes slightly smaller as we go from to in the first superfluid region. No further modifications are observed within our numerical errors for . between and (while no change is detectable above ). Hence, we can safely assume that for we are in the zero temperature regime.

The plot of Fig. 6 clearly shows a sequence of transitions from superfluid phases (SF) to Mott insulators (MI). The superfluid phase is characterized by a finite value of , which corresponds to a finite condensate fraction , while in the Mott insulator . In the limit of zero temperature in the MI regions the average number of bosons per site, is fixed to integer values, , , , . As a result, the compressibility of the system, vanishes identically. At small but finite temperature, we expect the compressibility to be exponentially small in the temperature, and indeed in Fig. 6 we observe that the compressibility at commensurate density is practically zero even at .

As shown in the inset of Fig. 6, the critical exponent for at the transition from the SF phase to the MI is consistent with the mean-field value . As in the classical case, the critical exponents found on the Bethe lattice are equal to the mean-field ones (see qcav_our () for a detailed discussion in the case of the quantum spin-1/2 ferromagnet on the Bethe lattice). Therefore, we also expect the other critical exponents to coincide with their mean-field values.

In Fig. 7 we plot the total energy, , as a function of the chemical potential for the same values of the parameters as in Fig. 6. We have also computed the free energy of the system, , but we do not show it since at these low temperatures the difference between the energy and the free energy is very small, and practically undetectable within our numerical precision, which confirms that we are effectively in the zero-temperature regime and the system can be considered to be in its ground state. In the inset of Fig. 7 we show the behavior of the kinetic energy, as a function of . The kinetic energy is, of course, lower in the SF phase than in the MI. Moreover we find that does not vanish in the MI phase, differently from the traditional mean-field approach, where hopping is completely suppressed in this phase. Note that both the kinetic energy and the potential energy show a discontinuity in the first derivative at the transition between the MI and the SF phase. On the contrary, the first derivative of the total energy (and of the free energy as well) is continuous at the transition; only its second derivative shows a discontinuity.

#### iv.2.2 Phase diagram

In Fig. 8 we present the “zero temperature” (recall that our method only works at finite temperature, hence we use quotes to remind that an extrapolation to is needed) phase diagram of the Bose-Hubbard model on the Bethe lattice of connectivity (left panel) and (right panel) in the plane. The phase boundaries have been obtained by computing as a function of at constant . In this case we used and (lower values of and higher values of have been used at higher ). As already discussed, at these values of we are effectively in the zero temperature regime. The lobelike shape of the phase boundary between MI and SF phases is qualitatively similar to the one found in the mean field approach. In Fig. 8 we also compare the results found with the quantum cavity method with the outcomes of Quantum Monte Carlo simulations (on the square lattice for  bh_qmc_2d () and on the cubic lattice for  bh_qmc_3d ()), and with the results of other analytical approaches. This comparison shows that the cavity method does a fairly good job in locating the lobelike contours between the MI and the SF phase and performs much better than the mean-field approach. It also performs slightly better than the Bosonic DMFT BDMFT_2 () for , although the difference between the cavity method and the B-DMFT is expected to become small as the connectivity is increased. It would be interesting in this respect to compare the cavity method and B-DMFT for .

In addition to the low temperature phase diagram we have computed the transition temperature line from the “insulating” phase (which we defined by , but has still a finite conductivity at finite temperature) to the superfluid one (defined by ) at unit filling (), for Bethe lattices of connectivity . This result is plotted in Fig. 9 and compared to the mean-field prediction and to the three-dimensional Monte-Carlo simulations of bh_qmc_3d ().

#### iv.2.3 Green functions and particle occupation statistics

As far as the results shown up to now are concerned, there are no qualitative differences between the Bethe lattice and the fully connected models (or equivalently the mean field approximation). We shall now present the results for the one-particle on-site imaginary time Green function, which exhibits a richer behavior with respect to the mean-field description. This quantity, defined in Eq. (6) as , is plotted in Fig. 10 for two set of values of the parameters , one in the MI phase, the other in the SF, both very close to the tip of the first insulating lobe.

Let us first remark that the decay of the Green function is independent on temperature in this temperature range, the effect of temperature being only to cut the decay around due to periodicity. Therefore we can safely assume that the Green functions reported in Fig. 10 are good estimates of their zero-temperature limit.

For the interpretation of these curves it is worth recalling the spectral representation of the Green function Negele () at zero temperature: for , </