Self-energy functional theory with symmetry breaking for disordered lattice bosons

We extend the self-energy functional theory (SFT) to the case of interacting lattice bosons in the presence of symmetry breaking and quenched disorder. The self-energy functional we derive depends only on the self-energies of the disorder-averaged propagators, allowing for the construction of general non-perturbative approximations. Using a simple single-site reference system with only three variational parameters, we are able to reproduce numerically exact quantum Monte Carlo (QMC) results on local observables of the Bose-Hubbard model with box disorder with high accuracy. At strong interactions, the phase boundaries are reproduced qualitatively but shifted with respect to the ones observed with QMC due to the extremely low condensate fraction in the superfluid phase. Deep in the strongly-disordered weakly-interacting regime, the simple reference system employed is insufficient and no stationary solutions can be found within its restricted variational subspace. By systematically analyzing thermodynamical observables and the spectral function, we find that the strongly-interacting Bose glass is characterized by different regimes, depending on which local occupations are activated as a function of the disorder strength. We find that the particles delocalize into isolated superfluid lakes over a strongly localized background around maximally-occupied sites whenever these sites are particularly rare. Our results indicate that the transition from the Bose glass to the superfluid phase around unit filling at strong interactions is driven by the percolation of superfluid lakes which form around doubly occupied sites.


Introduction
Ever since the seminal work by Giamarchi and Schulz [1,2] in one dimension and the extension to any dimension by Fisher et al. [3], the intricate interplay of disorder and interactions in bosonic lattice systems has been an active field of research. The advent of cold atom experiments in optical lattices [4,5], where disorder can be realized e.g. through the overlap of optical potentials with incommensurate wavelengths [6,7] or speckle-laser patterns [8,9,10], has further invigorated the interest in this class of systems. More recently, the field has moved to research frontiers such as manybody-localization [11] thanks to advances in monitoring real-time dynamics and state preparation.
The theoretical understanding of disordered and interacting lattice bosons has been primarily driven by numerical simulations. Exact diagonalization (ED) [12] can only be applied to relatively small finite system sizes, while the extension of the density matrix renormalization group (DMRG) [13,14,15,16] to disorder [17] is restricted to lowdimensional systems. In higher dimensions the state of the art method is path integral quantum Monte Carlo (QMC) with worm updates [18,19,20,21]. This algorithm provides numerically exact results for large but finite-sized bosonic lattice systems, while the disorder can be accounted for by averaging over many disorder realizations [21]. However, dynamical quantities such as the single-particle spectral function can only be determined by performing analytic continuation of imaginary-time propagators with stochastic noise [22,23]. The continuation is an inherently ill-posed problem, and cannot resolve sharp resonances. While the methods mentioned above excel with a high numerical accuracy, they rely on finite system sizes, which can represent a problem when rare disorder-driven fluctuations play an important role, which can only be captured once one approaches the thermodynamical limit. The available methods in the thermodynamical limit rely on approximations. The mean-field decoupling approximation [3] can be applied to disordered systems using an arithmetically averaged condensate. However, this overestimates the phase coherence and the extent of the superfluid phases, as locally condensed bosons are mistaken for a global condensate [24]. In fact, in mean-field methods with position-space resolution the geometric percolation of condensed regions appears to be a more accurate quantity to evaluate the global superfluid response [25]. A more accurate mean-field approach is the stochastic meanfield theory [26,27], where the condensate order parameter is treated as a disorderdependent quantity. However, mean-field methods are self-consistent only in terms of the condensate (i.e. the one-point propagator). This severely hampers the ability to describe uncondensed phases, which are simply approximated by the zero hopping limit (i.e. the atomic limit).
A non-perturbative method which includes also a self-consistency in terms of the two-point propagator is the dynamical mean-field theory (DMFT), originally formulated for fermions [28,29] and later generalized to bosons [30,31,32,33,34,35]. For fermions the formalism has been extended to disordered systems [36,37,38,39] by averaging the systems propagators over all disorder configurations. While an arithmetic averaging in this framework works well for weak disorder, it misses the essential physics in non-selfaveraging phases. In such phases, like the Anderson-localized regime [40,37], observables show broad tails in their disorder-distribution. An interesting idea for incorporating nonself-averaging effects is the typical medium theory [36,37,38,39], where the arithmetic average is replaced by a geometrical mean. However, it is not clear what the range of validity is for this approach. We would like to point out that the works above employing disorder and DMFT all study fermionic systems. As of today we are not aware of any works applying DMFT to disordered bosonic systems.
A more general theoretical framework for constructing non-perturbative approximations for interacting many-body systems is the self-energy functional theory (SFT) [41,42,43,44,45,46,47,48], from which DMFT can be derived as a certain constriction of the variational space. The formalism was first developed for fermions [41,42,43,44] and later extended to bosonic systems [45,46,47,48]. Our recent derivation [48] based on the bosonic Baym-Kadanoff functional [49,50] correctly includes U (1)-symmetrybreaking, and simplifies to bosonic DMFT [30,31,32,33,34,35] in a particular limit. Within SFT, non-perturbative approximations are readily constructed by restricting the self-energy domain of the original lattice system to the self-energies of a simpler exactly solvable reference system. This reduces the full complexity of the original problem to a search for stationary solutions in terms of the variational free propagators of the reference system. The generalization of SFT to systems with disorder has been developed for fermions in Ref. [51] and applied in a variational cluster approximation to bosons in the absence of U (1)-symmetry-breaking in Ref. [52].
The aim of this paper is to extend the bosonic SFT formalism of Ref. [48] to disordered lattice bosons including the possibility of U (1)-symmetry-breaking. As argued for fermions in Ref. [51], the geometrical mean used in the context of DMFT [36,37,38,39] is hard to reconcile with the variational SFT framework. We therefore derive an arithmetically averaged formalism, where, through the introduction of an appropriateT P V functional, the functional depends only on the self-energies of the arithmetically averaged propagators. Just as the version for clean systems [48], we find that SFT incorporates a disorder-averaged generalization of bosonic DMFT [30,31,32,33,34,35] in a certain limit. The resulting functional is, however, more general than DMFT by being amenable to a more general variational space.
The prototypical model for interacting disordered bosonic lattice systems is the Bose-Hubbard model (BHm) in the presence of local disorder (for a review, see Ref. [53]). In addition to the Mott insulating and superfluid phases of the clean system, the groundstate phase diagram exhibits a new phase: the Bose glass [1,2,3]. This is an insulating but gapless and compressible phase, which always intervenes between the Mott insulator and the superfluid phase at finite disorder [54,55,21]. While certain single-particle states can show a high (but not macroscopic) occupation, the disorder does not allow for global phase coherence as observed in the superfluid. The statistical fluctuations of the local potential, on the other hand, locally exceed the gaps of adding/removing a particle [21], creating gapless regions which induce a non-vanishing density of states at zero energy [3]. In the groundstate phase diagram of the disordered BHm on a cubic lattice the superfluid phase extends to surprisingly large values of the interaction U and the disorder strength D. For low/intermediate interactions and high disorder it can be argued that this is related to the percolation between localized states [21]. At stronger interaction and lower disorder, the phase diagram is characterized by the so-called "superfluid finger" which extends to much larger interactions than the critical value of the clean system and is characterized by an extremely low condensate fraction. The critical temperature at which the condensate vanishes is thus extremely low, making it very hard to access this regime in experiments [21].
We apply SFT to the BHm with local box disorder on the cubic lattice using the simplest imaginable reference system, comprising a single bosonic mode. This restriction to the minimal reference system produces a self-energy approximation with three variational degrees of freedom, which we will denote by SFA3. In the clean BHm the SFA3 approach has been shown to be in quantitative agreement with numerically exact QMC results [48]. In this work we investigate the disordered BHm in the vicinity of the superfluid finger, where the condensate density is extremely fragile, leading to a substantial shift in the phase boundaries even if the numerical error is very low. Nonetheless, we observe excellent agreement of the thermodynamic quantities computed with SFT and the QMC reference results.
Since the SFA3 reference system can be solved exactly, we can also evaluate the lattice spectral function and thereby obtain spectroscopic information not readily available from QMC. By systematically analyzing the local excitations of the SFA3 spectral functions, we find that the strongly-interacting Bose glass is characterized by different regimes, depending on which local occupations n are activated as a function of the disorder strength D. While local observables are described well by the atomic limit, we find that the particles delocalize into isolated superfluid lakes over the stronglylocalized background around highly-occupied sites whenever these sites are particularly rare. In particular, our results indicate that the transition from the strongly interacting Bose glass to the superfluid phase close to unit filling is driven by the percolation of superfluid lakes which form around doubly occupied sites. As D is further increased and the number of highly-occupied sites increases accordingly, the particles are localized by the increasing particle-number fluctuations and interaction energy, explaining the reentrant behavior of the superfluid finger at larger D.
We also present results deeper in the superfluid phase (i.e. at weaker interactions), showing excellent agreement with QMC for thermodynamical quantities at low disorder. When the disorder dominates both over the bandwidth and the interaction, the restricted variational subspace of our SFA3 reference system is however insufficient, as we no-longer can stabilize a stationary solution. Whether this can be remedied by a more general reference system construction is an open question. This paper is organized as follows. In Sec. 2 we derive the self-energy functional theory for disordered lattice bosons: starting from the free-energy functional (Sec. 2.1), we generalize the bosonic Baym-Kadanoff functional to the case of disorder (Sec. 2.2), perform a Legendre transform to the self-energy effective action (Sec. 2.3), average the effective action over all disorder configurations (Sec. 2.4), and finally arrive at the disorder-averaged self-energy functional (Sec. 2.5). In Sec. 3 we introduce the disordered BHm, discuss the SFA3 reference system used in the SFT calculations (Sec. 3.1) and derive analytic results in the atomic limit (Sec. 3.2). The numerical results are presented in Sec. 4, where we investigate the strongly-interacting Bose glass (Sec. 4.1), the stronglyinteracting superfluid phase transition (Sec. 4.2), and the superfluid phase (Sec. 4.3). Finally, Sec. 5 is devoted to the conclusion.

Self-energy functional theory for disordered lattice bosons
In this section we derive the self-energy functional theory for disordered lattice bosons. In analogy to the formalism for clean systems derived in Ref. [48], we do so by a series of Legendre transformations starting from the free-energy functional and introduce a simpler exactly solvable reference system sharing the same local interaction and disorder distribution. As was done in a previous work on fermions [51], we average over all possible disorder configurations, arriving at a functional which only depends on the self-energies of the arithmetically averaged propagators of the system.
Note that, in order to keep track of the various additional dependencies arising through the disorder, we introduce a slightly more complex notation than in our work on disorder-free bosons in Ref. [48], by denoting the explicit dependencies on system parameters as subscripts. Further, we will denote Nambu objects (i.e. matrices or vectors) as bold letters (e.g. O), scalars as simple letters (e.g. O), and functionals with a hat (e.g.Ô). Finally, for notational simplicity, we denote the one-point self-energy (formerly Σ 1/2 in Ref. [48]) as S.

Free-energy functional
Consider a lattice system of bosons in the presence of quadratic disorder, with creation where α is a superindex spanning both the site index i and the Nambu index ν, the HamiltonianĤ of the system can be written aŝ where repeated indices are summed over, F is an explicit symmetry-breaking field, t α β = t iη jν = t ij ⊗ 1 ην is the hopping, the quadratic disorder η α β = η iη jν = η ij ⊗ 1 ην has the probability distribution P (η), and the general three and four-body interactionV can be expressed asV ≡ V To keep the notation compact we will henceforth suppress the lattice-Nambu superindices.
At finite temperature T ≡ β −1 the free energy functional of the interacting system is given byΩ where the subscript 'V ' means that in addition to F and G −1 0 the functional also depends on the interaction vertexV , and S V is the imaginary-time action The free energy functionalΩ V [F, G −1 0 ] is equal to the free energy Ω FtηV of the lattice system in Eq. (1) with fix disorder configuration η, when evaluated at the symmetry breaking field F and the free single-particle propagator G tη0 , i.e.
where the non-interacting (V = 0) single-particle propagator of Eq. (1) is given by and the subscript means that it depends on the hopping t and the disorder configuration η only. By taking functional derivatives of the free energy functionalΩ V with respect to F and G −1 0 we obtain the two functionalŝ that reproduce the physical one-and two-point propagators (i.e. the condensate Φ FtηV and the connected Green's function G FtηV ) of the disordered interacting system in Eq.
(1) when evaluated at F and G −1

Baym-Kadanoff functional
When exchanging the functional dependence of the free energy functionalΩ V in Eq. (2), from F and G −1 0 to Φ and G by means of a Legendre transformation, one obtains the bosonic Baym-Kadanoff functional [49,50,48] Γ (BK) Here, the entire complexity of the many-body system is contained in the Luttinger-Ward functionalΦ which contains all two-particle irreducible diagrams [56,57]. For a more detailed discussion of the Luttinger-Ward functional in the context of SFT, see Ref. [48]. At the physical interacting one and two-point propagators, the Baym-Kadanoff functionalΓ and equal to the free energŷ The explicit functional derivatives take the form By identifying the variations of the Luttinger-Ward functionalΦ as the one and two-point self-energies [49,50,48] and applying the stationarity conditions [Eq. (10)] we find that the interacting propagators fulfill the two Dyson equations Consider now the result of substituting F and G −1 0 using Eq. (14) in the functionalŝ φ V andĜ V [Eqs. (6) and (7)]. This gives the highly non-linear coupled equationŝ For given self-energies S and Σ the concomitant solution of Eqs. (15) and (16) implicitly defines the functionalŝ depending solely on the self-energies S and Σ and the interactionV , producing the physical interacting propagators when evaluated at the physical self-energies, i.e.

Bosonic self-energy effective action
By means of a further Legendre transform the Baym-Kadanoff functionalΓ FtηV with functional dependence on Φ and G can be transformed into the self-energy effective actionΓ depending on the self-energies S and Σ, where the universal functionalF V [S, Σ] is the Legendre transform of the universal Luttinger-Ward functionalΦ , with variations (see Ref. [48] for details) The variations of the self-energy effective action give FtηV is stationary at the physical self-energies and equal to the free energŷ

Disorder-averaged self-energy effective action
While we up till now have treated a system with a single disorder realization η, we are interested in describing the averaged ensemble of systems with disorder probability distribution P (η) and its ensemble averaged free-energy In terms of the self-energy functional Eq. (19) the averaged free energy can be expressed as using Eq. (23). However, a direct application of the avaraged self-energy functional does not lend itself to the construction of approximations using disorder averaged propagators and self-energies.
To describe the combined effect of disorder and interaction we seek to construct an extended disorder-averaged functional parametrized by the disorder-averaged propagatorsΦ using the short-hand notationĜ η ≡Ĝ V [S η , Σ η ] andΦ η ≡Φ V [S η , Σ η ], where S η and Σ η denote the self-energies for the disorder configuration η, see Eq. (18). The corresponding average self-energiesS andΣ are defined through the Dyson equations where G t00 is the free propagator for the disorder-free system G t00 ≡ G tη0 | η=0 . By insertion of the averaged Dyson equations [Eq. (27)] in the definitions of the averaged propagatorsΦ andḠ [Eq. (26)] we obtain the relations where The concomitant solution of Eq. (28) and (29) implicitly defines the two universal functionalsΦ Using the universal averaged propagator functionals we define the extended averaged self-energy effective action whereT P V is a universal functional of the averaged self-energieŝ Here, we have extended the variational space from the fixed-disorder self-energies However, when evaluated at the physical self-energies,Γ

(SE)
FtP V takes the value of the disorder average of the self-energy functional Γ (SE) FtηV P , see Appendix A, and is thus equal to the disorder averaged free energy Ω FtP V by Eq. (25).
The variations ofT P V are derived in Appendix A and give, δS †T P V =Φ and 2δΣT P V =Ĝ, showing thatT P V is the analogue of theF V functional for the averaged self-energies, as by Eq. (20), δ S † ηF V =Φ η and 2δ ΣηFV =Ĝ η . The functional derivatives with respect to the self-energies at fixed disorder configuration η yield where Q η is defined in Appendix A and vanishes at the physical self-energies. The variations of the averaged self-energy effective actionΓ 2δ ΣηΓ the averaged self-energy effective actionΓ and equal to the average free energŷ The crucial part of the disorder-averaged self-energy effective action is that the functionalsT P V andF V are universal, in the sense that they do not depend on the non-interacting propagator G t00 or the symmetry-breaking field F (see Appendix B for the explicit derivatives ofT P V and Ref. [48] for the universality ofF V ), but only on the interaction V , the disorder probability distribution P (η), the disorder-dependent self-energies {S η , Σ η }, and the average self-energiesS andΣ. In the following we will make use of this property in order to derive consistent approximations ofΓ

Disorder-averaged self-energy functional theory
A versatile approach to non-perturbative approximations of the self-energy effective actionΓ FtηV is the self-energy functional theory (SFT) pioneered by Potthoff [41,42,43,44] for fermionic systems and later extended to bosonic systems [45,47,48]. The formalism for systems with disorder has been developed for fermions in Ref. [51] and applied in a variational cluster approximation (VCA) to bosons without symmetry breaking in Ref. [52]. Here, we generalize the bosonic case to also include U (1)symmetry-breaking and general reference systems.
We consider the general interacting bosonic system with quadratic disorder of Eq. (1), and introduce a second reference system with the same interactionV and disorder P (η) but with some arbitrary linear symmetry breaking field F , arbitrary free propagator and self-energy effective actionΓ Here, the free propagator G ∆η0 is parametrized by replacing the hopping t by a completely general matrix ∆(τ, τ ) ‡ . Now, since the self-energy effective actions of both systems contain the same universal functionalsF P V andT P V we can evaluateΓ The stationary condition in Eq. (41) now translates into 2 δΓ (46) ‡ In e.g. the context of dynamical mean-field theory ∆(τ −τ ) would represent a retarded hybridization of an impurity with a non-interacting bath, while in the case of e.g. an instantaneous ∆(τ, τ ) = t δ(τ − τ ) where t is diagonal in Nambu space, it can be considered to be a hopping amplitude.
If by an appropriate choice of ∆ and F the reference system can be made simple enough to be exactly solvable, one can go one step further and evaluate the original systems functionalΓ (SE) FtP V at the physical self-energies of the reference system, i.e. at S F ∆ ≡S F ∆P V ,Σ F ∆ ≡Σ F ∆P V , and {S F ∆ηV , Σ F ∆ηV }. This produces the selfenergy functional theory (SFT) approximation for the system and the SFT functional where we have used thatΓ is the self-energy effective action of the original systemΓ FtP V restricted to the domain of physical self-energies of the reference system. Note that by replacinĝ Γ with the scalar Ω F ∆P V , we eliminate all explicit dependencies on the fixed-disorder self-energies {S F ∆ηV , Σ F ∆ηV }, such that the disorder-averaged self-energy effective action now only depends on the average self- The domain ofΓ FtP V is therefore defined by the average physical self-energies of the reference system (S F ∆ andΣ F ∆ ) and parametrized by ∆ and F . By generalizing the variational principle of Eq. (41) to the restricted domain we obtain a thermodynamically optimal approximation when the self-energy variations are zero on the domain, i.e. we seek ∆ and F such that which by can be fulfilled if By the SFT approximation, the entire complexity of the original lattice system has therefore been reduced to finding stationary solutions of the functional in Eq. (47) in terms of the variational parameters F and ∆.
Note that by restricting ∆ to be local, but keeping its full imaginary-time dependence, SFT reduces to a disorder-averaged version of bosonic dynamical meanfield theory [30,31,32,33,34], see Appendix C.

Disordered Bose-Hubbard model
As a simple application of the formalism derived in Sec. 2 we study the disordered Bose-Hubbard model (BHm) on the cubic lattice with uncorrelated box disorder. The Hamiltonian has the form where b is the occupation number operator, i, j denotes summation over nearest neighbors, J is the hopping amplitude, U the on-site interaction, and µ the chemical potential. The local disorder potentials η i are uncorrelated and are assumed to have a flat probability distribution where D is the disorder strength. Thus, the free propagator [Eq. (5)] is given by In addition to the superfluid and Mott insulating phase of the clean (i.e. nondisordered) model, the groundstate phase diagram of the disordered BHm exhibits the Bose glass phase. It is an insulating compressible phase that always intervenes between the superfluid and the Mott insulator at finite disorder (D > 0) [21]. The phase is composed of local regions, including both strongly localized atomic levels, and isolated superfluid lakes which locally close the many-body gap. Since these superfluid lakes are spatially separated, global phase coherence is not reached, yielding a zero superfluid response. Hence, while the compressibility of the Bose glass is finite, the global condensate order parameter is zero as is the many-body gap. In our SFT formalism we therefore distinguish the Bose glass and the superfluid phase by the disorder-averaged condensateΦ. In the mean-field approach of Ref. [25] a different criterion was used, analyzing the spatial percolation of superfluid regions (i.e. isolated regions with non-zero quasi-condensates), allowing for a zero global superfluid response even thoughΦ = 0. As in this work we will treat only disorder-averaged translational-invariant quantities, such a real-space percolation of the condensate is not analyzed directly.

Minimal reference system
In this first application of SFT with symmetry breaking to the disordered BHm, we will make use of the simplest possible reference system, comprising a single bosonic mode per site. In this case the reference system Hamiltonian in the thermodynamic limit reads where the sum runs over an infinite number of independent single-site Hamiltonians The reference system is parametrized by three real translationally-invariant parameters F , ∆ 00 , and ∆ 01 , where F = (F , F ) and ∆(τ − τ ) = δ(τ − τ )∆ is instantaneous in imaginary time and site-local as This minimal reference system yields a non-perturbative self-energy functional approximation that we denote by SFA3. It has previously been shown to yield quantitatively correct results for the clean BHm, comparing with numerically exact QMC results [48].
In the case of uncorrelated disorder considered here, the disorder-averaging of observables [Eq. (24)] gives translationally invariant results, see Appendix D. The disorder-averaged free energy of the SFA3 reference system is therefore given by where N is the number of lattice sites, Ω i,F ∆η i U is the free energy of a single site in the reference system, and f (η) p ≡ dη p(η)f (η). Analogously, the propagators are obtained asḠ Hence, to evaluate the disorder-averaged quantities of the reference system it suffices to solve the single-site Hamiltonian of Eq. (57) for all possible values of η i and then average the result over the probability distribution p(η). The corresponding average self-energiesS F ∆P U andΣ F ∆P U of the reference system are then obtained from Eq. (27).
Physical solutions of the lattice system can be found by searching for stationary values of the SFT functional in Eq. (47) fulfilling Eq. (52) using a standard root solver to find the point with zero gradient. This procedure is identical to the algorithm detailed in Ref. [48]. Once a stationary solution is found, lattice quantities can be computed using the corresponding self-energies at stationarity as detailed in Appendix E.

Atomic limit
As in this work we mainly analyze the behavior of the disordered BHm at large interactions U/J 1, we want to compare to the analytic atomic limit of having decoupled sites, i.e. J = 0. In this section we analyze the properties of the infinite system in this limit.

Local occupations
We start by analyzing the local occupations as a function of disorder in the atomic limit. For J = 0, we can have a local occupation n i = n i at zero temperature if the local potential η i takes values η min (n i ) < η i < η max (n i ). In order to derive this, we turn to the local energy of the decoupled site i with occupation number n i and local potential η i , i.e.
The groundstate will have local occupation larger than where we used that the local occupation n i is bounded from below by zero, and therefore η min (0) = ∞. Additionally, in order to have a local occupation of n i we need to fulfill the condition E SS (n i , η i ) < E SS (n i + 1, η i ), resulting in As the minimum possible value of η i is −D [see Eq. (54)], this implies that the maximal possible local occupation n max is given by Furthermore, as the maximum value of η i is D and the local occupation is bounded from below by n i = 0, we have a minimal possible local occupation of We can use the information above to derive the probability of sites with occupation n in the infinite system. We denote this quantity by r n , defined as the number of sites with local occupation n divided by the total number of sites, which can be computed by r n = 0 if n < n min or n > n max , We can use the probabilities r n in order to derive expressions for the total density n and the interaction energy E int , given by Note that, while the values of r n depend on the disorder distribution P (η), the values of D where they become non-zero [and therefore the maximal and minimal possible occupations for a given disorder strength in Eqs. (64) and (65)] depend only on the maximal and minimal values of the local potential ±D (and the global parameters µ and U ). These are therefore universal, in the sense that they do not depend on the disorder distribution P (η) as long as it is uncorrelated and bounded [i.e. with p (|η i | > D) = 0].

3.2.2.
Local excitations On a single-site level, the process (n i → n i + 1) on site i leads to the energy difference As discussed in Sec. 3.2.1, in the groundstate we can have a local occupation of n i only if η min (n i ) ≤ η i ≤ η max (n i ), see Eqs. (62) and (63). If we average over all sites i, we therefore find, that the local processes (n l → n l + 1), with local groundstate occupations n l , span over the energy range given by We can further derive the energy difference for the opposite process in the same way, yielding The disorder-averaged local spectral function is defined as At zero temperature we have where GS is the groundstate, the sum runs over all other eigenstates, E n is the energy of eigenstate n, and ω + = ω + i with a small broadening parameter . Disorder-averaging over an infinite number of configurations therefore yields a translational invariant local Green's function whereñ (η min (n) < η < η max (n)) = n.
Using Eqs. (69) and (70) we therefore find that the resonances of the spectral function for the processes (n → n + 1) are bounded by Eq. (69), while the processes (n + 1 → n) are bounded by Eq. (70). Therefore, the effect of the disorder strength D on the spectral function in the atomic limit -which in the absence of disorder consists of sharp delta peaks -is to broaden the peaks to a width which is proportional to D. A consequence of this is that, apart from the process (0 → 1), all local resonances are bounded by −U ≤ ω ≤ U . Further, it can easily be shown that for D ≥ mU/2, with integer m, we have ω max (n → n + 1) > ω min (n + m → n + m + 1), leading to an overlap of the processes (n → n + 1) and (n + m → n + m + 1) (and equivalently for the reversed particle-removal processes)

Results
We analyze the BHm with box disorder using SFT with an SFA3 reference system, see Sec. 3.1. The calculations are compared to disorder-averaged path integral quantum Monte Carlo (QMC) [20,58,21] simulations on a finite cubic lattice of 8 3 sites. In the strongly-interacting case we further compare to analytic results in the atomic limit (i.e. the limit of zero hopping J = 0), detailed in Sec. 3.2. The resulting groundstate phase diagram computed with SFA3 at large interactions for fixed density n = 1 is shown in Fig. 1 together with the QMC results of Ref. [21]. The groundstate phases exhibited by the system are the superfluid, the Mott insulator, and the Bose glass. For the ordered BHm (η = 0), the SFA3 approximation showed remarkable agreement with exact QMC results [48]. The phase diagram in Fig. 1 shows that this remains true also for weak disorder D/J 30, where the SFA3 superfluid to Bose glass transition line shows excellent agreement with the QMC result.
For stronger disorder the situation changes, in particular in the so-called superfluid finger, i.e., the narrow region of the superfluid phase extending to large values of the relative interaction strength U/J. In the finger, the condensate density ρ c = 1 2Φ †Φ is extremely low, and therefore very hard to resolve experimentally [21]. Small deviations from numerically exact results in the SFA3 calculations therefore lead to a notable shift in the phase boundaries and an overestimation of the extent of the superfluid finger, as seen in Fig. 1. At even larger disorder when leaving the superfluid finger, the discrepancy between SFA3 and QMC results is reduced. The Mott insulator to Bose glass transition at fixed density n = 1 is very hard to resolve numerically (unlike the transition at fixed chemical potential discussed later), as the finite compressibility in the Bose glass close to the phase boundary is exponentially small [21]. Instead, in Fig. 1 we show analytic results on the phase boundary from Ref. [21].
Note that, while it is always possible to find a Bose glass stationary solution at strong disorder in mean-field approaches (by setting the condensate to zero), arithmetically averaged mean-field always finds a groundstate with a finite condensate order parameter and lower free energy [24,25]. In the context of SFT, the mean-field approximation can be understood as neglecting the kinetic contributions of uncondensed bosons in the self-energy functional [48]. By including these contributions in our SFA3 calculations, we are able to change the energy balance with respect to the mean-field approach, yielding a phase transition to the uncondensed Bose glass.

Strongly-interacting Bose glass phase
Using the local occupation probabilities r n of Sec. 3.2.1, it is possible to distinguish different regimes of the Bose glass in the atomic limit: the qualitative behavior of the Bose glass will change every time the disorder strength is large enough to activate a given local occupation n (i.e. if the probability of finding a site with local occupation n, r n , becomes non-zero as a function of D). Coming from the Mott-insulating groundstate at density n = 1 (where r n =1 = 0) and increasing the disorder strength D, as we enter the Bose glass one of the probabilities r 0 and r 2 becomes non-zero, as either empty or doubly-occupied sites are activated by the disorder depending on the value of the chemical potential. When the disorder is increased further, also higher occupancies are activated and other probabilities r n become non-zero every time we enter a new regime of the strongly-interacting Bose glass.
While the atomic limit shows sharp transitions between the different regimes (see the values of r n in Fig. 2d), for finite hopping J, the kinetic fluctuations turn the transitions into crossovers. However, as we will discuss in this section, in the case of strong interactions the qualitative behavior of local observables changes drastically also in our numerical results whenever a new regime is entered. For the sweep in disorder strength of Fig. 2 the results for local quantities such as the density (Fig. 2a) and the interaction energy (Fig. 2c) show perfect agreement between the analytic results in the atomic limit and both SFA3 and QMC results, except right at the transition/crossover between the different regimes. In fact, the kinetic energy (Fig. 2b) -which is the dominating additional contribution of the finite hopping in SFA3 and QMC, as compared to the atomic limit -is orders of magnitude smaller than the interaction energy at large disorder. In the following we will discuss these different strongly-interacting regimes in more detail, analyzing the qualitative behavior of the observables in Fig. 2 and extracting additional information from the corresponding local spectral functions A loc (ω) shown in Fig. 3.
We start at D = 0, i.e., in the non-disordered Mott insulator. As every site has the same local occupation n i = 1, the local spectral function (Fig. 3b) is characterized by the two Hubbard bands corresponding to the transitions (1 → 0) at negative frequencies and (1 → 2) at positive frequencies. While in the atomic limit these resonances would

HaL
HbL DêJ HcL correspond to delta-peaks, at finite hopping the shape of the spectral function depends on the non-interacting dispersion and its bandwidth W = 2zJ, where z = 6 is the coordination number of the lattice. In particular, the unit filling Mott insulator lower and upper Hubbard bands have the bandwidths W and 2W respectively, see Ref. [59] for a derivation. For weak disorder D < W the qualitative behavior remains the same. However, the Hubbard bands are broadened by the finite disorder strength D and the spectral weight at the center of the bands is reduced, see Fig. 3b. The situation changes when D > W (see Fig. 3c), where the spectral function is more similar to the one predicted by the atomic limit: as discussed in Sec. 3.2.2, the width of the Hubbard bands now is fully determined by the disorder strength D, and the dispersive features of the spectral function are lost. As we are still in the Mott phase, the spectral function shows a finite gap, defined as the minimal distance between the Hubbard bands and ω = 0. As the disorder strength D is increased, so does the kinetic energy (see Fig. 2b), due to increasing kinetic fluctuations, while the gap of the spectral function decreases (see Fig. 3a).
At D ≈ µ (D/J ≈ 56) the gap goes to zero, and we enter the Bose glass phase. The disorder activates empty sites (i.e. r 0 > 0, see Fig. 2d), and as a consequence the density drops (Fig. 2a), while the kinetic energy decreases (Fig. 2b). The lower Hubbard band now extends to ω = 0, and we find a finite spectral weight at small positive frequencies corresponding to the local excitation (0 → 1) of the unoccupied sites (Fig. 3d). In order to study trends in the spectral weight at zero frequency ω = 0, we introduce the spectral weight measure where δ = 0.002U . As shown in Fig. 3a, in this first regime of the Bose glass, the spectral weight ρ 0 for finite hopping is very close to the atomic limit result. In fact, the spectral function (Fig. 3d) only differs from the atomic limit result at the edges of the upper Hubbard band, corresponding to the excitation (1 → 2), indicating a strong localization around empty sites with large values of η i . The situation changes abruptly for D U − µ (D/J 84). As doubly occupied sites are activated by the disorder (see Fig. 2d), the density increases (Fig. 2a), and so does the kinetic energy (Fig. 2b), indicating an increase of non-local kinetic processes. The additional doublons lead to a substantial increase in interaction energy (Fig. 2c), which dominates over the kinetic energy. One would therefore naively expect a better agreement between the spectral functions computed with SFA3 and in the atomic limit. This is however not the case for the spectral weight around zero frequency ρ 0 which increases abruptly at D/J ≈ 84, see Fig. 3a, deviating markedly from the atomic limit prediction. The appearance of doubly-occupied sites in the atomic limit drives additional excitations (2 → 1) and (2 → 3) in the spectral function (Fig. 3e), which overlap with other excitations, leading to additional "bands" composed of multiple resonating excitations, see e.g. (1 → 0, 2 → 1) at low negative frequencies in Fig. 3e.
It is at the edges of these new "bands" that the spectral function is strongly peaked showing a considerable difference with respect to the atomic-limit spectral function, indicating delocalization of quasi-particles and quasi-holes in the vicinity of the rare sites with occupation n > 1 (i.e. occupation 2 in the atomic limit). However, in the Bose glass discussed here, there is no global superfluid response, as the sites contributing to these peaks are rare. Instead the physics is described by the notion of isolated superfluid lakes [21] around rare sites with particularly low local potential.
In this regime (denoted by the grey area in Fig. 2), the non-local Green's function of SFA3 develops a simple pole at zero Matsubara frequency, which can be integrated out when computing local quantities such as the local Green's function, see Appendix F for details. Whence, the self-energy functional and local observables can still be evaluated in this regime. In a homogeneous system, such a pole would indicate an instability towards spontaneous U (1)-symmetry-breaking and the particles would condense. Here, however, it is only the rare sites with n > 1 that contribute to the pole, not allowing for a global condensate. The pole therefore implies the presence of isolated quasicondensates on the lattice, which can have different U (1) phases and therefore do not allow for global phase-coherence (i.e. a finite superfluid response). These highly nonlocal processes in the vicinity of a superfluid phase transition cannot be expected to be fully captured by the self-energies of a local reference system with translationally invariant variational parameters, leading to a deviation in the SFA3 kinetic energy with respect to the numerically exact QMC data in Fig. 2b. This was also the case in close proximity to phase transitions in SFT [48] and BDMFT [33,34] calculations in the clean BHm.
For even stronger disorder, the situation changes when the number of doubly occupied sites in the atomic limit (proportional to r 2 in Fig. 2d) Fig. 2c) makes it harder for particles to delocalize. This can be observed in the kinetic energy of Fig. 2b, which decreases again as the particles localize. The same behavior can also be seen in the spectral function of Fig. 3f, where the bands involving highly occupied sites increase in width, but are much closer to the atomic limit results. The zero frequency spectral weight ρ 0 decreases accordingly, as shown in Fig.  3a.

of interaction energy in
When the disorder is strong enough to activate triplon occupancies, the behaviour changes once more. The kinetic energy (Fig. 2b) increases as the particles delocalize around the rare triply occupied sites, and so does the interaction energy (Fig. 2c). The number of doubly occupied sites on the other hand decreases and r 2 = r 1 (see Fig. 2d). This behavior arises naturally from the probabilities r n of Sec. 3.2.1: once D > η max (n) and D > −η min (n), the probability of finding a site with local occupation n > 0 is given by the particle-number-independent value r n = U/2D (as in this case r 1 and r 2 ).
In summary, our results show that at fixed interaction U/J (and chemical potential µ/U ) the strongly interacting Bose glass as a function of D is described by the subsequent activation of local occupations n. As these occupations accumulate, the interaction energy increases, driving the phase towards the atomic limit. This is however not the case whenever a particular higher occupation number is very rare (i.e. if 0 < r n 1 for some local occupation n > 0): in this case the particles tend to delocalize and form superfluid lakes [21] around these rare highly occupied sites.

Strongly-interacting phase transition
The regime where the Bose glass exhibits superfluid lakes around doubly occupied sites surrounds the superfluid finger at large interactions in the phase diagram of Fig. 1. In fact, the lower edge of the superfluid finger strongly correlates with the line where doubly-occupied sites are activated in the atomic limit (see blue dashed line in Fig. 1). The strongly-interacting phase transition at fixed chemical potential is illustrated in Fig. 4, where we show the superfluid to Bose glass phase transition as a function of U/J at fixed chemical potential µ/U = 0.25 and disorder strength D/J = 100.
With decreasing interaction U/J, these superfluid lakes percolate and resonances between the low-energy excitations (0 → 1) and (1 → 2) on neighboring sites (and between the corresponding particle-removal processes, see e.g. Fig. 3e) favour the spontaneous breaking of U (1)-symmetry through a homogeneous condensate. The particles therefore eventually condense, driving the transition to the superfluid phase. As the sites contributing to the resonating low-energy excitations remain relatively rare, the condensate fraction and the correction of the density with respect to the atomic limit close to the phase transition are extremely low (with a condensate density on the order of 10 −4 , see inset of Fig. 4a).
At density n ≈ 1 and larger disorder strength, the increase of highly occupied sites is compensated by a proliferation of empty sites (see e.g. Fig. 2d). Thus, the probability of having neighboring sites with resonating low-energy excitations decreases, making the spontaneous breaking of U (1)-symmetry less likely. The increased interaction energy and particle number fluctuations therefore suppress the superfluid phase, explaining the reentrant behavior of the superfluid finger at larger disorder (see Fig. 1).

Superfluid phase
We now turn to lower interactions, i.e., deeper into the superfluid phase away from the superfluid finger. If U/J is lower than the critical value of the clean system, the condensate density is much larger than in the superfluid finger, and the uncondensed background is no longer well described by the atomic limit.
In Fig. 5 we show a sweep of the thermodynamical observables as a function of D/J deep in the superfluid phase at U/J = 20 and µ/U = 0.35. At low disorder SFA3 shows excellent agreement with QMC, as the condensate density increases and the density decreases as a function of disorder (Fig. 5a). As a consequence of the larger condensate fraction, the magnitude of the kinetic energy increases as well (Fig. 5b). The interaction energy increases throughout the entire parameter range 0 ≤ D/J ≤ 25, indicating increasing spatial particle number fluctuations (Fig. 5c).
When the disorder becomes comparable to the single-particle bandwidth W = 2zJ, these fluctuations reverse the trend of the condensate density which starts to decrease as a function of D/J. It is at this point that also the kinetic energy starts to decrease and the non-local connected Green's function develops a pole at zero Matsubara frequency (see Appendix F). As in the Bose glass (see Sec. 4.1), the pole indicates the appearance of additional isolated quasi-condensates in the system: this is most likely related to the disorder inducing rare regions, explaining the decrease in condensate density and kinetic energy, and leading to a glassy behavior in the superfluid.

HaL
HbL HcL Eventually, deeper in the glassy regime of the superfluid where the disorder dominates over both the interaction and the single-particle bandwidth, our SFA3 approach of having translationally-invariant variational parameters on the reference system becomes too simple to fully capture the groundstate behavior. In fact, the SFA3 results start to deviate from the QMC results, see Fig. 5a. As shown in Fig. 5d, eventually at D/J ≈ 2W the variational parameters of the stationary SFA3 solution join with a metastable solution with higher free energy through a saddle-node bifurcation [60], vanishing for larger disorder.
A possibility to get around this problem, may be the introduction of a spatially modulated symmetry-breaking field on the reference system, as was also done in stochastic mean-field theory [26,27]. In SFT, this would however involve the inversion of a non-translationally-invariant connected Green's function, limiting us to very small system sizes, while we here want to analyze the thermodynamic limit.

Conclusion
In this work, we generalized the bosonic self-energy functional theory (SFT) to include quenched disorder. The derived formalism is a general framework for constructing non-perturbative approximations of disordered interacting bosonic lattice systems incorporating spontaneous U (1)-symmetry-breaking. We showed that the resulting SFT functional depends only on the self-energies of the disorder-averaged interacting one-and two-point propagator, the condensate and connected Green's function, respectively. The lattice self-energies can then be parametrized by the self-energies of a simpler exactly solvable reference system having the same interaction and disorder distribution as the original system. The resulting formalism is a general non-perturbative approach that contains disorder-averaged bosonic dynamical mean-field theory as a certain limit.
We applied SFT in combination with a simple SFA3 reference system, consisting of a single bosonic mode with only three variational parameters, to the Bose-Hubbard model with local box disorder. The SFA3 results were compared to numerically exact path integral quantum Monte Carlo (QMC) results and analytic calculations in the atomic limit.
Our results in the strongly-interacting regime close to unit filling, showing excellent agreement with QMC, indicate that the Bose glass phase is characterized by different regimes as a function of the disorder strength D. With increasing D, sites with local occupations n = 1 appear as predicted by the atomic limit, leading to crossovers between different regimes whenever a new local occupation n is activated by the disorder. While QMC has to resort to analytic continuation in order to compute dynamical quantities, through SFT we were able to compute spectral functions within SFA3 directly.
By systematically analyzing the local spectral function we observed that the bosons delocalize into superfluid lakes around highly occupied sites whenever these are particularly rare. In particular, we found that the transition from the stronglyinteracting Bose glass to the strongly-interacting superfluid phase (which extends to values of the interaction which are much larger than in the clean system) is driven by the percolation of superfluid lakes which form around doubly occupied sites, leading to a small condensate fraction over a strongly localized background. As D is further increased and the density of doublons increases accordingly, the particles are localized by the strongly interacting background, explaining the reentrant behavior of the superfluid phase.
Due to the extremely low condensate fraction in the strongly-interacting superfluid, even though the numerical error is small, the phase boundaries observed with SFA3 are shifted with respect to the QMC results. Deeper in the superfluid phase (i.e. at lower interactions), our SFT results show excellent agreement with the QMC data as long as the disorder is smaller or comparable to the non-interacting bandwidth. In the stronglydisordered weakly-interacting regime, however, the restricted variational subspace of the SFA3 reference system employed in this work is no-longer capable to find a stationary solution.
As opposed to QMC, SFT does not suffer from a general sign problem in the presence of e.g. gauge fields [61,62,63], or other complex Hamiltonian terms such as spin-orbit coupling [64,65,66]. The formalism derived in this work therefore represents a promising tool for future studies of such complex systems in the presence of disorder. In particular, an extension to real-time dynamics, as has been done for disorder-free fermionic systems [67,68], seems a promising route to study the elusive physics of many-body-localized systems and is left for future work.
In order to construct approximations using the disorder averaged propagators we now seek a functional that is equal to the disorder average of the self energy functional but that is defined in the extended space of both averaged and explicit self-energiesS,Σ, {S η , Σ η } and stationary at the physical solution in all self energies. An ansatz that fulfills equality at the physical self-energies is Eq. (32). Repeated application of the Dyson equations [Eqs. (27) and (14)] in Eq. (33) at the physical self-energies giveŝ and thereforeΓ whence the disorder averaged self-energy functional gives the physical disorder averaged free energy at stationarity. To show stationarity ofΓ FtP V we consider the variations of the universal functionalT P V

Appendix B. Canceling functional derivatives ofT P V
At stationarity the expression of theT P V functional in Eq. (A.1) contains expressions in terms of G tη0 and F, so one might wonder if there is no implicit dependence on the free propagators of the system. In order to check that this is not the case, we rewrite Eq. (A.1) in terms of G t00 and F aŝ (B.1) Self-energy functional theory with symmetry breaking for disordered lattice bosons 27 The variation in G −1 t00 yields which by the short-hand notations introduced in Eqs. (31), and (31) can be rewritten as where we have used that the trace and the arithmetic average commute, i.e. 1 2 TrG η P = 1 2 Tr G η P . Note that -as opposed to the arithmetical average -the geometrical average used in the context of fermionic DMFT [36,37,38,39] would not commute with the trace operator Tr in Eq. (B.3) and therefore break the universality of the functionalT P V . As pointed out in Ref. [51] for fermions, the geometrical average introduced in DMFT, therefore appears to be incompatible with SFT.
The variation ofT P V in F yields which using Eq. (31) can be rewritten as T P V is therefore completely independent of both F and G −1 t00 also at stationarity.
Appendix C. Disorder-averaged bosonic dynamical mean-field theory limit Disordered-averaged SFT for bosons has disorder-averaged bosonic dynamical meanfield theory (BDMFT) as a certain limit. In its simplest form, disorder-averaged BDMFT is restricted to site-local disorder η and site-local interactionV .
In this case, disorder-averaged BDMFT is obtained from SFT by restricting the reference systems free propagator to be site-local, i.e.
where i, j are the site-, and ν, ν the Nambu indices. The imaginary time retardation in ∆(τ, τ ), however, remains completely general. The reference systems local bare propagator G ∆η0 and interaction give rise to a purely local self-energy FtP V [Eqs. (45) and (46)] reduce to the disorderaveraged BDMFT self-consistency conditions which can be fulfilled exactly by the retarded ∆(τ, τ ), and can be simplified tō whereΦ andḠ ii are the disorder-averaged condensate and local connected Green's function of the lattice, whileΦ F ∆P V andḠ F ∆P V are the disorder-averaged condensate and connected Green's function of the reference system. This is therefore the standard BDMFT self-consistency condition of clean systems [33,34], where the propagators of the clean system have been replaced by their disorder-averaged counterparts, which for the case of uncorrelated disorder discussed in Appendix D are translationally invariant.
Appendix D. Uncorrelated disorder: translational invariance of the arithmetic average In the following we will specialize the formalism derived in Sec. 2 by assuming that the disorder is distributed according to an uncorrelated and translationally invariant probability distribution, i.e. where the product goes over all site-indices i, j, and the distribution p i−j depends only on the relative distance between the sites i and j. As we will see this enables us to simplify the reference system considerably due to the translational invariance of disorder-averaged observables.
The interacting propagators at a given disorder configuration η can be computed directly by where T is the time-ordering operator and . . . η means taking the expectation value with respect to the reference system with disorder configuration η.
Using Eq. (D.2) further enables the computation of the fixed-disorder self-energies through Eq. (14). The propagators G F ∆ηV , Φ F ∆ηV , and the corresponding self-energies Σ F ∆ηV , S F ∆ηV , are not translationally invariant and can therefore be very hard to handle numerically.
If we now assume that we average over an infinite number of disorder configurations, the reference system's propagators will be translationally invariant, since due to the translational invariance of the uncorrelated disorder probability distribution of Eq. (D.1) all values η ij will occur with the same weights for each pair of sites (i, j) with the same distance i − j, i.e.
with a translationally invariant condensatē According to Eq. (14) this implies that also the average self-energies will be translationally invariant with Finally, Ω F ∆P V = Ω F ∆ηV P can be computed directly from averaging over the fixeddisorder systems. As no fixed-disorder quantities are needed in order to evaluate the functional in Eq. (47), the evaluation of the self-energy functional has now the same complexity as the disorder-free case of Ref. [48], where the self-energies and propagators were translationally invariant by definition. The only difference lies in the treatment of the reference system, which has to be averaged over all disorder configurations η.

Appendix E. Lattice observables
Once a stationary solution fulfilling Eq. (52) has been found, the corresponding lattice observables can be computed using the self-energies Σ η ≈ Σ F ∆ηV , S η ≈ S F ∆ηV ,Σ ≈Σ F ∆ ,S ≈S F ∆ . (E.1) In particular, the disorder-averaged propagatorsḠ andΦ can be computed using the self-energiesΣ andS and the free propagator G t00 in the Dyson equations of Eq. (14). The fix-disorder propagators G η and Φ η , on the other hand, can be computed using Σ η and S η and the free propagator G tη0 in Eq. (14). As the latter are not translational invariant, however, they can only be computed on a finite sized lattice, as Eq. (14) requires the inversion of a matrix in position space. It is therefore preferable to use the translationally invariant averaged propagatorsḠ andΦ in the thermodynamic limit. As the arithmetic averaging is a linear operation, disorder-averaged observables of the lattice system which can be expressed as linear terms of one-and two-point quantities without any disorder-dependent prefactors, can be directly evaluated from the average propagatorsΦ andḠ. This is trivially the case for the disorder-averaged condensate through Eq. (31), while for the particle density we have n = 1 2βL Tr Tr[−Ḡ] +Φ †Φ , (E.2) where we have used Eq. (31) in the last step and the same definition of the trace operator Tr as in Ref. [48]. The same is true for the kinetic energy Tr t Ḡ −Φ †Φ .

(E.3)
The interaction energy, on the other hand, cannot be directly evaluated from the averaged propagators as [48] Tr[Σ η G η ] P = − 1 4βL Tr [ΣḠ], (E.4) However, as the SFT functional is equal to the free-energy at stationarity, we have direct access to the disorder-averaged free energy of the lattice Ω FtP V , from which we can compute the interaction energy by the numerical derivative boundary, see the grey regions in Figs. 2 and 5. In the clean system, such a pole signals an instability towards U (1)-symmetry breaking and arises only in the metastable Mott insulator phase. In the case studied here, which is no-longer homogeneous, as discussed in Sec. 4, the pole is related to the appearance of isolated quasi-condensates on the lattice.
Remarkably, although the poles make non-local quantities such as, e.g., n k = − n TrḠ(iω n , k)/2β diverge at certain values of k, the pole can be treated semianalytically in the computation of local quantities, as we will show in the following.
The central quantity where the lattice Green's function enters in the SFT functional of Eq. (47) is the trace-log term Tr ln −Ḡ −1 , which -as shown in Ref. [48] -is only defined up to a regularization factor C ∞ and can be evaluated as where R is the regularization function where D( ) is the single-particle density of states, and Q(iω n ) is the reguarlization function Q(iω n ) = ω 4 n , n = 0, β −4 , n = 0.
(F. 8) In order to evaluate the integral in Eq. (F.7) numerically, the dispersion is discretized on the energy grid = m . Using the linear interpolatioñ where the integrand is given by If the interval [ m , m+1 ] does not contain the poles p ± (iω n ), the mth summand of Eq. (F.10) can be straight-forwardly integrated analytically. Also in the presence of a