Dynamical mean-field theory for bosons

We discuss the recently developed bosonic dynamical mean-field (B-DMFT) framework, which maps a bosonic lattice model onto the selfconsistent solution of a bosonic impurity model with coupling to a reservoir of normal and condensed bosons. The effective impurity action is derived in several ways: (i) as an approximation to the kinetic energy functional of the lattice problem, (ii) using a cavity approach, and (iii) by using an effective medium approach based on adding a one-loop correction to the selfconsistently defined condensate. To solve the impurity problem, we use a continuous-time Monte Carlo algorithm based on a sampling of a perturbation expansion in the hybridization functions and the condensate wave function. As applications of the formalism we present finite temperature B-DMFT phase diagrams for the bosonic Hubbard model on a 3d cubic and 2d square lattice, the condensate order parameter as a function of chemical potential, critical exponents for the condensate, the approach to the weakly interacting Bose gas regime for weak repulsions, and the kinetic energy as a function of temperature.

Particles with bosonic statistics can macroscopically occupy a single mode at low enough temperature, even in the absence of correlation. This phenomenon is known as Bose-Einstein condensation, and leads to a variety of phases in strongly correlated bosonic systems. A typical example is 4 He, which exhibits normal, superfluid, solid, and possibly supersolid phases. Another example are dilute ultracold atomic gases, which are well described by a Bogoliubov Hamiltonian. Strong interaction effects can be induced by adding lasers producing an optical lattice. The resulting system is an essentially clean realization of the Bose-Hubbard model, which describes the competition between hopping and on-site repulsion. Commensurability effects in the lattice lead to a phase transition from a superfluid to a Mott insulator at integer fillings and strong enough interaction. 1 Both cold gases and 4 He can be controlled experimentally with great accuracy and are virtually free of impurities and disorder. Cold atomic gases have the additional flexibility of tunable interaction strengths, and provide the freedom of changing the mass by choosing different alkali atoms.
For both 4 He and ultracold bosonic gases in an optical lattice, powerful numerical methods exist for the strongly-interacting regime. Path integral simulations based on the worm algorithm can sample up to 10, 000 bosonic atoms at T = 1K for supersolid 4 He, enabling the study (see Ref. 2 for a review) of individual defects such as vacancies, 3 dislocations, 4 and grain boundaries. 5 For cold gases, up to 1.5 × 10 6 atoms were studied 6 and compared directly to experiment with excellent agreement in time-of-flight images. 7 Such simulations involve a stochastic evaluation of all connected and disconnected diagrams occurring in a high-temperature series expansion on a finite lattice, which is an expansion in the hopping or kinetic energy (over temperature) around the atomic limit. [8][9][10] The Monte Carlo simulation of fermionic lattice models suffers from the notorious sign problem, which prevents the study of large systems in the most interesting parameter regime. A computationally tractable approximate method to simulate these models is dynamical mean-field theory (DMFT). [11][12][13][14][15] In these calculations, the manybody selfenergy is approximated by all local skeleton diagrams involving local propagators only, which implies a selfconsistent determination of the selfenergy and the local propagators. Non-local contributions are neglected. This simplification is convenient because the approximated selfenergy can be evaluated efficiently from an appropriately defined impurity action. 12 By using sign-free (for single-site DMFT) efficient continuous-time Monte Carlo solvers, 16 one obtains the full Green's function as a solution to the effective impurity action in polynomial time. [17][18][19] The simplification of the diagrammatic structure 11 allows to define DMFT for arbitrary dimensions and lattice structures. A major success of DMFT is the understanding of the Mott metal-insulator transition it has provided. 13,15 DMFT has been extensively used to study model systems and -in conjunction with band structure techniques -to compute material properties for a wide range of compounds. 15,20 Several extensions make the approximation systematic and controlled: Cluster methods 14 like the dynamical cluster approximation 21 or the cellular DMFT 22,23 reintroduce momentum dependence by considering multi-site impurity clusters. Methods like DΓA 24 or dual fermions 25,26 systematically consider non-local diagrams beyond DMFT. In principle, a diagrammatic Monte Carlo evaluation of all neglected contributions to the selfenergy allows to estimate its accuracy, as, e.g., recently done for the Anderson localization problem. 27 The formulation of an analogous dynamical mean-field theory for bosonic lattice models has proven difficult. On the one hand, from the perturbative diagrammatic point of view, the idea of retaining an infinite subclass of bare diagrams may seem dubious since standard Feynman diagrammatic expansions in U fail notoriously for the Bose-Hubbard model due to Dyson's collapse argument: In the complex plane, the convergence radius is zero as all bosons would collapse to a single point for negative interactions, with an infinite negative energy in the thermodynamic limit. Hence, it seems impossible to think of meaningful diagrammatic expansions in U for bosons, and only a few techniques are known to deal with this problem such as Kleinert's variational perturbation theory, 28 the cutting off of large field contributions, 29 or the use of a sequence with appropriately chosen counter terms rendering an infinite convergence radius in the absence of phase transitions. 30 In the Baym-Kadanoff approach to DMFT, a subset of diagrams consisting of all contributions from local dressed propagators are retained. These skeleton diagrams may have radically different convergence properties than the bare series but those properties are essentially unknown. In the case of absolutely convergent series the bare and skeleton series are equivalent and physically meaningful, but due to Dyson's collapse argument for the bare series there is no guarantee that the skeleton approach will converge.
On the other hand, keeping knowledge of the series provenance and by using the effective action for a single site (ie., restoring the action on the basis of the series), may still be worthwhile. In the normal phase the usual DMFT formalism can be applied in the same way as it is usually done for fermions. By using an effective action on a single site (which can be solved with the method of choice), most of the aforementioned convergence issues can be sidestepped. It is in fact only the occurrence of a condensate (which happens in the single particle channel) in the broken symmetry phase that poses a challenge in the development of a B-DMFT formalism (cf. Ref 31). We will define our B-DMFT theory at the one-loop level beyond this selfconsistently defined mean-field (condensate). As shown in App. B, our effective action for the impurity problem is the most general action for an impurity with a broken U (1) symmetry and a one-loop correction. 32 The hybridizations are then determined selfconsistently with the underlying lattice problem. The approximations involved can still be understood in the language of diagrams, but a full interpretation in terms of Baym-Kadanoff functionals remains subjective in light of the asymptotic series expansions. We therefore prefer to use alternative derivation schemes in this work that do not rely on an expansion in the bosonic repulsion U .
In the limit of infinite coordination number the decoupling approximation for the Bose-Hubbard model becomes exact (see the Appendix of Ref. 1). This decoupling approximation is recovered in our action at the mean-field level, since the one-loop correction vanishes. The original B-DMFT paper by Byczuk and Vollhardt 33 (as well as subsequent work presented by Hu and Tong 34 ) is based on the assumption that in this limit the kinetic and potential energy contributions in the broken symmetry phase can still be comparable to each other, in apparent contrast to Ref. 1. Their work postulates a scaling Ansatz with different scalings for condensed and non-condensed bosons, leading to an effective action that is different from the one we shall describe. The authors of Ref. 35 performed an expansion in the inverse of the coordination number around the limit of infinite coordination number (i.e. the static mean-field result of Ref. 1), but treated the condensate in a perturbative way only valid on the Bethe lattice. After the publication of our results, 36 they corrected for this and also derived a fully selfconsistent (and general) version of the B-DMFT formalism. 37 The virtue of extending the DMFT framework to interacting Bose systems lies in the fact that certain model systems otherwise not amenable to bosonic simulations, e.g. models with frustrated interactions, or Bose-Fermi mixtures can now be studied numerically. The quality of the approximation Σ(k, ω) → Σ skel (ω)[G loc ] is system dependent and needs to be established on a case-by-case basis. We will see that the B-DMFT approximation is very good for the single-component Bose-Hubbard model in 3d, and thus may serve as a good starting point for more complicated systems.
In both the fermionic and bosonic versions of dynamical mean-field theory, the computationally challenging part is the solution of the impurity problem. For fermionic systems exact diagonalization, 38 semi-analytical resummation of diagrams, 39,40 quantum Monte Carlo, 41 and numerical renormalization group 42 methods are in wide use for single-orbital models. In recent years, significant progress has been made with the development of diagrammatic Monte Carlo impurity solver techniques, based on an expansion of the partition function in powers of the interaction 17,18 or the impuritybath hybridization, 19,43 allowing access to much larger impurity clusters, lower temperatures, and more general interactions. 16 In this paper we present a detailed derivation of the B-DMFT equations and show how the impurity-condensate coupling must be chosen to obtain a consistent theory. We quantify the errors introduced by the dynamical mean-field approximation for a system with finite coordination number by comparing with lattice Monte Carlo methods, and we describe the impurity solver proposed in Ref. 36 in such detail that the implementation of the method becomes straightforward.
The paper is organized as follows: Sec. II introduces the Bose Hubbard model, and in Sec. III we derive the B-DMFT formalism, which is summarized in Sec. IV. In section V we describe the diagrammatic Monte Carlo impurity solver. Section VI discusses solvable limits while B-DMFT results for interacting bosons on a Bethe and 3d simple cubic lattice are presented in Sec. VII. Sec. VIII provides a summary and outlook.

II. THE BOSE-HUBBARD MODEL
We consider a model of spinless bosons on a lattice, described by the Bose-Hubbard Hamiltonian in standard lattice notation, where t denotes the hopping amplitude, U the on-site interaction and µ the chemical potential. Unless otherwise written, our unit will be the hopping amplitude, t = 1. This model has three phases: (i) a normal phase at high temperature, (ii) a Mott insulating phase at zero temperature and commensurate filling for U ≫ zt, and (iii) a superfluid phase occurring at low temperature. In the following, we will develop an effective theory which goes beyond static mean-field theory by including temporal (in imaginary time) fluctuations, but which is restricted to a single site and therefore contains no momentum fluctuations.
The local terms of the action on a single site are The subscript 'int' denotes the internal degrees of freedom. This local action correctly describes the physics of the system in case of zero hopping (t = 0), when a factorization over each lattice site is exact. Hence it already contains the low-energy physics occurring deep in the Mott phase. At very high temperature the action is also accurate because, in terms of Feynman's path integrals, the world-lines describing the evolution of particles in imaginary time remain almost straight lines, i.e. few hopping processes are present. Our task is then to replace this action with an effective action that also describes the physics deep in the superfluid phase (when the interaction is weak, compared to the hopping amplitude), i.e., in the regime where the Bogoliubov theory of the weakly interacting Bose gas applies. In terms of the treatment of the weakly interacting Bose gas (WIBG) of Ref. 44, there is reason to believe that such an accurate DMFT-like effective action exists: To leading and subleading order (in the interaction), the normal and anomalous selfenergies are momentum independent at zero frequency, which is precisely the DMFT approximation. Note that in Ref. 44 an explicit small symmetry breaking field was added, which introduced a gap in the spectrum (and removed IR divergences in leading orders). Although this violates the requirement that the spectrum of a superfluid should be gapless, 45,46 it was argued in Ref. 44 that the leading orders of thermodynamic observables are found on short-range distances and provided by Beliaev's diagrammatic technique, while for long-range physics (ie, the long-range wavelength fluctuations of the order parameter) one has to resort to Popov's hydrodynamic theory. Similar considerations hold for our effective action where the gap is not fixed but depends on the condensate φ. The short-range physics can be computed, but the long-range fluctuations are not part of the theory. In particular, the explicit symmetry breaking still leads to a finite condensate density in two dimensions at finite temperature (see below, sec. VII B), which greatly simplifies the theory. However, while the explicit but fixed symmetry breaking field of Ref. 44 made the superfluidnormal transition first order, we expect in our case a second order-transition (see below, sec. VII A 2) with nonuniversal critical exponents different from static meanfield theory because of the dynamical corrections in the hybridization terms which form a one-loop correction to mean-field theory. 32

III. DERIVATION OF THE B-DMFT EQUATIONS
In this section we present the derivation of the B-DMFT equations (the action and the selfconsistency relations). For completeness we present alternative (but equivalent) derivations of the DMFT formalism such as a quantum cavity reasoning (Appendix A) and an effective medium approach (Appendix B) in the Appendix. Here we implement an expansion around the atomic limit, following almost literally the lecture notes by A. Georges,47 and consider B-DMFT as an approximation to the kinetic energy functional. This derivation can to a large extent also be found in the supplementary material accompanying our previous Letter, Ref. 36. The atomic reference system is interpreted as the impurity problem. We use a the coupling constant integration method and introduce source fields (Lagrange multipliers) to constrain the condensate field and the connected Green's function for the normal bosons to their physical values. By doing so, we avoid the collapse arguments associated with an expansion in U mentioned in the introduction of this paper.

A. Expansion parameter
We introduce a parameter α ∈ [0, 1] such that When α = 0, the atomic limit is recovered and the partition function factorizes over all sites. When α = 1 the full hopping is recovered, and this is the model we are ultimately interested in.

B. Source fields and constraining fields
Constraining the normal/anomalous Green's functions and the condensate to specified values can be done by introducing conjugate source fields (Lagrange multipliers) in the action. In order to constrain the condensate to Φ we introduce the source field J, and analogous for the connected Green's function G c with source field ∆. Throughout this document we use the Nambu notation in which Φ † = (φ * , φ), J † = (J * , J) and the individual components of G c and ∆ are given by and We can then explicitly write down the grand potential per site (there are N s sites) which is a functional of the source fields and also depends on the constraining fields, Here δb is the non-condensed part of the operator b given by b = b + δb.

D. Full model
The exact functional of the (local) Green's function and condensate are constructed by using the coupling constant integration method, starting from the atomic limit: By using the stationarity of Ω (α-derivatives of the Lagrange multipliers do not contribute) we get where ω n = 2πn/β are even Matsubara frequencies and ǫ k is the dispersion relation of the non-interacting bosons.
We now arrive at the formal expression for the exact functional Γ = Γ α=1 , with the kinetic energy functional given by Requiring stationarity (δΓ/δφ * i (τ ) = 0, δΓ/δφ j (τ ) = 0) determines the value of the source field conjugate to the condensate (assume a homogeneous condensate over the lattice), Since the condensate is time-independent (and taken real), we drop the τ dependence of J 0 as well. The other stationarity requirement (δΓ/δG c = 0, δΓ/δG c = 0) determines the hybridization function appearing in the impurity action: .
Note that for the case z = ∞, we have identically δb = 0 following Appendix A in Ref.1 and only static mean-field theory exists, which is physically clear: For a thermodynamic condensate to be gapless, the k = 0 component of the Green function should not decay in time. The infinite connectivity of the lattice removes any k dependence and in combination with bosonic statistics one then sees that the decoupling approximation is exact.
E. Approximation to the kinetic energy functional B-DMFT can now be considered as an approximation to the kinetic energy functional. With the single-particle Green's function of the Bose-Hubbard model in the presence of source fields and for arbitrary coupling constants, we can define a selfenergy where σ 3 is the Pauli matrix with ±1 on the diagonal.
The DMFT approximation consists in replacing the selfenergy Σ α for arbitrary α by the impurity model selfenergy Σ 0 . Hence, where we have used that the impurity selfenergy satisfies the Dyson equation and G −1 0 is the bare Green's function. Summing over k, and using the constraint on the local lattice Green's function, we obtain the following relation between G c and the hybridization function: with ζ = ∆ α − ∆ 0 + G −1 c . We used the non-interacting density of states D(ǫ) = 1 Ns k δ(ǫ − ǫ k ) and its Hilbert-transformD(z) = dǫD(ǫ)(z − ǫI) −1 . By introducing its inverse,D(R(g)) = g, the relation above can be inverted (αR(αG c ) = ζ = ∆ α − ∆ 0 + G −1 c ) and yields the hybridization function as a functional of the local Green's function, Inserting the above relation into (19), the lattice Green's function expressed as a functional of G c becomes Equation (13) can now be evaluated with G α c (k)| B-DMFT : An explicit expression for the B-DMFT approximation to K[Φ, G c ] therefore reads where the last term reduces to −ztφ * φ for a constant, homogeneous condensate.

F. Stationarity conditions
It immediately follows from Eq. (25) that the stationarity condition for the condensate is unaltered in the B-DMFT approximation (J 0 = ztΦ), while the stationarity condition for the connected Green's function (δΓ/δG c = 0, δΓ/δG c = 0) reads in the B-DMFT approximation (use R(αG)+αGR ′ (αG) = ∂ α [αR(αG)] and the cyclical properties of the trace), Here we have used again the Dyson equation for the second equality. ApplyingD(.) to both sides of Eq. (26) gives This equation defines the B-DMFT selfconsistency condition. B-DMFT maps the bosonic lattice problem to a selfconsistent solution of an impurity model, whose action (expressed in terms of the full operators b, therefore dropping the term ∂ τ ) now takes the form Assuming K = K * , φ = φ * and dropping all subscripts we can write the action in our final version as where the coefficient κ is given by The solution of the impurity problem yields the condensate Φ (Eq. (8)), the connected Green's function G c (Eq. (9)) and the selfenergy Σ imp of the impurity model. The right hand side of Eq. (27) then defines the local lattice Green's function, which is identified with the impurity Green's function and thus allows to obtain the new hybridization function for the next iteration by using Eq. (26).

IV. DMFT PROCEDURE
Solving the single site impurity model means computing the local Green's function and the condensate Φ for the impurity action given by Eq. (29). The hybridization function ∆(τ ) and the condensate order parameter Φ, which is constant in time, are calculated by a selfconsistent procedure starting from some arbitrary initial values (usually obtained from the non-interacting or static mean-field limit). We then solve the impurity problem and calculate the new updated parameters ∆(τ ) and Φ for the action via the Dyson equation. This procedure is repeated until convergence is reached. The selfconsistency equation for the condensate takes the simple form The new hybridization function ∆(τ ) is obtained in the following way: From Eq. (31) and Eq. (8) we obtain the connected Green's function via We then Fourier transform G c (τ ) to obtain G c (iω n ). In Appendix C we show how to obtain accurate Fourier transforms in spite of a finite number of measurement time-steps. With this we calculate the matrix selfenergy via the Dyson equation Here G −1 0 (iω n ) is the bare Green's function which is related to the hybridization function ∆(iω n ) via The parameterμ = µ− ǫ is chosen such that ∆(iω n ) → 0 in the limit ω n → ∞. As we only consider symmetric density of states ( ǫ = 0) here, we write µ from now on. Eq. (35) allows us to rewrite the Dyson equation as Employing the DMFT approximation that the lattice selfenergy coincides with the impurity selfenergy, i.e. the selfenergy loses its momentum dependence: Σ(k, iω n ) = Σ(iω n ), we can calculate the k-summed (or local) lattice Green's function with where ǫ k is the dispersion of the lattice. For some dispersions ǫ k it may be more convenient to transform the summation over wave vectors k into a integration over the density of states D(ǫ).  (39) and the new value for κ: After an inverse Fourier transform we obtain ∆(τ ) and the selfconsistency loop is closed. We then solve the impurity problem again with the updated action, until convergence is reached.

A. Convergence
In this section we discuss two problems which may occur in the iteration process but can easily be overcome. The first problem is that in the first few iterations, when the solution is still far from the converged result, the Hilbert transformation may diverge for some Matsubara frequencies. The origin of the problem can be understood by writing Eq. (38) for the individual components where ζ = iω n + µ − Σ(iω n ), Σ (Σ) are the diagonal (offdiagonal) components of the matrix selfenergy and we have used thatΣ =Σ * . Obviously, the denominator of Eq. (41) can vanish for certain values of Σ andΣ. If the divergence is due to the statistical errors on Σ andΣ, the problem can be cured by running more accurate simulations. However, another reason for a divergent integral can be that the initial approximations for ∆(τ ) and Φ are unphysical ("too far away" from the converged solution). In almost all cases this problem can be avoided by choosing suitable initial hybridization functions ∆(τ ) and condensates Φ. We found that a good choice for the initial Φ was the static mean-field value, corresponding to ∆(τ ) = 0. For this value of Φ we then calculated G c (τ ) and obtained our initial hybridization function ∆(τ ) using the lattice Hilbert transform. Should the integral still be divergent, which only happens in the first few iterations, one can shift (increase) the value of Σ such that the integral becomes well defined. Another strategy to fix this problem is to apply a damping procedure to the selfenergy in the first few iterations. The second problem is related to the convergence properties of the condensate Φ. Sufficiently far away from the SF-Mott transition we typically need about 10 to 20 iterations to reach convergence. Close to the SF-Mott transition, however, the convergence of Φ is very slow, and sometimes up to 500 iterations are needed to reach a converged solution. With a runtime of a few minutes per iteration this translates into a large numerical effort.
Such slowing down of the convergence close to the phase transition can be overcome by using overrelaxation of the condensate Φ. The updating for Φ is changed from Eq. (8) to where Φ old is the condensate from the previous iteration and α > 1.

V. MONTE CARLO METHOD
To formulate our Monte Carlo Method we start by writing the impurity action Eq. (29) in non-matrix form: With this action we expand the partition function Z = Tr[T e −Simp ] in powers of the hybridization functions F , K, K * and the source fields φ and φ * . This leads to the series where |n denotes the eigenstate with n bosons and m = (m F , m K , m K * , m φ , m φ * ). We can now sample individual diagrams with weight These contributions, illustrated in Fig. 1, can be represented by a collection of m F + 2m K * + m φ = m F + 2m K + m φ * creation and annihilation operators on the imaginary time interval [0, β). Hybridization functions F connect m F pairs of creation and annihilation operators, off-diagonal hybridization function K (K * ) connect m K (m K * ) pairs of creation (annihilation) operators, while m φ (m φ * ) creation (annihilation) operators are linked to source fields φ (φ * ). The configuration is fully specified by additionally giving the occupation n of the impurity at times τ = 0, which in combination with the collection of operators determines n(τ ). 1: (Color online) Diagram corresponding to perturbation orders mF = 1, mK = 1, mK * = 1, m φ = 2, m φ * = 2 and n(τ = 0) = 2. The hybridization function F determines the amplitude for transitions of bosons from the impurity into the normal reservoir, while operators coupling to source fields φ and φ * represent transitions between the impurity and the BEC reservoir. The off-diagonal hybridization functions K and K * , present only in the BEC phase, represent the amplitudes for creating or annihilating two bosons at different times.

A. Updates and detailed balance
An ergodic sampling of all possible diagrams requires the following updates: In order to improve the efficiency additional moves can/should be used. Updates which may improve the efficiency are shifts of operator positions, moves which reconnect the hybridization lines and moves which increase/decrease n by one. For example, in the normal or Mott phase, where φ = K = 0, or close to the SF transition, where φ is small, it is useful to have an additional move which reconnects two F -lines. The reason for this is that update 1 is the most time consuming and no F -lines can be reconnected via source fields φ. Deep in the Mott phase the insertion of F -lines is strongly suppressed. In order to have an efficient sampling, we now need an update which increases/decreases the occupation in the whole imaginary time interval without inserting new operators. A graphical representation of the updates is shown in Fig. 2.
In order to satisfy detailed balance, we decompose the transition probability to go from state i to state j as where p prob (p acc ) is the probability to propose (accept) this move.
The only move which changes the number of operators in the imaginary time interval, and therefore the total perturbation order m, is the move to insert or remove an F -line. If we choose to insert an operator pair (Fline) at random times τ and τ ′ drawn from a uniform distribution [0, β) and propose to remove this pair with probability 1/(m F + 1) the proposal probability becomes The factor 2 comes from the fact that we have to decide if the occupation number is changed in the interval between the two operators (with length |τ − τ ′ |) or in the interval which winds around τ = β (with length β − |τ − τ ′ |). If we choose the interval which winds around τ = β this will change the occupation n at time τ = 0. Denoting the factor n| . . . |n in Eq. (45) by w T r (n; τ F 1 , . . . , τ F mF , τ ′F 1 , . . . , τ ′F mF ; . . .), the detailed balance condition for inserting/removing a pair where n ′ and n are different only if the new F -line spans over time τ = 0 (and τ = β).
For moves which increase/decrease the occupation of the impurity by one we have p prop = 1, so the acceptance probabilities satisfy Moves which interchange hybridization lines are easiest, since the operator trace term w T r does not change. The probability to propose such a move is given by the inverse of the number of possibilities to choose the line (or lines) to exchange. For moves which exchange For moves which change an off-diagonal hybridization line into two condensate lines we have to consider that the K-lines do not have a defined direction like the Flines. Therefore the probability to choose one K-line in a state with m K K-lines is given by p prob = 2/m K . Similarly, the probability to choose two φ-lines in a state with m φ φ-lines is given by p prob = 2/(m φ (m φ − 1)). The detailed balance condition for moves which exchange and similarly for the complex conjugate.
For the condensate the functional derivative yields where A MC means that the quantity A should be averaged over all configurations obtained in the Monte Carlo sampling. Due to the time independence of φ this reduces to One can see from Eq. (54) that one needs a condensate in order to measure a condensate. Therefore one always chooses φ = 0 as initial value in the simulation. The diagonal and off-diagonal Green's function can be measured in two different ways. Either by taking the functional derivative with respect to F (τ − τ ′ ), K(τ − τ ′ ) and K * (τ − τ ′ ) or by differentiating with respect to both φ * (τ ) and φ(τ ′ ). This yields in terms of the hybridization and in terms of the condensate with and similarly for the adjoint. The end point G(β − ) can be accurately measured through the density, and G(0 + ) = G(β − ) − 1. In the condensate phase the expansion order of φ and φ * is much higher than the expansion order of the hybridization, see Sec. V D. Therefore, the Green's functions are measured by using the condensate. Close to the SF transition, where φ is very small, it is more efficient to measure the diagonal Green's function according to Eq. (55). In the normal or Mott phase, where φ = 0, the diagonal Green's function can only be calculated in this way. In practice the off-diagonal Green's function is always measured with φ, even close to the SF transition. This is because the expansion order of K and K * is always lower than m φ . The kinetic and potential energy can also be evaluated easily and very accurately. The potential energy is given by and can be computed directly in the Monte Carlo simulation. As shown in Sec. III D the kinetic energy in the DMFT approximation is given by E kin = dΓ dα α=1 . Equation (24) in the presence of a constant, homogeneous condensate then yields By using Eq. (22) and ∆ α=1 = 0 we obtain Going back to imaginary time, we can write this in the individual components as or in terms of the full Green's function By plugging Eq. (54) and Eq. (55) into the above expression we see that the kinetic energy is directly related to the expansion order via where m tot = m F + m K + m K * + (m φ + m φ * )/2 is the total number of operator pairs in the interval [0, β). This algorithm is a bosonic version of the hybridization expansion method of Ref. 19. Aside from the condensate term and the anomalous hybridization, the main difference lies in the fact that for bosons we don't use any summation of diagrams. In the fermionic case, the analytical summation of all diagrams corresponding to a given sequence of creation and annihilation operators (i.e., all diagrams obtained by linking these operators in different ways by hybridization functions) results in a determinant, and was the essential step which allowed to suppress the sign problem and to formulate an efficient algorithm. In a bosonic model, a similar summation of diagrams would lead to a matrix permanent whose evaluation is #P-complete. The only efficient way of evaluating permanents is their stochastic sampling, which is exactly our algorithm of individually sampling the different diagrams (as illustrated in Fig. 1) instead of summing them up expicitly.
Since the off-diagonal hybridization function K is negative, some diagrams have negative weights which leads to a sign problem. However, as we will show later, this becomes an issue only at very low temperatures in the presence of a condensate and does not prevent an accurate computation of phase diagrams and dynamical quantities.

C. Solver test
To check if the diagrammatic sampling and measurement procedure have been implemented correctly we compare the QMC-result with exact diagonalization. In a Hamiltonian formulation, the impurity model can be written as The hybridization functions F and K in Eq. (29) are related to the hybridization parameters V l and W l through For the test, we represent the bosonic bath by a finite number n bath of levels with creation (destruction) operators a † l (a l ) and energies ǫ l . By Fourier transforming Eq. (66) we get For the comparison we choose one orbital (n bath = 1) and the following parameters: t = 1, β = 1, µ = 20, U = 20, κ = 6, φ = 1, ǫ = 1, V = 1 and W = −0.2. In Fig. 3 we show the diagonal and off-diagonal Green's function calculated with exact diagonalization (ED) and Monte Carlo simulations. The density and condensate order parameter are n ED = 1.6762 and φ ED = 1.07976 for ED, and n MC = 1.67619(3) and φ MC = 1.07974(9) for the Monte Carlo simulations. Note that for this set of parameters the average sign s in the Monte Carlo simulation is s = 0.40695 (7). The perfect agreement with the ED result shows that the diagrammatic sampling and measurement procedure have been implemented correctly. Since the model we considered here contains both diagonal and off-diagonal hybridization terms and since a bath with a finite number of levels is as difficult to treat as any other bath, this serves as a nontrivial test for the Monte Carlo solver.

D. Perturbation order and average sign
In this section we show how the expansion order of the hybridization function scales with interaction and how the sign problem, caused by the off-diagonal hybridization K scales at the Mott insulator (or normal phase) to SF transition. Since we sample the operator configuration in the imaginary time interval [0, β) the perturbation order grows roughly proportional to β in all phases.
In Fig. 4 we show how the mean perturbation order grows as one goes from the normal phase to the SF phase at filling n = 1. In the normal phase we have only contributions from the diagonal hybridization F and the perturbation order m F decreases with increasing U . In the SF phase the perturbation order grows rapidly with decreasing interaction U/t, mainly because of the condensate contribution m φ and the off-diagonal hybridization contribution m K .  At the phase boundary where the condensate vanishes (and therefore K vanishes) the sign problem disappears. Therefore the sign problem does not prevent an accurate computation of phase diagrams. The sign problem is only severe deep in the condensate phase.

A. Non-interacting bosons
For an ideal Bose gas (U = 0), the total number of particles is given by where zt is the half-bandwidth of the lattice. Just like for an ideal Bose gas in continuous space, the chemical potential µ has to be lower than the bottom of the band, µ ≤ −zt in order to keep the number of particles finite. For temperatures T > T c , the condensate vanishes, |φ| 2 = 0 and the above equation determines the density n as a function of the chemical potential µ. For T < T c the chemical potential must be pinned at the lower edge of the band, µ = −zt in order to have condensation. The total number of particles is then a function of the condensate density |φ| 2 . In the non-interacting limit the offdiagonal hybridization functions K and K * vanish and the B-DMFT equations equations become exact. The impurity action now takes the simple form This is a quadratic action that can be solved analytically. The solution for the non-interacting Green's function and hybridization function is given by

B. Static mean-field
One obtains a selfconsistent mean-field theory (the decoupling approximation) by substituting into our Hamiltonian defined by Eq. (1). If one drops the term −φ 2 , which is just a constant shift in energy one obtains the following Hamiltonian (72) where κ = ·,j t = zt is the hopping term summed over the nearest neighbors. For the 3d cubic lattice (z = 6) this just gives the half-bandwidth κ = 6t. This Hamiltonian can be expressed as a matrix in the occupation number basis (truncated at some maximum occupation) and solved by exact diagonalization. One chooses an initial value for the condensate φ and determines φ iteratively by solving φ = b until convergence is reached.
By using B-DMFT we can reproduce the static meanfield results by setting the hybridization function to zero, i.e. ∆(τ ) = 0. Since there is no hybridization Eq. (30) reduces to κ = zt.

VII. NUMERICAL RESULTS
In this section we present results for the Bose-Hubbard model on a 3d cubic and 2d square lattice obtained with B-DMFT. All our results are compared to other methods like static mean-field theory, worm-type quantum Monte Carlo (QMC) simulation on a lattice of up to 40 3 sites 8,10 and with a recently developed numerically exact method on the Bethe lattice. 49 Since the Bethe lattice can be also directly simulated with B-DMFT we will show a direct comparison between B-DMFT and the numerically exact solution.  A. 3d cubic lattice

Kinetic and total energy
In this section we present the results for the 3d cubic lattice. Since DMFT can be considered as an approximation to the kinetic energy functional, Eq. (25), it is interesting to see how the kinetic energy per site obtained from B-DMFT compares to results from lattice Monte Carlo simulations 48 and static mean-field theory. In Fig. 6 we show the kinetic energy as a function of temperature as one goes from the SF to the normal phase. In the case of static mean-field theory, the kinetic energy is just given by E kin = −ztφ 2 , where φ = b is the condensate order parameter. In the ground state regime this gives a good approximation of the kinetic energy. For T > T c , where φ = 0, the kinetic energy vanishes since hopping of the normal bosons is completely neglected. The agreement of the B-DMFT results with the exact QMC results is excellent over the whole temperature range. Only close to T c there is some small deviation. For the QMC simulation we used 10 3 lattice sites, except for temperatures in the range 4 ≤ T /t ≤ 6 where 40 3 sites where used. In Fig. 7 we show the total energy as a function of temperature. The remarkable accuracy of the B-DMFT result for the total energy at all temperatures implies that entropies may also be computed reliably using B-DMFT.
In Fig. 8 we show the finite temperature phase dia-gram (top panel) and the ground state phase diagram (bottom panel) for the first and second Mott lobe of the Bose-Hubbard model on a 3d cubic lattice and compare results obtained with B-DMFT to exact results from lattice QMC simulations, the exact solution for the Bethe lattice with coordination number z = 6, 49 and to static mean-field results. 36 The SF phase is characterized by a finite value of φ = b , while we have φ = b = 0 in the Mott insulating and normal phase. For the calculation of the ground state phase diagram we used βt = 2, which is shown in Fig. 8 to be a sufficiently low temperature. Based on simulations done at βt = 4 and βt = 8 we found that any systematic error is smaller than the statistical error. The excellent agreement between our B-DMFT results and the full solution of the Bose-Hubbard model shows that the Mott-transition is a local phenomenon, well described by a momentum-independent selfenergy and that the condensed bosons are accurately described by a uniform condensate.

Critical exponents
We now examine the critical exponent of the order parameter φ. In Ref. 50 fermionic DMFT was employed to study the superfluid state. The temperature dependence of the condensate order parameter ∆ goes as ∆ ∼ |(T c − T )/T c | β in the vicinity of the critical temperature T c , with β = 1/2 being the mean-field exponent. The symmetry breaking happens on the meanfield level in the two-particle channel given by a spin-up and a spin-down fermion. For any fermionic operator we have that c ≡ 0 such that condensation can occur the earliest in the two-particle channel, e.g., c † ↑ c ↓ or c ↑ c ↓ may be non-zero. For a bosonic operator, a single operator b may already have a non-zero expectation value, b = 0, leading to condensation at the one-body level. The B-DMFT action allows for this, while at the two-body level (i.e., the selfconsistently determined hybridization terms) the one-loop correction to the condensate will change the critical exponents from their meanfield values. The exponents we obtain are therefore not universal but depend on the parameters t, µ, U , β and the lattice dispersion ǫ k . Whenever κ = const, we recover the mean-field exponents. This happens of course in the static mean-field limit (κ = zt) and for non-interacting bosons (κ = −G −1 0 (iω n = 0)), where the mean-field exponent β = 1/2 is exact. In Ref. 49 mean-field exponents were also found for the exact solution on the Bethe lattice.
In Fig. 9 (main panel) we show the condensate order parameter φ as a function of temperature obtained by B-DMFT and compare the results to lattice QMC and static mean-field theory. Close to the critical point the condensate goes as φ ∼ |(T c − T )/T c | β for T < T c . The QMC results are obtained on a lattice with 40 3 sites, which still leads to a rounding of the data compared to the thermodynamic limit. They are based on the k = 0

The weakly interacting Bose gas regime
When bosonic field theories are expanded in U (cf. the weakly interacting Bose gas of Ref. 44), the effect of the chemical potential is non-perturbative. The chemical potential is negative for the ideal Bose gas, and has to change sign in the presence of repulsive interactions. There is an implicit relation between the condensate density n 0 = |φ| 2 and the chemical potential, which follows from the condition of thermodynamic equilibrium, 52 leading to µ = 0| ∂Ĥint ∂n0 |0 , withĤ int the interacting twobody terms and |0 the ground state (finite temperature extensions also exist 44 ). In any expansion order this remains valid and can be worked out to yield which is the famous relation for a gapless spectrum first derived by Hugenholtz and Pines. 53 The Hugenholtz-Pines relation does not hold in our B-DMFT formalism. To start the discussion, let us consider a mean-field approach in which the condensate is determined selfconsistently but the hybridization functions are not. Then it is easy to see that the Hugenholtz-Pines relation is shifted to µ = Σ(ω = 0) −Σ(ω = 0) − zt. If, however, the hybridization functions are also determined selfconsistently (as in the full B-DMFT scheme), then they will also depend on and influence the condensate density (and contain µ in them), leading to the breakdown of the Hugenholtz-Pines relation. Indeed, although Ref. 44 uses a skeleton approach for the description of the WIBG, the discussion of the Hugenholtz-Pines relation is provided in terms of bare Green functions in order not to mix the expansion orders.
However, it turns out that the deviations from the (shifted) Hugenholtz-Pines theorem are minimal. This is shown in Fig. 10 where we plot the relative deviation from the Hugenholtz-Pines relation at a value of U/t = 10 across the SF to normal transition. One clearly sees that the deviation grows as one moves deeper into the SF phase. In the Mott and normal phases, where φ = 0, the Hugenholtz-Pines relation is not valid.
As already mentioned in the introduction, it was shown in Refs. 45,46 that gapless field theories have selfenergies that should be momentum dependent and, in particular, that the anomalous selfenergy has to go to zero for small momenta at zero frequency. In Beliaev's diagrammatic approach of the weakly interacting Bose gas, the first order selfenergies are momentum and frequency independent, and this condition is not satisfied, questioning Beliaev's approach. Similarly, in the B-DMFT scheme the selfenergies are momentum independent. In Fig. 11 one sees that the anomalous selfenergy is non-zero for zero frequency, as could be expected from the explicit sym- metry breaking terms in the B-DMFT action. On the other hand, the arguments of Ref. 44 for the weakly interacting Bose gas also hold for B-DMFT. In particular, the authors of Ref. 44 argued that the system reaches its thermodynamic limit for properties such as the energy and the entropy for small system sizes. In Figs. 6, 7 we have already seen that the (kinetic) energy is reproduced remarkably well over the entire temperature range in 3d, except in the fluctuation region near the normalsuperfluid phase transition point.
From the viewpoint that DMFT and B-DMFT sum up all skeleton diagrams built with local propagators only (which is an (asymptotic) expansion in U ), we expect that our B-DMFT theory should recover the physics of the weakly interacting Bose gas (Bogoliubov Hamiltonian), which is known to be an excellent approximation at low values of nU . It is the only limit for which it is not obvious that the B-DMFT formalism will be successful, while the other regimes of strong U or high temperature are well captured by a single-site action. In the skeleton approach to the WIBG theory, 44 the selfenergies are to leading and subleading order given by Σ = 2 n U andΣ = bb U (see Fig. 15). To this order of accuracy, the chemical potential is given by µ (1) = 2 n U − bb U ≈ n 0 U . The latter approximation is made to keep the equations simple (it is not fundamental) since we don't want to go into technical details of the WIBG theory here 44 but it allows us to simply relate the condensate n 0 to the chemical potential µ (1) ≈ µ in this limit, after also taking the shift −6t from the band edge into account. Because the chemical potential in B-DMFT is an input parameter and fully renormalized from the start, this is another way of understanding the deviations from the Hugenholtz-Pines relation.
We now discuss the numerical results for the weakly interacting Bose gas limit. In order to see some depletion effects, we choose parameters such that n = 0.5 and βt = 4. The main result can be seen in Fig. 12, where one sees that the selfenergy and the anomalous selfenergy become frequency independent and approach their WIBG values, which are precisely the leading terms in the high-frequency expansions used in B-DMFT: For the normal selfenergy it is Σ = 2 n U , while for the anomalous selfenergy it is U bb .
In Fig. 13 the chemical potential is calculated such that the density is n = 0.5. The chemical potentials found by B-DMFT are in excellent agreement with the ones found by quantum Monte Carlo worm-type simulations. This shows again that the deviation from the Hugenholtz-Pines relation in B-DMFT is minimal, and goes to zero as U/t → 0. We also compare the condensate density for these theories in Fig. 14, a more challenging quantity than the total density. The agreement with quantum Monte Carlo simulations and the agreement with the WIBG theory for U/t < 0.5 is remarkable, and the difference with static mean-field theory is notable. Technically, calculations in this limit are hampered by the low values of the selfenergies and the proximity of the band edge due to which very precise calculations are needed. For the density n = 0.5 the sign problem is not the lim- iting factor (the average sign is approximately s = 0.7), but the situation deteriorates quickly for higher densities at these temperatures.
To conclude this section, we briefly comment on the difference between the anomalous selfenergy found here (Σ (1) = bb U ) and the one found in Ref. 44 (Σ (1) = n 0 U ) for a system in continuous space. The difference between the two is the diagram given by the convolution of the anomalous Green function with the interaction line, at finite momenta (that is the second diagram forΣ (1) in Fig. 15). In Ref. 44 this diagram was not present in first order for the anomalous selfenergy since it was argued that interactions can to leading order be replaced by their k = 0 values (Eq. (3.44), for which this diagram is zero. The diagram does contribute to the chemical potential however for non-zero momenta, see Fig. 9 on page 22 of Ref. 44. In case of B-DMFT however, we cannot make the distinction that the strength of the interaction is different for small or large momenta. Hence, to leading order there is still a contribution of the anomalous Green function to the anomalous self energy; in other words U bb is the correct value for the selfenergy in the WIBG limit for B-DMFT to leading order.

B. 2d square lattice
In this section we present results for the 2d square lattice. As in the previous section data from B-DMFT are compared to lattice QMC, 54 static mean-field theory and an exact solution of this model on a Bethe lattice with connectivity z = 4. 49 In Fig. 16 we show the finite temperature phase diagram (top panel) and the ground state phase diagram (bottom panel) for the first and second lobe on the 2d square lattice. The ground state phase diagram was calculated at βt = 8 which is sufficiently low as can be seen from the top panel of Fig. 16. As one might expect, the agreement of B-DMFT with exact Monte Carlo lattice simulations is not as good in 2d as in 3d. We note here that in Fig. 16 we compare the condensate transition obtained by B-DMFT to the SF tran- sition obtained by Monte Carlo lattice simulations. Since B-DMFT is just a dynamical extension to static meanfield theory we obtain (as in the static case) a finite condensate φ which is not present in a 2d system except at T = 0 according to the Hohenberg-Mermin-Wagner theorem. Similar as in the discussion of Ref. 44 on the weaklyinteracting Bose gas theory, Beliaev's diagrammatic technique is useful in computing thermodynamic observables since they converge on short-range distances. The longrange fluctuations of the phase that ultimately destroy the condensate can be described at long wavelengths with a separate hydrodynamical action. The diagrammatic technique simplifies hence enormously by working with explicit symmetry breaking and a genuine condensate in the calculation of thermodynamic observables. 44 These ideas are illustrated for the kinetic energy in Fig. 17, where at high and low temperatures the agreement with quantum Monte Carlo worm-type simulations is excellent, in contrast to static-mean field theory. Only in the vicinity of the transition point relatively small deviations are found. As in 3d, B-DMFT predicts a kink in the kinetic energy curve at the transition point, and is hence not well equipped to describe Kosterlitz-Thouless phase transitions. Still, it gives a substantially improved finite temperature phase diagram compared to static meanfield theory (top panel of Fig. 16) and a ground state phase diagram (bottom panel) in good agreement with the QMC result. tivity z the density of states is given by with V the volume of the system. The continuum part D c reads where t is the hopping between the sites. The isolated state at ǫ = −zt matters for condensation and prevents the existence of Goldstone modes in the symmetry broken phase. 55 In Fig. 18 we show a comparison of the density n and the condensate order parameter φ between B-DMFT and the exact solution for a Bethe lattice with connectivity z = 4. As one can see in Fig. 19, the agreement is excellent except for the condensate φ very close to the SF to Mott transition. This comes from the fact that B-DMFT does not reproduce the mean-field exponents like the exact solution 49 but has a smaller exponent in this case, i.e. the transition is much sharper. This excellent agreement of the B-DMFT condensate φ with the exact solution also shows that the definition of the condensate used both in this work and in our previous Letter, Ref. 36, is correct, unlike the condensate defined in the Comment 56 on our previous work, 36 which corresponds to the blue diamonds in Fig. 19.

Limit of infinite connectivity
For systems without symmetry breaking, the isolated state in the non-interacting density of states is negligible. It is customary in (fermionic) DMFT to rescale the  hopping in the limit z → ∞ as t =t/ √ z in order for the density of states to remain finite. Equation (76) reduces then to the semi-circular density of states given by For this density of states the selfconsistency equations simplify since the Hilbert-transformation can be performed analytically. This allows us to formulate the whole selfconsistency loop in τ -space, i.e. we do not have to go to Matsubara frequencies. We obtain a simple relation between the hybridization and the Green's function given by For bosons, however, the isolated state in the noninteracting density of states is not negligible. As soon as z > 2 the condensation temperature is non-zero, and one has then to keep track of the uniform mode. 55 Integer scaling t =t/z is then the only possibility, and the static mean-field theory 1 becomes exact for z → ∞ (for all temperatures). This integer scaling is also compatible with Eq. (30) such that the hybridization terms become identically zero for z → ∞, again in line with Refs. 1,55 but in contrast with the different scalings postulated in Ref. 33.

VIII. CONCLUSIONS
In summary, we have derived the B-DMFT formalism in a number of different ways: (i) by approximating the kinetic energy functional, (ii) by considering a cavity method or 1/z expansion, (iii) by using an effective medium approach. We highlighted similarities and differences with previous approaches. [33][34][35]37,49 We discussed the Monte Carlo impurity solver in great detail including detailed balance of the updates, technical difficulties one may encounter, and the measurement procedure. We have shown results for the 3d Bose-Hubbard model, where the zero temperature phase diagram and the finite temperature phase diagram are in excellent agreement with quantum Monte Carlo worm algorithm results. We compared with the theory of the weakly interacting Bose gas in the regime of weak repulsion. The kinetic energy is reproduced at the 1 − 2% level over the full range of temperatures, except in the very close vicinity to the superfluid-to-normal transition where deviations are larger. We studied the critical exponents of this transition on the superfluid side, and found non-universal values depending on the parameters of the Bose-Hubbard model that are different from the mean-field values and from the universal ones found in the 3d XY universality class. The reason is that the contribution of the hybridization terms can be understood as a one-loop correction to the condensate. In two dimensions, the ground state phase diagram is still in good agreement with quantum Monte Carlo worm simulations, but at finite temperature the superfluid-to-normal transition is less accurately reproduced, because of the Kosterlitz-Thouless nature of the phase transition: in B-DMFT there is still a genuine condensate at finite temperature in 2d, and we miss the long-range fluctuations of the phase. However, quantities that are local such as the energy are still remarkably accurately reproduced over the entire parameter regime, just as in 3d. We also compared B-DMFT with the exact solution of Ref. 49 on the Bethe lattice, and found again good agreement except in the very close vicinity of the phase transition points. The agreement over the whole parameter regime (except in the vicinity of phase transitions) was argued to be the result of fast diagrammatic convergence of thermodynamic quantities with system size, similarly as for the weakly interacting Bose gas system. 44 Our approach constitutes thus a natural extension of the Bogoliubov mean-field theory of the weakly interacting system (which is unable to find a Mott insulator) to a dynamical mean-field theory successful over the entire parameter regime of the Bose-Hubbard model. It shows that, when technical difficulties such as the asymptotic behavior of the series are properly taken care of or sidestepped, elusive diagrammatic expansions for bosonic systems are in fact promising, in line with the spirit of (bold) diagrammatic Monte Carlo 57 studies, whether or not in combination with DMFT. 27 The virtue of developing a DMFT formalism for the Bose-Hubbard model lies in the possible extension of the formalism to models where known numerical methods fail. Prime examples are Bose-Fermi mixtures in two or three dimensions, [58][59][60][61] where previous calculations treated the spinless fermions within DMFT and the bosons at the mean-field level. 62 The formalism can also be generalized to mixtures with spinful fermions. Cluster expansions are possible along the lines of Ref. 43, but it remains to be seen whether the sign problem remains tolerable in the condensed phase. Developing the formalism for real-time applications can be considered along the same lines as for fermionic DMFT, 63 but here, too, it remains to be seen how the sign (as well as a complex condensate) behaves.
The connected part of the two-point correlator satisfies (A5) We now plug Eq. (A4) into Eq. (A3) and obtain In the limit of large connectivity, one has to scale the hopping as t ∼ 1/z, hence zt ∼ 1 and t 2 z ∼ 1/z.
We now proceed to the equations for the true fields, which cannot be found in Ref. 49. 65 To this end, the true marginal η(b) = 1 Z w(b)e zΓ(tb) is needed. It is still expressed in terms of the cavity field (and not the true field) and is obtained by replacing z − 1 by z in Eq. (A6). We try to find now suitable expressions for the true connected Green's function and condensate, valid for the impurity problem 'imp' instead of for the cavity problem 'cav'. The action for the impurity problem is valid up to order 1/z. Since the prefactor of the cavity Green's function is already of O(1/z), we can identify the impurity Green function with the cavity Green function without loss of accuracy, where the average can now be taken with respect to the impurity action. The prefactor of the condensate in Eq. (A6) is of order unity, hence we need the impurity condensate to O(1/z). We have The condensate is now found for the total impurity action S imp , which can be inverted (accurate up to this order), By plugging this equation into Eq. A8 a closed set of equations for the true condensate is found, which is tantamount to B-DMFT. The iteration scheme proceeds in the usual way. Note that the last step does not rely on the scaling in Eq. (78) and Eq.(77), but does use the recurrence relation of a tree. We note that the derivation presented in Ref. 37 is conceptually equivalent to the one presented in Ref. 49, on which this appendix is based.

Appendix B: Effective medium approach
In this section we will expand the full Bose-Hubbard action around a single site, identify the low-energy modes, and re-exponentiate them. The derivation in this section is similar to the derivation found in Ref. 33, but differs in the treatment of the condensate. We provide a microscopic RG-like picture for the condensate field and hybridization fields in the effective action. The total action of the full Bose-Hubbard model is split into three parts, with the first one being the local part defined in Eq. (2), We had that S = S int +S ext +∆S, but would like to derive an action S ≈ S imp + S ext by expanding the exponential function with respect to the action ∆S, and perform the functional integral over the external degrees of freedom, and incorporate this into a new effective internal action S imp and a decoupled, external action S ext . We assume that 'ext' is a thermodynamic system that can spontaneously break the symmetry, if needed, and develop a condensate. If so, we decompose where b ext = φ ext represents the condensate wave function. It is a classical field (c-number, taken to be real) that breaks the symmetry, and its presence forces us to consider anomalous averages such as b ext b ext , which are nonzero. The commutation relations are now [δb (j) ext , δb †(k) ext ] = δ j,k and zero otherwise. The impurity part 'int' is by itself small and cannot spontaneously break the symmetry and develop a non-zero expectation value for its local operators. In our perturbative expansion in ∆S an atom from the 'int' part will hop to the 'ext' part and acquire an expectation value through the full correlators in 'ext'. A non-zero expectation value corresponds to the condensate mode, and is not possible for the bosons outside the condensate. We do not attempt here to build in all lowenergy modes via 'ext', only the classical field contributions from the condensate will suffice. Let us anticipate in advance on such a possibility and allow We also have the commutation relations [δb (j) int , δb †(k) int ] = δ j,k and [δb (j) int , δb †(k) ext ] = 0. In second order, we then also have to allow for b int b int to become large. The validity of Eq. (B5) can then be checked a posteriori.
The coupling between the internal and external degrees of freedom, expressed in Eq. (B2), can then be written as where we have omitted terms in φ int φ ext since the condensate is treated as a c-number and such terms involve an arbitrary shift of the condensate energy. There are four terms in ∆S, two of which have the condensate directly in them, and two that don't. Since we are interested in adding the physics of the WIBG to the action, we have to treat the condensate as a large contribution compared to the small contributions coming from the depleted atoms, while the third term will be contributing to one (and higher) loop corrections. We can now immediately write down the contribution of the parts with the condensate lines to the action, S 1 = ztφ ext δb(τ ) + δb † (τ ) = ztφ ext b(τ ) + b † (τ ) (B7) and we keep only the following terms in ∆S (the other one factorizes), which we treat as a small term, The advantage of our approach is that we have first added the condensate to the impurity action, before looking at depletion (fluctuation) terms. In this picture, the stationarity and time-independence of the condensate is immediately built in. We are now in a position to expand in ∆S. This results in an infinite series, We start with the evaluation of the first order term, δb † int δb ext ext + h.c., (B10) which is zero, as could have been expected from symmetry.
The second order term is the lowest order term that survives, The notation (j) and (k) denotes different sites that are coupled to the impurity through the hopping term in the Bose-Hubbard Hamiltonian. The anomalous terms survive in the presence of a condensate, which tells us that there are four terms in second order. Let us analyze S 1 in more detail, This describes the anomalous process of a depleted boson hopping from the impurity site to the environment, propagating there according to the full (but unknown to us) anomalous two-particle Green's function, and finally being annihilated on the impurity site. We will now perform a cumulant re-exponentiation of this term in the presence of a condensate, for which only the connected diagrams can be kept, where K will be a function which we treat as unknown and which originates from the connected diagrams only. The DMFT selfconsistency iteration scheme will determine K selfconsistently. In a numerical algorithm it is practical to work with full operators b instead of δb, so we will write b int (τ ) − φ int instead of δb int . The other terms S j , j = 2, 3, 4 are treated in the same way as S 1 , but we will not write this out explicitly. After convergence, we wish that our site is in equilibrium with the surroundings. We thus impose a first selfconsistency relation The meaning of this relation is that in every iteration we have to compute b and equate the condensate with this value. It can always be taken real and time-independent.
We collect now all terms needed for our effective impurity action, where the elements of the hybridization matrix are (dropping the subscript 'c' for connected). The remaining hybridization term can now be treated similarly as the hybridization term in fermionic DMFT, and has hence a similar selfconsistency relation 13 (see Sec. IV).