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

###### pacs:

71.10.Fd, 67.85.Hj## I Introduction

Cold atoms in optical lattices provide a fascinating new class of interacting quantum many-particle systems.Anglin and Ketterle (2002); Bloch et al. (2008) Due to the unprecedented precision of experimental techniques in this field it is now possible to simulate and experimentally test theoretical models.Greiner et al. (2002); Bloch (2005); Lewenstein et al. (2007); Jördens et al. (2008); Giorgini et al. (2008); Bloch et al. (2008) In particular, experiments with bosonic atoms have revived the theoretical interest in the properties of the bosonic Hubbard model.Fisher et al. (1989); Jaksch et al. (1998); ref () 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. Fisher et al. (1989) within a static mean-field theory derived from the atomic limit. With the formulation of the bosonic dynamical mean-field theory (B-DMFT) Byczuk and Vollhardt (2008, 2009) 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.Byczuk and Vollhardt (2008, 2009) In general, the self-consistent equations have to be solved numerically or by employing approximate analytical methods. The experience with the fermionic DMFT Metzner and Vollhardt (1989); Kotliar and Vollhardt (2004); Georges et al. (1996) 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,Hu and Tong (2009) and Hubener, Snoek, and HofstetterHubener et al. (2009) employed exact diagonalization (ED), and Anders et al.Anders et al. (2010, 2011) 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)Metzner (1991); Bartkowiak and Chao (1993) around the atomic limit. The method is analogous to the fermionic strong-coupling solver developed by Dai, Haule, and KotliarDai et al. (2005) 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 modelKampf and Zimanyi (1993); Sengupta and Dupuis (2005); dos Santos and Pelster (2009); Bradlyn et al. (2009); Freericks et al. (2009) 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 Hubener et al. (2009) and CT-QMC,Anders et al. (2010, 2011) respectively, from numerically exact evaluations on a Bethe lattice,Semerjian et al. (2009) and from numerical results obtained by direct Monte Carlo simulations of the bosonic Hubbard model.Capogrosso-Sansone et al. (2007) 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.

## Ii Cumulant 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. Byczuk and Vollhardt, 2008. 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

(1) |

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.

### ii.1 Local action of the B-DMFT

In the B-DMFT the -dimensional lattice problem (1) 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.Byczuk and Vollhardt (2008) The time evolution of bosons on a particular site is represented by the local Green function

(2) |

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

(3) |

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

(4) |

The impurity problem is defined by the local action

(5) |

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

(6) |

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

(7) |

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

(8) | |||||

and

(9) |

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 (4) by

(10) |

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

(11) |

and using the lattice Hilbert transform

(12) |

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

(13) |

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

### ii.2 Cumulant 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 (5) is further divided into two parts

(14) |

where

(15) | |||||

and

(16) |

The partition function of the impurity problem is thus written as

(17) |

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

(18) | |||||

The series is then re-exponentiated with the help of cumulants (i.e., connected -particle Green functions) Kubo (1962); Metzner (1991)

(19) | |||||

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

The above approximation is in the spirit of other strong-coupling expansions Freericks et al. (2009); Bradlyn et al. (2009) 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.Byczuk and Vollhardt (2008); Metzner (1991)

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

(20) |

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

(21) |

and

(22) |

Furthermore, the local BEC order parameter is given by

(23) |

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.cm2 ()

## Iii Application 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.Semerjian et al. (2009)

### iii.1 Bethe lattice with coordination number

In Fig. 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,Semerjian et al. (2009) the B-DMFT solution with ED by Hubener, Snoek, and Hofstetter,Hubener et al. (2009) and the static mean-field solution of Fisher et al. Fisher et al. (1989) 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 Fig. 1, where below the curve is practically vertical. For this reason we conclude that the phase diagram presented in the lower panel of Fig. 1 is essentially the ground state phase diagram.

The results shown in the lower panel of Fig. 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. Hubener et al., 2009). 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.Dai et al. (2005)

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.Fisher et al. (1989) 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).

### iii.2 Cubic lattice

#### iii.2.1 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 Fig. 2. These results are compared with the lattice quantum Monte Carlo (QMC) results.Capogrosso-Sansone et al. (2007) 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. Anders et al., 2011 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.

Fig. 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 Capogrosso-Sansone et al. (2007) 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.

#### iii.2.2 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 Fig. 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. (13)). Therefore, implicitly determines the momentum distribution. The plots in Fig. 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. Kato et al. (2008) within QMC solution. The lower panel in Fig. 3 shows a comparison between obtained for the same parameters using different methods. As pointed out by Freericks et al. Freericks et al. (2009) 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. Freericks et al., 2009.

#### iii.2.3 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. 4 and 5). 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.cm1 () The spectral functions presented in Figs 4 and 5 were obtained by analytic continuation with Padé approximants. Calculations of bosonic spectral functions were also done with the functional renormalization group Sinner et al. (2009, 2010) and in the variational cluster approach (VCA).Knap et al. (2010, 2011)

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. 4) 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. 4). 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. 4) 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. 5). Upon increasing the chemical potential at constant interaction, the spectral function is shifted as a whole to lower frequencies, simultaneously forming a gap.

### iii.3 The Bethe lattice with

The phase diagram for the Bethe lattice is presented in Fig. 6. At sufficiently high temperatures (e.g., as in Fig. 6) 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,cm3 () . The results for are shown in Fig. 6, 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.

## Iv Summary

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 Hubener et al. (2009), continuous-time quantum Monte Carlo Anders et al. (2011), and direct lattice QMC calculations Capogrosso-Sansone et al. (2007) 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).Costi et al. (1996) 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.

###### Acknowledgements.

This research was done when Anna Kauch worked at the University of Augsburg. Partial support by the Deutsche Forschungsgemeinschaft through TRR 80 is gratefully acknowledged. Krzysztof Byczuk acknowledges support by the grant No. N N202 103138 of the Polish Ministry of Science and Education.## References

- Anglin and Ketterle (2002) R. J. Anglin and W. Ketterle, Nature 416, 211 (2002).
- Bloch et al. (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008).
- Greiner et al. (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, Nature 415, 39 (2002).
- Bloch (2005) I. Bloch, Nature Phys. 1, 23 (2005).
- Lewenstein et al. (2007) M. Lewenstein, A. Sanpera, V. Ahufinger, B. Damski, A. Sen, and U. Sen, Adv. Phys. 56, 243 (2007).
- Jördens et al. (2008) R. Jördens, N. Strohmaier, K. Günter, H. Moritz, and T. Esslinger, Nature 455, 204 (2008).
- Giorgini et al. (2008) S. Giorgini, L. P. Pitaevskii, and S. Stringari, Rev. Mod. Phys. 80, 1215 (2008).
- Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
- Jaksch et al. (1998) D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, Phys. Rev. Lett. 81, 3108 (1998).
- (10) For a detailed list of references concerning the investigation of the bosonic Hubbard model in condensed matter theory prior to the research on cold atoms in optical lattices, see Ref. Byczuk and Vollhardt, 2008.
- Byczuk and Vollhardt (2008) K. Byczuk and D. Vollhardt, Phys. Rev. B 77, 235106 (2008).
- Byczuk and Vollhardt (2009) K. Byczuk and D. Vollhardt, Ann. Phys. 18, 622 (2009).
- Metzner and Vollhardt (1989) W. Metzner and D. Vollhardt, Phys. Rev. Lett. 62, 324 (1989).
- Kotliar and Vollhardt (2004) G. Kotliar and D. Vollhardt, Physics Today 3, 53 (2004).
- Georges et al. (1996) A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Rev. Mod. Phys. 68, 13 (1996).
- Hu and Tong (2009) W.-J. Hu and N.-H. Tong, Phys. Rev. B 80, 245110 (2009).
- Hubener et al. (2009) A. Hubener, M. Snoek, and W. Hofstetter, Phys. Rev. B 80, 245109 (2009).
- Anders et al. (2010) P. Anders, E. Gull, L. Pollet, M. Troyer, and P. Werner, Phys. Rev. Lett. 105, 096402 (2010).
- Anders et al. (2011) P. Anders, E. Gull, L. Pollet, M. Troyer, and P. Werner, New J. Phys. 13, 075013 (2011).
- Metzner (1991) W. Metzner, Phys. Rev. B 43, 8549 (1991).
- Bartkowiak and Chao (1993) M. Bartkowiak and K. A. Chao, Phys. Rev. B 47, 1616 (1993).
- Dai et al. (2005) X. Dai, K. Haule, and G. Kotliar, Phys. Rev. B 72, 045111 (2005).
- Kampf and Zimanyi (1993) A. P. Kampf and G. T. Zimanyi, Phys. Rev. B 47, 279 (1993).
- Sengupta and Dupuis (2005) K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 (2005).
- dos Santos and Pelster (2009) F. E. A. dos Santos and A. Pelster, Phys. Rev. A 79, 013614 (2009).
- Bradlyn et al. (2009) B. Bradlyn, F. E. A. dos Santos, and A. Pelster, Phys. Rev. A 79, 013615 (2009).
- Freericks et al. (2009) J. K. Freericks, H. R. Krishnamurthy, Y. Kato, N. Kawashima, and N. Trivedi, Phys. Rev. A 79, 053631 (2009).
- Semerjian et al. (2009) G. Semerjian, M. Tarzia, and F. Zamponi, Phys. Rev. B 80, 014524 (2009).
- Capogrosso-Sansone et al. (2007) B. Capogrosso-Sansone, N. Prokof’ev, and B. Svistunov, Phys. Rev. B 75, 134302 (2007).
- Weaire and Thorpe (1971) D. Weaire and M. F. Thorpe, Phys. Rev. B 4, 2508 (1971).
- Eckstein et al. (2005) M. Eckstein, M. Kollar, K. Byczuk, and D. Vollhardt, Phys. Rev. B 71, 235119 (2005).
- Kubo (1962) R. Kubo, J. Phys. Soc. Jap. 17, 1100 (1962).
- (33) Results presented here were usually computed with the maximal occupation number equal to 5, which introduces an error much below the resolution of the rest of the calculations.
- Kato et al. (2008) Y. Kato, Q. Zhou, N. Kawashima, and N. Trivedi, Nature Phys. 4, 617 (2008).
- (35) Our experience is that the convergence is much faster in the imaginary time implementation, so we have dropped the real frequency approach for the moment.
- Sinner et al. (2009) A. Sinner, N. Hasselmann, and P. Kopietz, Phys. Rev. Lett. 102, 120601 (2009).
- Sinner et al. (2010) A. Sinner, N. Hasselmann, and P. Kopietz, Phys. Rev. A 82, 063632 (2010).
- Knap et al. (2010) M. Knap, E. Arrigoni, and W. von der Linden, Phys. Rev. B 81, 024301 (2010).
- Knap et al. (2011) M. Knap, E. Arrigoni, and W. von der Linden, Phys. Rev. B 83, 134507 (2011).
- (40) It follows from the definition of the diagonal part of the bosonic Green function in imaginary time that it has to be concave.
- Costi et al. (1996) T. A. Costi, J. Kroha, and P. Wölfle, Phys. Rev. B 53, 1850 (1996).