# Strong-coupling solution of the bosonic dynamical mean-field theory

## Abstract

We derive an approximate analytical solution of the self-consistency equations of the bosonic dynamical mean-field theory (B-DMFT) in the strong-coupling limit. The approach is based on a linked-cluster expansion in the hybridization function of normal bosons around the atomic limit. The solution is used to compute the phase diagram of the bosonic Hubbard model for different lattices. We compare our results with numerical solutions of the B-DMFT equations and numerically exact methods, respectively. The very good agreement with those numerical results demonstrates that our approach captures the essential physics of correlated bosons both in the Mott insulator and in the superfluid phase. Close to the transition into the superfluid phase the momentum distribution function at zero momentum is found to be strongly enhanced already in the normal phase. The linked-cluster expansion also allows us to compute dynamical properties such as the spectral function of bosons. The evolution of the spectral function across the transition from the normal to the superfluid phase is seen to be characteristically different for the interaction driven and density driven transition, respectively.

## 1Introduction

Cold atoms in optical lattices provide a fascinating new class of interacting quantum many-particle systems.[1] Due to the unprecedented precision of experimental techniques in this field it is now possible to simulate and experimentally test theoretical models.[3] In particular, experiments with bosonic atoms have revived the theoretical interest in the properties of the bosonic Hubbard model.[8] This model describes the quantum mechanical competition between the kinetic energy of lattice bosons, which is responsible for their Bose-Einstein condensation, and the repulsive interaction, which favors localization of the particles. The phase diagram of the bosonic Hubbard model was first calculated by Fisher *et al.* [8] within a static mean-field theory derived from the atomic limit. With the formulation of the bosonic dynamical mean-field theory (B-DMFT) [11] a comprehensive investigation scheme for correlated lattice bosons in the thermodynamic limit has become available, which allows one to calculate also dynamical properties such as spectral functions of the interacting bosons. The B-DMFT is a thermodynamically consistent, non-perturbative many-body approach which is applicable for all values of the input parameters, e.g., the interaction, density, and temperature. It leads to a set of nonlinear equations which need to be solved self-consistently. An exact solution can be found only in special cases, e.g., for the Falicov-Kimball model.[11] In general, the self-consistent equations have to be solved numerically or by employing approximate analytical methods. The experience with the fermionic DMFT [13] shows that both numerically exact (but computationally expensive) methods *and* approximate analytical methods are important to gain insight into the solution of the complicated self-consistency equations. So far solutions of the B-DMFT equations had to be obtained fully numerically. Hu and Tong,[16] and Hubener, Snoek, and Hofstetter[17] employed exact diagonalization (ED), and Anders *et al.*[18] made use of continuous-time quantum Monte Carlo (CT-QMC) to solve the B-DMFT equations. Analytical or semi-analytical solutions of the B-DMFT equations did not exist up to now.

In this paper we present an analytical strong-coupling solution of the B-DMFT derived by a linked-cluster expansion (LCE)[20] around the atomic limit. The method is analogous to the fermionic strong-coupling solver developed by Dai, Haule, and Kotliar[22] for the fermionic DMFT. While in the fermionic case the strong-coupling expansion is unable to capture the low temperature Fermi liquid physics due to the existence of a characteristic low energy (Kondo) scale, there is no such limitation in the bosonic case. Our approach differs from previous strong-coupling expansions to the bosonic Hubbard model[23] since they performed the expansion in the hopping amplitude.

The paper is organized as follows. We first introduce the B-DMFT and its self-consistency equations. Then we formulate the linked-cluster expansion and thereby derive a strong-coupling approximation to the B-DMFT equations. This is then applied to the Bethe lattice and the cubic lattice, both with coordination number , and to the Bethe lattice with . The phase diagrams of the bosonic Hubbard, model calculated in this way are compared with those obtained from numerical solutions of the B-DMFT computed with ED [17] and CT-QMC,[18] respectively, from numerically exact evaluations on a Bethe lattice,[28] and from numerical results obtained by direct Monte Carlo simulations of the bosonic Hubbard model.[29] The momentum distribution functions and spectral functions of correlated lattice bosons in the normal and the Bose-Einstein condensed phase are also calculated. Finally we discuss possible extensions of the approach.

## 2Cumulant expansion in the bosonic dynamical mean-field theory

The B-DMFT is the bosonic counterpart to the well-established DMFT for lattice fermions described by the Hubbard model. Its derivation is described in detail in Ref. . Here we focus on a single species of bosons. The expansion presented below is easily generalized to the case of more than one type of boson.

The bosonic Hubbard model is given by the Hamiltonian

where and are creation and annihilation operators, respectively, for a boson at a lattice site , is the hopping between lattice sites and , is the local interaction, and is the number operator of the local occupation. In this paper we consider nearest-neighbor hopping, i.e., for the nearest-neighbor sites ,, and otherwise. In the following we set the Boltzmann constant and the lattice spacing equal to unity.

### 2.1Local action of the B-DMFT

In the B-DMFT the -dimensional lattice problem is replaced by an effective single-site (“impurity”) problem in which the local interaction remains unchanged, but the rest of the lattice is replaced by two dynamical mean fields (“baths”) corresponding to bosons in the normal state and in the Bose-Einstein condensate, respectively.[11] The time evolution of bosons on a particular site is represented by the local Green function

where we used the imaginary time, finite temperature formalism and Nambu notation with

and the Bose-Einstein condensate (BEC) is described by the local order parameter

The impurity problem is defined by the local action

where is the chemical potential, is a lattice dependent parameter, and

is the condensate wave function, i.e., a dynamical mean field. The dynamical mean field corresponding to bosons in the normal state is represented by the hybridization function

The dynamical mean fields , , and are determined by the self-consistency equations

and

Here the notation indicates that the thermodynamic average is performed on a lattice with a cavity, i.e., with one site removed. We note that in equilibrium is constant. For finite dimensional lattices is related to the local BEC order parameter by

The self-consistency loop is closed by introducing the self-energy in the Matsubara frequency representation through the -integrated Dyson equation

and using the lattice Hilbert transform

The latter equation links the local Green function to the self-energy for a specific lattice described by the non-interacting density of states . The momentum dependent lattice Green function is then given by

where is the dispersion relation of the non-interacting system and is the self-consistent solution of equations -.

For a Bethe lattice with infinite connectivity () [30] the self-consistency conditions reduce to the simple expressions and In general, e.g., for a cubic lattice, the self-consistency equations - need to be solved numerically.

### 2.2Cumulant expansion

In order to solve the impurity problem defined above we use the cumulant (linked-cluster) expansion in the dynamical mean fields and . The action is further divided into two parts

where

and

The partition function of the impurity problem is thus written as

where is the partition function for the system described by , and denotes the thermodynamic average with respect to the action .

Now the exponential function appearing in the average is expanded, leading to an infinite series

The series is then re-exponentiated with the help of cumulants (i.e., connected -particle Green functions) [32]

Here the superscript indicates that only the connected part of the averages with respect to is included. Now the partition function can be calculated to the desired order in .

The above approximation is in the spirit of other strong-coupling expansions [27] and becomes exact in the atomic limit (, ). However, it should be stressed that it is *not* an expansion in the hopping amplitude but rather in the dynamical mean fields and . The fact that these fields are obtained self-consistently implies that all orders of the hopping amplitude contribute.[11]

In the following we perform the cumulant expansion to second order in and in the partition function . Since the Green function is determined by the functional derivative

the diagonal element and off-diagonal element are then of first order in and :

and

Furthermore, the local BEC order parameter is given by

The thermodynamic averages are performed as , with . The trace is calculated over the eigenstates of , which are obtained by an exact diagonalization of the Hamiltonian matrix which is represented in the occupation number basis. Since the local Hilbert space of for the bosonic impurity problem is infinite dimensional, the diagonalization has to be performed numerically which, in principle, implies a further approximation. The Hilbert space has to be cut off in the occupation number of the impurity. The error introduced thereby can be controlled by performing calculations with different values of the cut-off and choosing the smallest cut-off value such that the results do not differ within the required accuracy.[33]

## 3Application of the linked-cluster expansion to various lattices

In the following we apply the results of the LCE to the Bethe lattice and the cubic lattice, both with coordination number , as well as to the Bethe lattice with infinite connectivity (). Our results for the Bethe lattice with coordination number can benchmarked by the exact numerical solution based on the cavity method.[28]

### 3.1Bethe lattice with coordination number

In Figure 1 we show the results obtained with the LCE for the interaction dependence of the Bose-Einstein condensation temperature , as well as for the phase diagram *vs.* at . We also compare them with results from other methods: the exact numerical evaluation (cavity method) by Semerjian, Tarzia, and Zamponi,[28] the B-DMFT solution with ED by Hubener, Snoek, and Hofstetter,[17] and the static mean-field solution of Fisher *et al.* [8] The static mean-field and the ED results were calculated at , whereas the results of the cavity method were obtained for . The phase transition line *vs.* only weakly depends on at such low temperatures as can be seen in the upper panel of Figure 1, where below the curve is practically vertical. For this reason we conclude that the phase diagram presented in the lower panel of Figure 1 is essentially the ground state phase diagram.

The results shown in the lower panel of Figure 1 demonstrate that the agreement between the two B-DMFT solutions is excellent. Namely, the blue circles (LCE, this work) are seen to lie practically on the red line (ED from Ref. ). Apparently the transition from the Mott-insulator to the superfluid is well described by the LCE approximation, which expands to first order in the dynamical mean field . This is different from the case of the fermionic DMFT where the low temperature physics of the Hubbard model close to the metal-insulator transition can not be described by the strong-coupling approximation.[22]

The value of the transition temperature obtained by the B-DMFT and the cavity method, respectively, is significantly lower than the results obtained by the static mean-field theory.[8] Since the B-DMFT captures local fluctuations exactly we conclude that they are responsible for the lowering of and the associated increase of the size of the Mott lobes.

For strong interactions the system is a Mott insulator for most values of the chemical potential . Upon lowering the interaction the system enters the superfluid phase with an order parameter . For the values of the chemical potential between the Mott lobes the superfluid phase persists up to very large values of . Since the LCE calculations were performed at a low but finite temperature (), there is no superfluid phase at below (not discernible in the figure).

### 3.2Cubic lattice

#### Phase diagram

The phase diagram of the Bose-Hubbard model for the cubic lattice obtained from the B-DMFT with the LCE and with CT-QMC, respectively, is presented in Figure 2. These results are compared with the lattice quantum Monte Carlo (QMC) results.[29] The LCE results are shown for two different temperatures ( and ). It is evident that the size of the Mott lobes decreases with decreasing temperature. Upon lowering the temperature the computation of the phase boundary using the B-DMFT with the LCE was found to become more elaborate. As already noted in Ref. for the CT-QMC solver, the convergence of the DMFT cycle close to the phase transition is very slow and the initial guess of and has to be carefully chosen.

Figure 2 shows that there is a small quantitative difference between the results obtained by different methods. It is unlikely that these differences can be explained by the different temperatures used in the computations (the lattice QMC calculations [29] were performed at , which is lower than the temperature used in the B-DMFT calculations). Indeed, at such low temperatures the temperature dependence of the phase diagram is very weak, as discussed earlier for the Bethe lattice. Nevertheless, the overall agreement between the results obtained from the three different methods is clearly very good. As in the case of the Bethe lattice the local dynamical fluctuations described by the B-DMFT lead to an increase of the size of the Mott lobes compared to the static mean-field solution.

#### Momentum distribution

The momentum distribution function of the normal phase, calculated at , is found to have an interesting behavior close to the transition to the superfluid phase. As shown in the upper panel and the inset of Figure 3 the distribution is strongly peaked at already in the normal phase. In the B-DMFT the momentum dependence of the momentum distribution is expressed only through the non-interacting dispersion relation (cf. Eq. ). Therefore, implicitly determines the momentum distribution. The plots in Figure 3 show and for different values of upon approaching the phase transition at constant density .

The peak in the momentum distribution in the normal phase close to the phase transition was noted previously by Kato *et al.* [34] within QMC solution. The lower panel in Figure 3 shows a comparison between obtained for the same parameters using different methods. As pointed out by Freericks *et al.* [27] the increase in the occupation at is an effect which is only partially described by a strong-coupling expansion in the hopping amplitude. The B-DMFT does capture this enhancement, and our LCE results are in a very good agreement with the lattice QMC data of Ref. .

#### Spectral functions

The B-DMFT approach with the LCE solver also allows one to investigate the behavior of the -integrated spectral function across the phase transition from the superfluid to the Mott phase (Figs. ? and ?). Since in our current implementation of the LCE the computations are performed on the imaginary time or imaginary frequency axes, spectral functions at real frequencies have to be calculated by analytic continuation.[35] The spectral functions presented in Figs ? and ? were obtained by analytic continuation with Padé approximants. Calculations of bosonic spectral functions were also done with the functional renormalization group [36] and in the variational cluster approach (VCA).[38]

Here we focus on three distinct cases: (i) the interaction driven phase transition at the tip of the Mott lobe, keeping the ratio constant; (ii) the interaction driven transition at the bottom of the lobe, also with constant; and (iii) the density driven transition at constant interaction . Due to the approximation introduced by the analytic continuation one can draw only qualitative conclusions about the spectral density in the region close to the chemical potential (e.g., one cannot reliably estimate the size of the gap). Nevertheless the qualitative behavior and the spectral weight transfer is well illustrated and the difference between the three cases considered here is clearly visible. At the tip of the Mott lobe (case (i), left panel of Fig. ?) an increase of the interaction leads to a symmetric shift of the spectral weight on both sides of the chemical potential. At the same time a Mott gap opens and two Hubbard bands are formed (see the bottom plot in the left panel of Fig. ?). The shape of the bands vaguely resembles the non-interacting density of states for the cubic lattice. Away from the tip (case (ii), right panel of Fig. ?) the shift of the spectral weight is not symmetric with respect to the chemical potential. The lower Hubbard band resides close to the chemical potential, whereas the upper Hubbard band is shifted to higher frequencies. A different behavior is observed in the density driven transition (case (iii), Fig. ?). Upon increasing the chemical potential at constant interaction, the spectral function is shifted as a whole to lower frequencies, simultaneously forming a gap.

### 3.3The Bethe lattice with

The phase diagram for the Bethe lattice is presented in Figure 4. At sufficiently high temperatures (e.g., as in Figure 4) the LCE gives convergent results both for the superfluid and normal phases near the phase transition. However, at temperatures below we have not been able to find a convergent solution in the superfluid phase around the tip of the second Mott lobe. The iterations converge either to (normal phase), or to a solution with but with a non-concave, and hence unphysical,[40] . The results for are shown in Figure 4, where the solution at around the first Mott lobe converges both in the normal and the superfluid phase, thus making it possible to calculate the phase boundary. In the range a convergent solution was only obtained in the normal phase (), i.e., it was not possible to determine the phase boundary of the second lobe completely. As the temperature is lowered, the range of the chemical potentials for which we did not obtain a superfluid solution increases. For temperatures below we did not even obtain solutions with non-zero superfluid order parameter around the tip of the first Mott lobe. Upon further lowering the temperature, the region of convergence of the method in the superfluid phase is reduced to the values of near the edges of the lobes. At this moment it is not clear whether the absence of a solution in the superfluid phase in the Bethe lattice for some chemical potentials at low temperatures is a consequence of the strong-coupling approximation to the B-DMFT, or the B-DMFT itself. This is an open question which needs to be answered in the future. Such problems did not occur for the other lattices investigated here.

## 4Summary

We developed an analytical approximation scheme to solve the B-DMFT equations for correlated lattice bosons in the strong-coupling limit. The solution makes use of a linked-cluster expansion in the hybridization function of normal bosons around the atomic limit. Explicit results were obtained for the Bose-Hubbard model on the cubic lattice and the Bethe lattice with connectivity and , respectively. Remarkably good agreement with numerical solutions of the B-DMFT equations obtained with exact diagonalization [17], continuous-time quantum Monte Carlo [19], and direct lattice QMC calculations [29] was found. This agreement demonstrates that the strong-coupling solution derived here provides a correct description of the physics of correlated bosons. The method is computationally inexpensive and, with a good choice of the initial guess of the parameters, usually leads to a fast convergence of the iteration of the self-consistency equations. The Bethe lattice with infinite connectivity is an exception which still requires further investigation.

We also employed the linked-cluster expansion to calculate the momentum distribution function of normal bosons close to the phase transition as well as the bosonic spectral function in the normal and superfluid phase.

The approximation scheme presented in this paper can, in principle, be systematically improved by the inclusion of higher order terms. However, the non-interacting limit can only be reached if terms up to infinite order are included, e.g., by an appropriate resummation. This has been achieved for fermions by the non-crossing approximation (NCA).[41] The fundamental problem of the NCA, namely its failure to describe the low temperature Fermi liquid regime adequately owing to the existence of a characteristic coherence scale (the Kondo temperature), may be absent in the case of bosons where such a coherence scale does not exist. For that reason it should be clarified whether it is possible to construct a renormalized expansion for correlated bosons which is applicable for all temperatures and interaction strengths.

### References

- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemStop
- bibitemNoStop
- bibitemStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemNoStop
- bibitemStop
- bibitemNoStop