Exact solution of the BoseHubbard model on the Bethe lattice
Abstract
The exact solution of a quantum Bethe lattice model in the thermodynamic limit amounts to solve a functional selfconsistent equation. In this paper we obtain this equation for the BoseHubbard 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 BoseEinstein 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 BoseHubbard 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 condensedmatter physics.
From a more theoretical point of view the BoseHubbard model has been studied with various techniques, notably meanfield 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 meanfield 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 meanfield level. As often the assumptions underlying the meanfield 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 MeanField 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 meanfield treatment of FiFi (), by treating the BoseHubbard 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 selfconsistent functional equation. By writing this equation in the coherent state basis one recovers the BDMFT 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 meanfield category, yet their finite connectivity makes them richer than the (fullyconnected) meanfield 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 nontrivial for the Bethe lattice version of the model, at variance with the meanfield 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 meanfield 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 meanfield 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 BDMFT 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 BoseHubbard 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 BDMFT 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 meanfield approach
ii.1 Definition of the model and of the observables
We shall consider the BoseHubbard model for spinless bosons on a graph of vertices (sites) defined by the following Hamiltonian,
(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 onsite 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 multipoint imaginary time correlations), in the following we will be mainly interested in standard observables such as the mean density, kinetic and onsite energy, condensate fraction, and Green functions. For the sake of clarity we recall now their definitions. The partition function at temperature is defined by
(2) 
where (we set ), the corresponding free energy is ; the freeenergy per site is . The thermodynamic average of an operator is defined as
(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
(4) 
where denotes the thermodynamic average in presence of the perturbation.
The imaginarytime evolution of an operator is given by ; we then define the twotimes correlation functions as Negele ()
(5) 
The Green function is defined, for , by
(6) 
where is the usual timeordering operation, which should not be confused with the temperature. Note that the cyclic property of the trace implies that , hence
(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 oneparticle density matrix .
ii.2 Review of meanfield 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 meanfield treatments FiFi (); Gutzwiller () of the problem. Let us briefly recall the various points of view and methodology that the meanfield approach usually encompasses. It can first be seen as an approximation to finitedimensional 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,
(8) 
the being here free parameters. Note that breaks the particle conservation symmetry. The freeenergy at inverse temperature of the original model, , can be upperbounded as , where is the freeenergy 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 freeenergy per site can then be upperbounded as , with
(9) 
where the trace is over the Hilbert space of a single site. The physical properties of this variational freeenergy are wellknown 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 freeenergy per site of the BoseHubbard Hamiltonian. The usual meanfield 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 sitedecoupled 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 selfconsistency equations leading finally to the same expression (9) of the freeenergy 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 freeenergy of the system. Finally, another reasoning yielding this meanfield result consists in devising a model which has exactly (9) as its freeenergy per site in the thermodynamic limit, instead of taking it as an approximation for the finitedimensional case. As could be expected this model corresponds to the fullyconnected 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 freeenergy of this fullyconnected model converges indeed to (9).
The above described meanfield 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 BoseHubbard model, for which it forbids the description of the Bose Glass phase predicted in FiFi (). In order to cure this defect of the meanfield treatment one has to account in some way for the correlations between neighboring sites. Thinking in terms of classical spin models, the meanfield 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 BoseHubbard 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 meanfield approximation was exact for the fullyconnected 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 treelike, 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 finitedimensional structure gives them a meanfield character, but their finite connectivity leads to a richer content than the fullyconnected models. The freeenergy of Bethe lattice models is selfaveraging 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 socalled 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,
(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 GibbsBoltzmann probability measure is
(11) 
with the partition function ensuring its normalization. The goal is to compute the freeenergy 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,
(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
(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 GibbsBoltzmann 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 treelike (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 selfconsistent boundary conditions and avoids the illdefined 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 freeenergy 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 selfconsistent equation
(14) 
with a normalization constant. A pictorial representation of this equation can be found in Fig. 3. The local magnetization is then computed as
(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,
(16) 
solutions of
(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 spinspin 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 fluctuationdissipation 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 SuzukiTrotter formalism and quantum cavity method
We have just seen how the cavity method deals with classical spin models on locally treelike graphs. The extension of the method to quantum models is based on the SuzukiTrotter 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 BoseHubbard model this decomposition splits the Hamiltonian between the hopping term and the onsite potential energy . The noncommutativity of can be cured with the application of the SuzukiTrotter formula suzuki (),
(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 timedependent classical degrees of freedom. In particular, if the interactions of the quantum model are defined on a treelike 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 SuzukiTrotter 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 SuzukiTrotter 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 BDMFT 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 SuzukiTrotter slices leads, in the continuous time limit , to the coherent state path integral expression of the partition function of the BoseHubbard model Negele ():
(19)  
(20) 
Here and in the following we use bold symbols to denote imaginarytime dependent quantities. and are two (formally conjugate) complex numbers indexing the coherent state at site and imaginary time , with the pathintegral 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 selfconsistency equation:
(21) 
with the onsite weight of the path given by
(22) 
This selfconsistent 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: onsite 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 twopoint 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 BDMFT
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 BDMFT studies BDMFT_1 (); BDMFT_2 (). For completeness we shall explain in this subsection how to recover the BDMFT 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
(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 superfluidinsulator 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 selfconsistently as averages over the nonGaussian .
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 freeenergy (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 nexttoleading order in the cumulant expansion of (24) gives the BDMFT equations recently derived in BDMFT_1 (); BDMFT_2 (). Note that the truncation at this twopoint level was also used in the context of spin models in qcav_scardi (). Plugging the expansion (24) in (23), we obtain to order :
(26) 
Then, and have to be computed selfconsistently as averages with the local action . This set of equations correspond exactly to the BDMFT 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 SuzukiTrotter 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 SuzukiTrotter formula (18) leads to an expression of the partition function of the BoseHubbard 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 BoseHubbard model defined on a tree (or on a locally treelike 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 selfconsistent equation on ,
(27) 
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 BoseHubbard 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 (27) ^{1}^{1}1the 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
(28) 
where the weights of the trajectories are normalized according to
(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 selfconsistent equation (27) as
(30) 
where we have defined
(31) 
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 singlesite effective Hamiltonian.
Suppose now a representation of is available under the form (28); one can plug it into the righthandside of Eq. (30) and construct a new set of trajectories and weights corresponding to its lefthandside, 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 condition^{2}^{2}2The 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 imaginarytime 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,
(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 ,
(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
(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 .
iv.2 Results
In the following we present the results obtained for the BoseHubbard 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 threedimensional 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 cutoff 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 meanfield treatment described in Sec. II.2 and the bosonic Dynamical Mean Field Theory (BDMFT)^{3}^{3}3 Note that we compare our results for a Bethe lattice with finite connectivity with those obtained in BDMFT_2 () using the semicircular 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 superfluidinsulator transitions, where a little temperature effect is detected^{4}^{4}4 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 meanfield value . As in the classical case, the critical exponents found on the Bethe lattice are equal to the meanfield ones (see qcav_our () for a detailed discussion in the case of the quantum spin1/2 ferromagnet on the Bethe lattice). Therefore, we also expect the other critical exponents to coincide with their meanfield 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 zerotemperature 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 meanfield 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 BoseHubbard 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 meanfield approach. It also performs slightly better than the Bosonic DMFT BDMFT_2 () for , although the difference between the cavity method and the BDMFT is expected to become small as the connectivity is increased. It would be interesting in this respect to compare the cavity method and BDMFT 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 meanfield prediction and to the threedimensional MonteCarlo 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 oneparticle onsite imaginary time Green function, which exhibits a richer behavior with respect to the meanfield 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 zerotemperature limit.
For the interpretation of these curves it is worth recalling the spectral representation of the Green function Negele () at zero temperature: for , </