---
abstract: 'We calculate the momentum dependent spectral function of the Bose-Hubbard model on a simple cubic lattice in three dimensions within the bosonic dynamical mean-field theory (B-DMFT).  The continuous-time quantum Monte Carlo method is used to solve the self-consistent B-DMFT equations together with the maximum entropy method for the analytic continuation to real frequencies.  Results for weak, intermediate, and strong interactions are presented.  In the limit of weak and strong interactions very good agreement with results obtained by perturbation theory is found.  By contrast, at intermediate interactions the results differ significantly, indicating that in this regime perturbative methods fail do describe the dynamics of interacting bosons.' 
author:
- Jaromir Panas
- Anna Kauch
- Jan Kuneš
- Dieter Vollhardt
- Krzysztof Byczuk
title: ' Numerical calculation of spectral functions of the Bose-Hubbard model using B-DMFT '
---

Introduction
============

During the last few years impressive progress was made in the experimental investigation of ultracold atomic gases in optical lattices.[@bloch_short; @bloch_long] It is now possible to measure not only density profiles and static correlation functions but even dynamical quantities such as spectral functions of the trapped particles, using Bragg spectroscopy.[@bragg_rmp_2005] Thereby momentum resolved bosonic spectral functions were obtained in the case of condensed bosons[@sengstock_nphys; @bragg_prl_2008; @bragg_3d_njp_2010] and across the phase transition to a Mott insulator.[@bragg_MI_prl_2009; @bragg_fabbri_prl_2012] By contrast, only rather few calculations of bosonic spectral functions have been performed so far for the Bose-Hubbard model.  Early studies were based on the strong-coupling approximation to the Bose-Hubbard model.[@sengupta] Later the weak-coupling limit at zero temperature was thoroughly analyzed with the functional renormalization group.[@dupuis_prl_2009; @dupuis_pra_2009; @kopietz_prl_2009; @kopietz_pra_2010] Other approaches include the variational cluster approximation,[@knap_sf_2011] which was used to investigate systems in one[@knap_1d_2010] and two[@knap_2d_2010; @knap_sf_2011] dimensions, the quantum rotor approximation,[@zaleski] which was recently applied to three-dimensional systems,[@zaleski3d] and the linked-cluster expansion (LCE).[@kauch_lce] In particular the LCE, which was employed to approximately solve the equations of the bosonic dynamical mean-field theory (B-DMFT),[@byczuk_bdmft] allows one to study the Bose-Hubbard model for strong interactions and near the phase boundary between the superfluid and the Mott-insulating phase.  The B-DMFT treats local correlations in time exactly and includes spatial correlations on the mean-field level.  It is applicable for all values of the Hubbard interaction, density and temperature.  The B-DMFT was derived subsequent to the DMFT for lattice fermions[@metz_voll; @kotliar_rmp] and was applied to various bosonic problems[@anders_dmft; @hofstetter_2species_2009] as well as to mixtures[@byczuk_bf; @anders_bf] of bosons and fermions.  It has been extended and applied also to inhomogeneous situations[@hofstetter_2species_trap_2011] and to bosonic systems in non-equilibrium[@eckstein_noneq_2014]. 

In this paper we solve the Bose-Hubbard model within the B-DMFT framework on a simple cubic lattice in three dimensions using a continuous-time quantum Monte Carlo (CT-QMC) solver.[@gull_ctqmc] This approach is known to give excellent results for the phase diagram and static properties[@anders_dmft_short; @anders_dmft] of interacting bosons.  We compute the momentum resolved and the momentum integrated spectral functions and the dispersion relation of interacting bosons in, both, the superfluid and the Mott insulating phase.  The following questions will be addressed and answered: (i) How well can strong-coupling approaches, which are known to capture the phase diagram of correlated bosons very well, describe dynamical properties such as spectral functions? ( ii) How does the presence of the superfluid influence the spectral properties of normal bosons? ( iii) How is the dispersion relation modified by the interaction?  The paper is organized as follows: In Section \[sect2\] we introduce the Bose-Hubbard model and recapitulate the main steps of the B-DMFT framework.  We also introduce an improved method for calculating the self-energy which makes use of two-particle Green functions and discuss the method for the numerical analytic continuation.  In Section \[sect3\] we present the results for the momentum resolved spectral function in the limit of weak, intermediate, and strong interactions, respectively.  Finally, Section \[sect4\] concludes the paper with a summary. 

Model and investigation method {#sect2}
==============================

We consider spinless bosons on a lattice described by the Bose-Hubbard Hamiltonian $$\label{bh_ham}
\hat{H}=-\sum_{ij} t_{ij} \hat{b}^{\dag}_{i} \hat{b}_j - \mu \sum_{i} \hat{b}^{\dag}_{i} \hat{b}_i  + \frac{U}{2}\sum_{i} \hat{b}^{\dag}_i \hat{b}^{\dag}_i \hat{b}_i \hat{b}_i,$$ where $ \hat{b}^{\dag}_i $ ($ \hat{b}_i $) is a bosonic creation (annihilation) operator on a lattice site $ i $, $\mu$ is the chemical potential, $U$ is the local interaction strength, and $ t_{ij} $ is the hopping amplitude.  We assume nearest neighbour (NN) hopping, i.e., $t_{ij}=t>0$ if sites $i$ and $j$ are NN and $0$ otherwise.  Our calculations were performed for a simple cubic lattice with coordination number $z=6$ and $t=0.5$ in arbitrary units. 

The B-DMFT action
-----------------

In the B-DMFT the lattice problem is replaced by a single-site (“impurity”) problem with self-consistency conditions.  A detailed derivation can be found in Ref. .  The impurity action reads $$\label{loc_action}
\begin{aligned}
S_{loc}= & \int_0^{\beta} \d \tau b^{\ast}(\tau)\left(\partial_{\tau}-\mu \right) b(\tau)\\
& +\frac{U}{2}\int_0^{\beta} \d \tau b^{\ast}(\tau)b^{\ast}(\tau)b(\tau)b(\tau)\\
& -\kappa\int_0^{\beta} \d \tau \mathbf{\Psi}^{\ast}\mathbf{b}(\tau)\\
& +\frac{1}{2}\int_0^{\beta} \d \tau \int_0^{\beta} \d \tau' \mathbf{b}^{\ast}(\tau)\mathbb{\Delta}(\tau-\tau')\mathbf{b}(\tau'). 
\end{aligned}$$ Here $\beta=1/T$ is the inverse of the temperature $T$, $\kappa=\sum_{i} t_{ij}=zt$ is a geometrical parameter depending on a lattice type, and $\tau$ is the imaginary (Matsubara) time.  We use the Nambu vector notation for operators $\hat{b}$ and complex variables $b$, i.e., $$\mathbf{\hat{b}}=\left(
\begin{array}{l}
\hat{b}\\
\hat{b}^{\dag}
\end{array}
 \right), \ \ \
 \mathbf{b}=\left(
\begin{array}{l}
b\\
b^{\ast}
\end{array}
 \right).$$ Then the connected Green functions are defined by $$\label{local_green}
\begin{aligned}
\mathbb{G}(\tau)&=-\langle T_{\tau} \mathbf{\hat{b}}(\tau) \mathbf{\hat{b}}^{\dag}(0) \rangle^c\\
&=-\langle T_{\tau} \mathbf{\hat{b}}(\tau) \mathbf{\hat{b}}^{\dag}(0) \rangle+\langle \mathbf{\hat{b}}(\tau)\rangle \langle\mathbf{\hat{b}}^{\dag}(0) \rangle,
\end{aligned}$$ where $\langle\ldots\rangle$ denotes the equilibrium average in the grand canonical ensemble, and $\mathbb{G}(\beta+\tau)=\mathbb{G}(\tau)$.[@Negele] The last two terms in Eq.  represent the coupling of the site to two types of external mean-fields: (i) the Bose-Einstein condensate (BEC) which is represented by a static mean-field $\mathbf{\Psi}$, and (ii) the dynamical mean-field of normal bosons represented by the matrix $\mathbb{\Delta}(\tau-\tau')$, the elements of which are hybridization functions. 

The hybridization functions are related to the Green functions $\mathbb{G}^{(0)}$ for a lattice with a cavity at site $0$ (i.e., where site $0$ is removed) by[@byczuk_bdmft] $$\mathbb{\Delta}(\tau-\tau')=\sum_{i,j\neq 0} t_{i0}t_{j0} \mathbb{G}^{(0)}_{ij}(\tau-\tau').$$ Similarly, the condensate field $\mathbf{\Psi}$ for a lattice with a cavity at site $0$ reads $$\mathbf{\Psi}=\langle \mathbf{\hat{b}}\rangle^{(0)}.$$ The connection between the condensate field on a full lattice, $\mathbf{\Phi}= \langle \hat{\mathbf{b}}\rangle$, and the condensate field on a lattice with a cavity, $\mathbf{\Psi}$, is expressed by $$\label{self_con_phi}
\mathbf{\Psi}=\left[\mathbf{1}
+\frac{1}{\kappa}\int_0^{\beta} \d \tau \mathbb{\Delta}(\tau)\right]\mathbf{\Phi}.$$

Once the mean fields $\mathbf{\Psi}$ and $\mathbb{\Delta}(\tau)$ are known it is possible to solve the impurity problem and calculate the local Green functions and local condensate fields.  This can be achieved by exact diagonalization[@hofstetter_2species_2009; @dmft_ed], the LCE approximation[@kauch_lce], or CT-QMC.[@anders_dmft_short; @anders_dmft] In this paper we use the CT-QMC method to solve the impurity problem.  Our code is based on the work of Anders *et al.*[@anders_dmft] We employ an improved calculation of the self-energy via two-particle Green functions, to be discussed in Section \[imp\_SE\] and in the Appendix.  Tests of our code show numerical agreement with the results of Ref. . 

Self-consistency
----------------

Once the single-particle Green functions have been obtained by solving the single-site problem, the local self-energy can be calculated from the Dyson equation $$\label{dyson_loc}
\mathbb{\Sigma}(\i\omega_n)=\begin{pmatrix}
\i\omega_n + \mu & 0\\
0 & -\i\omega_n+\mu
\end{pmatrix}
-\mathbb{\Delta}(\i\omega_n) - \left[ \mathbb{G}(\i\omega_n) \right]^{-1},$$ where $\omega_n$ are even Matsubara frequencies for bosons ($\omega_n=2\pi n/\beta$).  In the B-DMFT the self-energies are local, i.e., momentum independent.  The momentum resolved Green functions are then given by $$\label{dyson_glob2}\begin{aligned}
 & \mathbb{G}(\mathbf{k},\i\omega_n)=\\
 & \left[\begin{pmatrix}
\i\omega_n+\mu-\epsilon_\mathbf{k} & 0\\
0 & -\i\omega_n+\mu-\epsilon_\mathbf{k}
\end{pmatrix}-\mathbb{\Sigma}(\i\omega_n)\right]^{-1},
\end{aligned}$$ where $\epsilon_\mathbf{k}$ is the dispersion relation for noninteracting bosons.  To close the self-consistency equations we calculate the local Green functions from the $\mathbf{k}$-integrated Dyson equation $$\label{dyson_glob}\begin{aligned}
 & \mathbb{G}(\i\omega_n)= \sum_{\mathbf{k}} \mathbb{G}(\mathbf{k},\i\omega_n)=\\
& \int \d \epsilon D(\epsilon)\left[\begin{pmatrix}
\i\omega_n+\mu-\epsilon & 0\\
0 & -\i\omega_n+\mu-\epsilon
\end{pmatrix}-\mathbb{\Sigma}(\i\omega_n)\right]^{-1},
\end{aligned}$$ where $D(\epsilon)$ is the noninteracting density of states (DOS). 

Efficient calculation of self-energies {#imp_SE}
--------------------------------------

The self-energies can be calculated by solving the single-impurity model with the CT-QMC method and then using Eq.  .  However, this calculation is subject to stochastic uncertainties at high positive and negative frequencies for the following reason: In the first step the impurity problem is solved with the CT-QMC from which one obtains the Green function $G(\i\omega_n)$, which has the form $G(\i\omega_n)=\frac{1}{\i\omega_n} + O(\frac{1}{(\i\omega_n)^2})$ for large $|\omega_n|$. The numerical error of this result is almost independent of $\omega_n$. Therefore relative error grows linearly with $\omega_n$. In the second step the Green function is inverted and substituted into Eq.  .  The linear term cancels out, and only terms of the order of $(\i\omega_n)^0$, $(\i\omega_n)^{-1}$ and smaller remain.  This results in a progressively larger relative error, i.e., the error increases quadratically with frequency $\omega_n$ (see Fig.  \[self\_energy\]). 

To improve the accuracy of the calculation of the self-energies we adapt the method proposed by Snoek and Hofstetter,[@dmft_snoek] and use the two-particle Green functions.  The strategy goes as follows: Finding equations for the self-energies requires solving the equations of motion for the Green functions in the bosonic single-impurity model .  This is similar to the method proposed for fermions.[@bulla] The final expression for the self-energy reads $$\label{imp_est}
\mathbb{\Sigma}(\i\omega_n)=\left[U\mathbb{F}(\i\omega_n) + (\kappa+\mu)\begin{pmatrix}
\phi\phi^{\ast} & \phi\phi \\
\phi^{\ast}\phi^{\ast} & \phi^{\ast}\phi
\end{pmatrix} \delta_{n0} \right] \mathbb{G}^{-1},$$ where $\mathbb{F}(\i\omega_n)$ is the Fourier transform of a matrix of two-particle, *disconnected* Green functions $\mathbb{F}(\tau)$ given by $$\mathbb{F}(\tau)=-\bigg\langle T_\tau
\arraycolsep=1.4pt\def\arraystretch{1.5}
\begin{pmatrix}
[\hat{b}^{\dag} \hat{b} \hat{b}](\tau) \hat{b}^{\dag}(0) & [\hat{b}^{\dag} \hat{b} \hat{b}](\tau) \hat{b}(0) \\
[\hat{b}^{\dag} \hat{b}^{\dag} \hat{b}](\tau) \hat{b}^{\dag}(0) & [\hat{b}^{\dag} \hat{b}^{\dag} \hat{b}](\tau) \hat{b}(0)
\end{pmatrix}\bigg\rangle.$$ Details of this derivation are presented in Appendix \[AppA\].  In Fig.  \[self\_energy\] we see that the improved method significantly reduces the error, i.e., it now grows only linearly in frequency. 

Analytic continuation
---------------------

The main goal of this paper is to calculate the momentum resolved spectral function $A(\mathbf{k},\omega)$ and the momentum integrated spectral function $A(\omega)$, respectively, of the Bose-Hubbard model .  From the CT-QMC impurity solver we obtain the Green functions and the self-energies in imaginary time or in Matsubara frequencies.  Since these functions are analytic in the upper complex plane we can analytically continue them to the real axis.  This involves an inversion of the Hilbert transform: $$G(\i\omega_n) =\int_{-\infty}^{\infty}\d\omega \frac{A(\omega)}{\i\omega_n-\omega},$$ or $$G(\tau) =\int_{-\infty}^{\infty}\d\omega \frac{A(\omega)\e^{-\tau\omega}}{1-\e^{-\beta\omega}}.$$ It is well-known that this is a numerically ill-posed problem, in particular for noisy data obtained with CT-QMC.  To make the analytic continuation tractable we use the maximum entropy (MaxEnt) method.[@skilling_bryan; @bryan; @jarrell] We found that the MaxEnt procedure sometimes fails to converge[^1], especially in calculations of sharply peaked momentum resolved spectral functions $A(\mathbf{k},\omega)$. In these cases we use the *historic* rather than the *Bryan* version of MaxEnt.[^2]

Results {#sect3}
=======

In the following we present results for spectral functions of the Bose-Hubbard model on a simple cubic lattice for weak, intermediate, and strong interaction strengths.  For illustration we indicate the sets of parameters at which the calculations for strong and intermediate coupling, respectively, were performed in the phase diagram (Fig.  \[phase\_diag\]).  The diagram was calculated within the static Fisher mean-field theory [@fisher] and the B-DMFT,[@byczuk_bdmft; @anders_dmft_short] which was solved by CT-QMC[@gull_ctqmc; @anders_dmft] and LCE.[@kauch_lce]

Weak-coupling limit {#non_lim}
-------------------

Our results in the weak interaction limit were obtained for $T=0.5$, $U=0.25$, and $\mu = -2.875$. For these parameters the average occupation per site is $\langle n \rangle=0.4968\pm 0.0005$ and the average number of condensed bosons per site is $\langle b\rangle^2 = 0.4542 \pm 0.0004$. Since $U\langle n\rangle \approx 0.1242 $ is small in comparison to the bandwidth $2zt=6 $, the weak-coupling Bogoliubov approximation[@bogoliubov; @shi; @wibg] can be expected to be applicable.  The momentum resolved spectral function $A(\mathbf{k},\omega)$ obtained within the B-DMFT is presented in Fig.  \[spectr\_weak\]a.  There are two bands: with positive energies for particle addition and with negative energies for particle removal.  Most of the spectral weight is concentrated in the upper band.  The widths of the peaks for specific $\mathbf{k}$-points represent the mean lifetime of the quasiparticles.  The weight of the lower band is orders of magnitude smaller, and therefore its exact shape and position is not determined reliably by MaxEnt.  The dispersion relation $\epsilon(\mathbf{k})$, obtained from $A(\mathbf{k},\omega)$ according to the definition $$\label{disprel_form}
\epsilon(\mathbf{k})=\max_\omega A(\mathbf{k},\omega),$$ is shown in Fig.  \[spectr\_weak\]b.  The dispersion relation obtained from the Bogoliubov approximation is found to be in very good agreement with the B-DMFT result, except near the $\Gamma$ point, where the B-DMFT dispersion does not go to zero.  This is attributed to the fact that the Hugenholtz-Pines theorem[@hugenholtz] is not fulfilled in the B-DMFT as already reported by Anders *et al.*[@anders_dmft]

Strong-coupling limit {#sc_lim}
---------------------

Calculations in the limit of strong interactions were performed for $T=0.5$, and $U=50$. In Fig.  \[spectr\_SCL\] we present the spectral function $A(\mathbf{k},\omega)$ for the Mott insulating phase at $\mu=0.4U$ and the superfluid phase at $\mu=0.9U$ (see Fig.  \[phase\_diag\]).  A striking difference between the spectra in the Mott insulating and superfluid phases, respectively, is the number of bands and their width at specific $\mathbf{k}$-points. 

In the Mott insulating phase there are two bands, separated by a gap whose width is approximately given by $U-3zt$. This value becomes exact in the atomic limit ($t\to 0$).  The factor $3zt$ corresponds to the sum of half-widths of the upper and lower bands.  The shape of this dispersion relation is almost the same as that for a noninteracting band.  In the Mott insulating phase with average occupation $\langle n\rangle\approx 1$ one particle is frozen on every lattice site such that an extra particle or hole can move almost freely through the system.  Since the energy for creating an excitation is large, the density of excitations is very low, and therefore one can neglect their interaction.  As a result quasiparticle and quasihole excitations behave almost like noninteracting particles, i.e., have almost noninteracting dispersion relations.  The width of the peaks is therefore small, implying that the mean lifetime is large.  In fact, we checked that the width corresponds to the uncertainty of the analytic continuation rather than to the mean lifetime. 

The bandwidth of particle excitations is wider than that of hole excitations by approximately $2zt$. In the strong interaction limit the bandwidths of the hole and particle excitations in the Mott insulating phase with integer filling $\langle n \rangle$ are $2zt\langle n\rangle$ and $2zt(\langle n \rangle+1)$, respectively.  This is a quantum effect related to particle indistinguishability, which is simple to derive by starting with a state in which each site is occupied by $n$ particles and treating the hopping term in the Hamiltonian  as a small perturbation. 

In the superfluid phase this is no longer valid.  As the chemical potential increases, the energy for creating particle excitations decreases.  Therefore the interaction between the quasiparticles needs to be taken into account.  Indeed the spectrum of the superfluid is significantly different from that in the Mott insulating phase.  In particular, we find four rather than two bands.  Only the positions of the three lowest bands are determined reliably by MaxEnt.  The four bands arise from two processes: (i) the splitting of a single band due to the interaction, similar to the Mott insulating phase; (ii) the creation of Bogoliubov quasiparticles due to the mixing of particle and hole excitations.  We checked that, in contrast to the Mott insulating phase, the width of the peaks for small $|\omega|$ is robust with respect to MaxEnt parameters as well as re-sampling and therefore represents the mean lifetime of quasiparticles and not the accuracy of analytic continuation. 

In Fig.  \[disprel\_SCL\] we show the dispersion relations $\epsilon(\mathbf{k})$ obtained within the B-DMFT scheme according to Eq. .  We compare them with the dispersion relations obtained from the self-energies calculated within the static Fisher mean-field approximation.[@fisher] The results of both methods are in good agreement in the Mott insulating phase (Fig.  \[disprel\_SCL\], top panel). 

In the superfluid phase (Fig.  \[disprel\_SCL\], bottom panel) this comparison is presented for the three dominant bands.  The results are also in good agreement with the static Fisher mean-field results.  In both approaches the high-energy (negative) band is similar to the band of a noninteracting hole.  The remaining two dispersions are linear for small values of $\mathbf{k}=0$ in the Fisher mean-field theory and correspond to massless Bogoliubov quasiparticles.  In the B-DMFT we also see the linear behaviour except for the vicinity of the $\Gamma$-point.  The dispersion relation around $\omega=0$ does not go to zero, since the Hugenholtz-Pines theorem is not obeyed in the B-DMFT (see also Section \[non\_lim\]).  The slight deviation of the CT-QMC result from that of the static Fisher mean-field theory[@fisher] in the middle of the band around the R point is attributed to the finite resolution of MaxEnt. 

Intermediate interaction and comparison with the strong-coupling solver {#interm}
-----------------------------------------------------------------------

The most interesting regime is that of intermediate interactions where both the Bogoliubov approximation and the static Fisher mean-field theory[@fisher] are no longer valid.  We computed the spectral functions in this regime within B-DMFT and compare our CT-QMC results with those obtained within the LCE[@kauch_lce].  We used the same set of parameters as in Ref. ,  i.e., each set corresponds to a point in a {$\mu /U,t/U$} parameter space for $T=1$. The selected points allow us to study the evolution of the spectral functions through the phase transition.  Here we consider the phase transitions driven by the change of the interaction and by the change of the chemical potential (Fig \[phase\_diag\], left and right panel, respectively).  In Fig.\[lce\_ctqmc\] we present the B-DMFT results and compare with the LCE.  There is good agreement regarding the widths and the positions of the bands in the superfluid phase.  However, in the Mott insulating phase the LCE spectral functions develop a two-peak structure for both positive nad negative energies.  This is not supported by the CT-QMC results.  The origin of this feature in the LCE may be an overfitting of the numerical analytic continuation, since in the LCE the error is small but unknown.  We now focus on the momentum dependence of the spectral functions and the dispersion relations for the Mott insulating and superfluid phases (Figs.  \[Ak\], \[disprel\_int\]).  The most striking feature is the difference of the width of the peaks for specific $\mathbf{k}$-points in these two phases.  Deep in the Mott insulating phase the width represents the uncertainty of the MaxEnt fit.  This is not the case in the superfluid phase in which the width represents the mean lifetime of quasiparticles.  Effectively, deep in the insulating phase we can describe the particles as almost free, whereas approaching the phase transition the quasiparticles obtain a mean lifetime.  In the superfluid phase we observe high-energy excitations, which are shown in Fig.  \[Ak\]c, resembling a structure found in the strong-coupling regime.  Within the accuracy of our method a more precise calculation of their position is not possible.  The dispersion relations shown in Fig.\[disprel\_int\] are compared with the strong-coupling dispersions obtained for $U=18$ and $U=12.5$. In the Mott insulating phase the gap obtained by the static Fisher mean-field theory[@fisher] is much smaller than in the B-DMFT.  Thereby the phase transition is shifted to larger interactions, see Fig \[phase\_diag\].  In the superfluid phase we observe that, in contrast to the weak and strong-coupling limits, the dispersion relation obtained at intermediate coupling within the B-DMFT is not symmetric with respect to zero energy.[^3]

In Fig.  \[spectr\_weak\_U10\] we present the results obtained in the dilute gas regime.  The parameters are $T=0.5$, $U=5$, and $\mu=-2.625$. This corresponds to $\langle n \rangle = 0.12713\pm 0.00001$ and $\langle b \rangle = 0.30301\pm 0.00001$. As in the previous cases (superfluid phase) we observe that the peaks in the spectral function are wide, which means that the mean lifetime of excitations is finite.  The negative energy band has a weight which is orders of magnitude smaller, and therefore its exact shape and position is again not determined reliably by our method of analytic continuation.  The results at the $\Gamma$-point are not included, since MaxEnt failed to converge for very small $\mathbf{k}$.

Although the particle density is low and the majority of the particles constitute the condensate, the interaction is not weak.  Therefore one is not in the parameter range where the Bogoliubov approximation is applicable.  Indeed, as seen in Fig.  \[spectr\_weak\_U10\]b, the dispersion relation is not reproduced by this approximation.  A similar discrepancy was reported in the experiment[@sengstock_nphys], where at high momenta the Bogoliubov approximation was found to overestimate the measured excitation energy.  A further effect of the interaction is the appearance of an additional band at high energies, Fig.  \[spectr\_weak\_U10\]a.  This high energy band appears to repel the one with low energies.  As a result the low energy band is narrower.  This explains the discrepancy between the B-DMFT and the Bogoliubov results, since the high energy band is absent in the latter. 

Summary {#sect4}
=======

We presented a computational method which allows one to calculate spectral functions for bosonic systems described with the Bose-Hubbard Hamiltonian.  This approach has the great advantage of being applicable for arbitrary values of $U/t$, the ratio of the interaction strength and the hopping amplitude.  It reproduces well the results in the limiting cases of large and small values of the interaction.  Results obtained by this method were shown for parameters where perturbative methods fail, thus providing new insights into the properties of correlated lattice bosons in previously inaccessible parameter ranges.  The results show that while in the insulating phase particle and hole excitations tend to behave like free particles, in the superfluid phase the system is described by Bogoliubov quasiparticles with a finite mean lifetime.  Finally, we discussed an alternative method for calculating the local self-energy in the CT-QMC solver which employs two-particle Green functions, and which produces accurate data at large Matsubara frequencies. 

This work was supported by the Foundation for Polish Science (FNP) International Ph.D. Projects Programme co-financed by the EU European Regional Development Fund.  KB acknowledges support by the Foundation for Polish Science (FNP) through the TEAM/2010-6/2 project, co-financed by the EU European Regional Development Fund.  This research was also supported in part by the Deutsche Forschungsgemeinschaft through TRR 80 (DV). 

Calculation of bosonic self-energy from equations of motion {#AppA}
===========================================================

We consider the Hamiltonian representation[@dmft_snoek] of the impurity model $$\label{app_ham}
\begin{array}{rl}
H_{imp}=&-\mu \hat{b}^{\dag} \hat{b} +\frac{U}{2} \hat{b}^{\dag}\hat{b}^{\dag} \hat{b} \hat{b} +\\
& \sum_l\left(\mathbf{\hat{b}}^\dag \mathbb{V}_{l} \mathbf{\hat{a}}_l + \epsilon_l \hat{a}_l^{\dag}  \hat{a}_l\right) - \kappa\mathbf{\Psi}^\dagger \mathbf{\hat{b}}
\end{array}$$ where $\mathbf{\hat{a}}_l^\dag = (\hat{a}_l^\dag \ \hat{a}_l)$ and $\mathbf{\hat{b}}^\dag$ are vectors in Nambu notation.  Here $\hat{a}_l$ ($\hat{a}^\dag_l$) annihilates (creates) a particle from the bath state $l$, $\hat{b}$ ($\hat{b}^\dag$) annihilates (creates) a particle on the impurity, $\mathbb{V}_{l}$ is the Nambu matrix of couplings between the impurity and a bath state $l$, $\epsilon_l$ is the energy of the bath state $l$, and the vector $\mathbf{\Psi}$ represents the condensate field to which the impurity is coupled.  The bath state energies $\epsilon_l$ and couplings $\mathbb{V}_{l}$ are chosen such that $$\label{delta_def}
\mathbb{\Delta}(\i\omega_n) = \sum_l \mathbb{V}_l (\i\omega_n\bbsigma_3-\epsilon_l\mathbb{1})^{-1} \mathbb{V}_l^\dagger,$$ where $\mathbb{\Delta}(\i\omega_n)$ is the dynamical mean field from Eq.  . 

For this Hamiltonian we consider the disconnected Green functions $$\begin{aligned}
\label{app_G}
\mathbb{G}_{00}^d(\tau)=&-\langle T_\tau \mathbf{\hat{b}}(\tau) \mathbf{\hat{b}}^\dag (0) \rangle , \\
\mathbb{G}_{l0}^d(\tau)=&-\langle T_\tau \mathbf{\hat{a}}_l(\tau) \mathbf{\hat{b}}^\dag (0) \rangle .\end{aligned}$$ Here the indices $0$ and $l$ refer to the impurity and bath states, respectively.  The relation between $\mathbb{G}_{00}^d(\tau)$ and the connected impurity Green function $\mathbb{G}_{00}(\tau)$ defined in is given by $$\mathbb{G}^{d}_{00}(\tau)=\mathbb{G}_{00}(\tau)-\langle\mathbf{\hat{b}}(\tau)\rangle\langle\mathbf{\hat{b}}^{\dag}(0)\rangle.$$

The time evolution of the operators is governed by the Hamiltonian $H_{imp}$ (for brevity we will drop the index henceforth), and we can calculate the derivatives of the Green functions with respect to imaginary time.  For the local Green function we have $$\label{app_dtG00}
\arraycolsep=1.4pt\def\arraystretch{1.5}
\begin{array}{rcl}
\partial_\tau \mathbb{G}^d_{00}(\tau)& =& -\delta(\tau)\bbsigma_3 -\left\langle \begin{pmatrix}
[\hat{H},\hat{b}](\tau) \\ [\hat{H},\hat{b}^\dag ](\tau)
\end{pmatrix}
\begin{pmatrix}
\hat{b}^\dag (0) & \hat{b} (0)
\end{pmatrix}\right\rangle \\
&=&-\delta(\tau)\bbsigma_3+\mu \bbsigma_3 \mathbb{G}^d_{00}(\tau) -U\bbsigma_3 \mathbb{F}(\tau)\\
&&-\bbsigma_3\sum_l \mathbb{V}_{l}\mathbb{G}^d_{l0}(\tau)- \kappa\bbsigma_3\mathbf{\Psi}\mathbf{\Phi}^\dagger,
\end{array}$$ where $[ \cdots,\cdots ]$ denotes the commutator, $\bbsigma_3$ is the diagonal Pauli matrix, $\mathbf{\Phi}^\dagger=(\langle b^\dag\rangle \ \langle b\rangle)$ is a Nambu vector, and the Nambu matrix $\mathbb{F}$ is defined as $$\label{app_F}
\begin{aligned}
& \mathbb{F}(\tau)= \\ 
& -\left\langle T_\tau
\arraycolsep=1.4pt\def\arraystretch{1.5}
\begin{pmatrix}
\hat{b}^\dag (\tau) \hat{b}(\tau) \hat{b}(\tau) \hat{b}^\dag (0) & \hat{b}^\dag (\tau) \hat{b}(\tau) \hat{b}(\tau) \hat{b}(0) \\
\hat{b}^\dag(\tau) \hat{b}^\dag(\tau) \hat{b} (\tau) \hat{b}^\dag (0) & \hat{b}^\dag(\tau) \hat{b}^\dag(\tau) \hat{b} (\tau) \hat{b} (0)
\end{pmatrix}\right\rangle;
\end{aligned}$$ a similar relation holds for $\mathbb{G}^d_{l0}$: $$\label{app_dtGl0}
\arraycolsep=1.4pt\def\arraystretch{1.5}
\begin{array}{rl}
\partial_\tau \mathbb{G}^d_{l0}(\tau) &= -\left\langle \begin{pmatrix}
[\hat{H},\hat{a}_l](\tau) \\ [\hat{H},\hat{a}^\dag_l ](\tau)
\end{pmatrix}
\begin{pmatrix}
\hat{b}^\dag (0) & \hat{b} (0)
\end{pmatrix}\right\rangle \\
&=-\epsilon_l\bbsigma_3\mathbb{G}^d_{l0}-\bbsigma_3 \mathbb{V}^\dagger_{l}\mathbb{G}^d_{00}. 
\end{array}$$ To handle the imaginary time derivative we perform a Fourier transform to Matsubara frequencies and obtain $$\label{app_g_omn00}
\arraycolsep=1.4pt\def\arraystretch{1.5}
\begin{array}{rl}
\i\omega_n\bbsigma_3\mathbb{G}^d_{00}(\i\omega_n) =&  \mathbb{1}-\mu \mathbb{G}^d_{00}(\i\omega_n) + \sum_l \mathbb{V}_{l} \mathbb{G}^d_{l0} (\i\omega_n)\\
& +\delta_{n0}\kappa\mathbf{\Psi}\mathbf{\Phi}^\dagger +U\mathbb{F}(\i\omega_n),
\end{array}$$ $$\label{app_g_omn10}
\arraycolsep=1.4pt\def\arraystretch{1.5}
\begin{array}{rl}
\i\omega_n\bbsigma_3\mathbb{G}^d_{l0}(\i\omega_n) =
\epsilon_l\mathbb{G}^d_{l0}(\i\omega_n)+\mathbb{V}^\dagger_{l}\mathbb{G}^d_{00}(\i\omega_n). 
\end{array}$$ Using (\[app\_g\_omn10\]) it is easy to find an expression for $\mathbb{G}^d_{l0}$ in terms of $\mathbb{G}^d_{00}$: $$\label{app_g_omn10_2}
\mathbb{G}^d_{l0}(\i\omega_n)=\left( \i\omega_n\bbsigma_3 - \epsilon_l\mathbb{1} \right)^{-1}\mathbb{V}^\dagger_{l}\mathbb{G}^d_{00}(\i\omega_n).$$

We concentrate now on the case $n\neq 0$, i.e., nonzero Matsubara frequencies.  Inserting into one obtains $$\label{app_int1}
\begin{aligned}
& \mathbb{G}^d_{00}(\i\omega_n)= \\
& \Big[ \i\omega_n\bbsigma_3  +\mu\mathbb{1}
- \sum_l \mathbb{V}_l (\i\omega_n\bbsigma_3-\epsilon_l\mathbb{1} )^{-1} \mathbb{V}_l^\dagger \Big]^{-1}(\mathbb{1}+U\mathbb{F}(\i\omega_n)). 
\end{aligned}$$ Using one can write $$\label{app_int2}
\mathbb{G}^d_{00}(\i\omega_n)= ( \i\omega_n\bbsigma_3  +\mu\mathbb{1}
-\mathbb{\Delta}(\i\omega_n))^{-1}(\mathbb{1}+U\mathbb{F}(\i\omega_n)).$$ Combining with the Dyson equation we arrive at for non-zero Matsubara frequencies (for $n\neq0$ one has $\mathbb{G}^d_{00}=\mathbb{G}_{00}\equiv\mathbb{G}$). 

By inserting into one finds in the case $n=0$: $$\label{app_int1_n0}
\mathbb{G}^d_{00}(0)= \Big(\mu\mathbb{1}
+ \sum_l  \epsilon_l^{-1} \mathbb{V}_l\mathbb{V}_l^\dagger \Big)^{-1}(\mathbb{1}+\kappa\mathbf{\Psi}\mathbf{\Phi}^\dagger+U\mathbb{F}(0)).$$ Starting from for $n=0$ and replacing the disconnected Green function with the connected one leads to $$\label{app_int2_n0}
\mathbb{G}_{00}(0)-\mathbf{\Phi}\mathbf{\Phi}^\dagger= (\mu\mathbb{1}
- \mathbb{\Delta}(0) )^{-1}(\mathbb{1}+\kappa\mathbf{\Psi}\mathbf{\Phi}^\dagger+U\mathbb{F}(0)).$$ In the last step the self-consistency condition , i.e., $\kappa\mathbf{\Psi}=(\kappa\mathbb{1}
+ \mathbb{\Delta}(0)) \mathbf{\Phi}$, is inserted into .  After some regrouping one obtains $$\begin{aligned}
& \mathbb{G}_{00}(0)-\mathbf{\Phi}\mathbf{\Phi}^\dagger= \\
&  (\mu\mathbb{1}
- \mathbb{\Delta}(0) )^{-1}[\mathbb{1}+(\kappa+\mu)\mathbf{\Phi}\mathbf{\Phi}^\dagger+U\mathbb{F}(0)]-\mathbf{\Phi}\mathbf{\Phi}^\dagger. 
\end{aligned}$$ After the above transformations one finally finds $$\mathbb{G}_{00}(0)=(\mu\mathbb{1}
+ \mathbb{\Delta}(0) )^{-1}[\mathbb{1}+(\kappa+\mu)\mathbf{\Phi}\mathbf{\Phi}^\dagger+U\mathbb{F}(0)],$$ which, in combination with the Dyson equation , reduces to for $n=0$.

[45]{}ifxundefined \[1\][ ifx[\#1]{} ]{}ifnum \[1\][ \#1firstoftwo secondoftwo ]{}ifx \[1\][ \#1firstoftwo secondoftwo ]{}““\#1””@noop \[0\][secondoftwo]{}sanitize@url \[0\][‘\
12‘\$12 ‘&12‘\#12‘12‘\_12‘%12]{}@startlink\[1\]@endlink\[0\]@bib@innerbibempty @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [**]{} (, , ) @noop [****,  ()]{} “,”  (, ) Chap.  , pp.   @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{} @noop [****,  ()]{}

[^1]: We use the code which was implemented by Jarrell and Gubernatis following Skilling and Bryan,[@skilling_bryan] and Bryan[@bryan]. 

[^2]: In MaxEnt we use Bayesian inference to specify the probability for finding a spectral function under a constraint of the data.  We then search for a maximum of this probability.  At the same time there exists a parameter which controls the weight with which the entropic prior enters the equations.  In *historic* MaxEnt we decrease the value of this parameter until the spectral function reproduces the data within the uncertainty of measurement.  In *Bryan* MaxEnt we use Bayesian inference to estimate the most probable value of this parameter. 

[^3]: This is not visible in the weak-coupling results, since the negative energy band has a very low weight, such that its position is not determined reliably.  However, the investigation of the self-energy[@anders_dmft] shows that it becomes static in the limit $U\to 0$, which leads to a symmetric dispersion relation.