{'fraction_non_alphanumeric': 0.03482429559947538, 'fraction_numerical': 0.009248790194925603, 'mean_word_length': 4.432800982800983, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variational cluster approach for strongly correlated lattice bosons in the superfluid phase

We extend the variational cluster approach to deal with strongly correlated lattice bosons in the superfluid phase. To this end, we reformulate the approach within a pseudoparticle formalism, whereby cluster excitations are described by particlelike excitations. The approximation amounts to solving a multicomponent noninteracting bosonic system by means of a multimode Bogoliubov approximation. A source-and-drain term is introduced in order to break U(1) symmetry at the cluster level. We provide an expression for the grand potential, the single-particle normal and anomalous Green's functions, the condensate density, and other static quantities. As a first nontrivial application of the method we choose the two-dimensional Bose-Hubbard model and evaluate results in both the Mott and the superfluid phases. Our results show an excellent agreement with quantum Monte Carlo calculations.


I. INTRODUCTION
Cluster approaches have been proven to be very useful for the numerical investigation of strongly correlated many-body systems. These approaches consist in embedding finite size clusters, for which a numerical exact solution is available, within a lattice of infinite size. The embedding is done by introducing additional fields to the cluster Hamiltonian, in order to take into account the coupling to the rest of the lattice in some appropriate dynamical mean-field way. We will term these fields Weiss fields, since they play an analogous role as in Weiss meanfield theory of ferromagnetism (cf. Ref. 1). Different cluster embedding techniques, such as cluster perturbation theory 2,3 (CPT), variational cluster approach 4-6 (VCA) , cellular dynamical mean field theory 7 (C-DMFT), and dynamical cluster approximation, 8 differ by the nature of the Weiss fields and of the mean-field treatment which fixes their optimal value. In the present paper, we consider VCA, which has been applied to a large variety of fermionic 4,6 and bosonic [9][10][11] systems. VCA can be understood in a more general framework called self-energy functional approach 5,12 (SFA), in which the grand potential of the physical system is expressed as the stationary point of a particular functional of the self energy. Here, we will adopt an alternative approach to VCA in which single-particle excitations are expressed in terms of "pseudoparticles," which are similar to Hubbard operators, 13 and external fields are "added" to the cluster Hamiltonian and "subtracted perturbatively." 14 We discuss in Sec. II the advantages of this alternative approach.
Strongly correlated lattice bosons are currently in the focus of research due to seminal experiments on ultracold gases of atoms. [15][16][17] In these experiments quantum mechanical interference effects can be observed on a macroscopic scale. In particular, ultracold gases of atoms on a lattice undergo a quantum phase transition from the Mott phase, in which particles are localized on individual lattice sites, to the delocalized superfluid phase, in which U (1) symmetry is broken and a finite fraction of the particles forms a Bose-Einstein condensate (BEC).
Up to now, bosonic VCA has been formulated for the normal phase only. The principal aim of this paper is to extend this formalism to the symmetry-broken, superfluid phase. The theoretical framework, developed in the next two sections, is applicable to a large class of lattice boson systems in the Mott insulating as well as in the superfluid phase. In particular, besides the widely studied Bose-Hubbard model, 15,18 the method can be straightforwardly extended to include disordered systems or multiple components containing, for example, fermionboson mixtures. The extended VCA theory can be applied even to the U (1) broken, superfluid phase of lightmatter systems, where photons are confined in coupled, nonlinear quantum-electrodynamics cavities. 19,20 In order to achieve the extension to the U (1) broken, superfluid phase, it proves convenient to reformulate VCA in terms of a pseudoparticle approach, whereby singleparticle excitations within a cluster are approximately mapped onto particlelike excitations. We show that this approach, first applied to normal bosons, quite naturally suggests the extension to the superfluid case. In a following publication, 21 we show that the results obtained from the pseudoparticle formalism in the superfluid phase can be equivalently obtained within an appropriate extension of the SFA taking into account condensed bosons. One of the aims of the present paper is to illustrate the advantages of the pseudoparticle formalism, which can be used to extend VCA to a large variety of problems with strongly correlated lattice systems.
The pseudoparticle formalism is in some aspects related to the standard basis matrix operator method developed by Haley and Erdös in Ref. 22 and to the Hubbard-operator approach, see for instance Ref. 13. The idea is to introduce pseudoparticle operators b µ and b † µ , which connect the ground state |ψ 0 with singleparticle excited states |ψ ν of a Hamiltonian describing disconnected clusters in the lattice. In the VCA language the cluster Hamiltonian is termed reference systemĤ ′ . Of course, the b µ and b † µ do not have the properties of ordinary single-particle creation and annihilation operators. The crucial point is that by treating them as such, one recovers the very same results as obtained from CPT and from VCA, as has been shown for fermionic systems in Ref. 14 (see appendix therein). In Sec. II we prove the same result for the bosonic (normal) case, which is somewhat more subtle, as it requires a multimode Bogoliubov transformation. In this picture, excited states |ψ µ are treated as pseudoparticle excitations with the properties Within the VCA approximation, pseudoparticles are regarded as noninteracting particles. We stress that, while this may seem a rather crude approximation, it is equivalent to CPT and VCA. Furthermore, with appropriate extensions it becomes equivalent to C-DMFT. It is straightforward to show (see Sec. II) that within this approach the original bosonic operators a i and a † i can be expressed as linear combinations of the pseudoparticle operators b µ . This makes it possible to write the coupling of the cluster to the rest of the lattice, which in VCA consists of intercluster hopping terms, as a quadratic form in the b µ . In combination with the fact that the cluster Hamiltonian is by construction quadratic in these operators as well, one finally obtains a Hamiltonian which is completely quadratic in the pseudoparticle operators, and can, thus, be solved exactly.
Our paper is organized as follows: In Sec. II we first show that the standard VCA results for the Green's functions and for the grand potential Ω are recovered within the pseudoparticle formalism applied to normal bosonic systems. In order to be able to treat the superfluid phase we then extend the theory in Sec. III by introducing Weiss fields in the form of "source-and-drain" terms which explicitly break the U (1) symmetry of the reference HamiltonianĤ ′ . The main result of Sec. III is the expression for the grand potential N c Ω of the physical system (N c is the number of clusters). Within our extension of VCA to the superfluid phase we obtain 23 The first three terms on the right hand side are essentially identical to those which are also present in standard VCA expressions. Particularly, Ω ′ is the grand potential (per cluster) of the reference system, and G ′ and G are the connected Green's functions of the reference and of the physical system, respectively. However, they are expressed in the Nambu representation, which explains the additional factor 1/2 and the fourth term in comparison with previous results. [9][10][11] The suffix (0) used in the second line of Eq. (1) means that the corresponding Green's functions are calculated for q = 0 and ω = 0, where q is the superlattice vector associated to the cluster tiling, and ω is the Matsubara frequency. As usual within VCA theory, the two Green's functions share the same selfenergy. The expectation values A ′ and A are the corresponding condensate densities, again in Nambu (vector) notation. The latter are connected by the relation where the vector F describes the strength of the sourceand-drain term which is introduced in the reference system in order to explicitly break U (1) symmetry. The value of F [see Eq. (23)] has to be determined from the variational principle. Details for the notation are provided in Sec. II, and III. In addition to the formula for the grand potential Ω, we evaluate expressions for other quantities, which are useful for describing the superfluid phase. In particular, we derive expressions for the normal and anomalous Green's functions, the particle density, and the condensate density. In Sec. IV this extended VCA theory is applied to the two-dimensional Bose-Hubbard (BH) model in the superfluid phase. Finally, we summarize and conclude our findings in Sec. V.

II. PSEUDOPARTICLE APPROACH
In this section we reformulate CPT/VCA within the pseudoparticle approach for bosonic systems. In principle, one may argue that the formulation of CPT/VCA using pseudoparticles is complicated and in the case of the normal phase (i. e. Mott phase) CPT/VCA can be obtained from simpler approaches, as, for example, from Dyson's equation (see, e.g. Ref. 3) and the SFA. 5,9 The reason why we present this alternative formulation here is that this approach, while not as rigorous as SFA, provides useful hints on how to deal with more complicated situations, like the superfluid phase discussed in this work (see Sec. III). In addition, it gives insight on other properties. For example, in the case of normal bosons the pseudoparticle approach is useful in order to understand the occurrence of noncausality of the Green's function in cases, where the chosen reference system is not suitable to describe the phase of the physical system, as we point out below. Thus the aim of this section is to derive the principal theoretical framework of the pseudoparticle approach, for the normal phase, reproducing the known result for the grand potential Ω, which has to be optimized. The extension to the superfluid phase is the subject of the next section.
The physical system of interacting particles is described by a grand-canonical HamiltonianĤ, which is related to the canonical Hamiltonian in the usual way by the additional single-particle term −µN . The Hamiltonians, which can be treated by the extended VCA theory, generally have the formĤ =Ĥ t +Ĥ U , whereĤ t consists of arbitrary one-particle terms andĤ U of local twoparticle terms. The physical system is defined on a large or even infinite lattice with periodic boundary conditions. The underlying lattice is now tiled into N c clusters each one containing L orbitals (sites). We split the Hamiltonian into a cluster partĤ cl , which only describes processes within the various clusters, and the residual part T , containing the intercluster processes, which consist of single-particle terms only, so that CPT amounts to first solving for the HamiltonianĤ cl and then carrying out a perturbation expansion in the intercluster HamiltonianT . Of course, within CPT one is free to add an arbitrary single-particle Hamiltonian −∆ to the cluster HamiltonianĤ cl provided it is then subtracted fromT so thatĤ remains unchanged. This defines a new cluster HamiltonianĤ ′ The physical HamiltonianĤ, given in Eq. (3), can now be expressed in terms of the new cluster Hamiltonian leading to a new "perturbation"T ≡∆ +T . The CPT expansion is now carried out in this new "perturbation". While ideal exact results should not depend on∆ (this occurs, for example, in the noninteracting case), in practice results do depend on∆ due to the approximate nature of the expansion. The idea is to fix the parame-ters∆ by an optimization prescription, which amounts to finding the stationary point of the grand potential Ω obtained from the perturbative expansion. The optimization prescription is put on a rigorous framework within the SFA. 12 It is straightforward to show that this procedure is equivalent to the standard VCA prescription, wherebyĤ ′ is the corresponding reference system. 24 In the following, we consider N c identical disconnected clusters, and denote the sites (orbitals) within a cluster by i. The position of each cluster on the large, physical lattice is specified by a lattice vector R. Accordingly, we denote by a i,R the annihilation operator for a boson on site i of cluster R, and similarly for creation operators a † i,R . In order to keep a compact notation we combine the annihilation operators of a given cluster R into a column vector of operators a R = (a 1,R , a 2,R , . . . a L,R ) T , and correspondingly, the creation operators are row vectors a † R = (a R ) † . Using these expressions we rewrite the intercluster Hamiltonian aŝ where t(R−R ′ ) is a matrix describing the hopping terms from cluster R ′ to cluster R, with the property t(R − R ′ ) = t(R ′ − R) † . Here we have assumed translation invariance by a cluster translation vector. Similarly, we can express∆ in terms of an intracluster hopping matrix h∆ such thatT , defined in Eq. (5), can be written as Eq. (6) with the replacement As explained above, the reference system consists of a sum of Hamiltonians acting on independent clusters R Again considering translation invariance, allĤ ′ (R) are identical. Thus it suffices to determine numerically the ground state |ψ 0 , R , as well as single particle or singlehole excited states |ψ µ , R of a single cluster Hamiltonian H ′ (R), with corresponding eigenenergies E ′ 0 and E ′ µ , respectively. The key idea of the approach, to be presented here, is to introduce pseudoparticle operators b † µ,R and b µ,R ≡ (b † µ,R ) † , which are defined by their matrix elements In other words, the pseudoparticle operator b † µ,R applied to the exact many-body groundstate |ψ 0 , R of a cluster creates the exact excited many-body state |ψ µ , R . In this respect, it is of course forbidden to apply a second pseudoparticle creation operator on the excited state. This leads to the supplementary hard-core con- To neglect this hard-core constraint and to restrict to single-particle and singlehole excitations within each cluster is the approximation made here. We show below that this approximation, combined with the variational procedure discussed above, gives the same results as VCA. In particular, we obtain the same expression for the grand potential Ω, and for the Green's function. It should be mentioned, however, that within the pseudoparticle approach there is no known rigorous variational principle for Ω. One can simply heuristically state that the "best" solution is the one that "minimizes" the energy, although, as we know from VCA, the variational solution is not always a minimum. Also for parameters, such as the chemical potential, for which Ω turns out to be a maximum, one can argue that the stationary condition is a kind of "constraint" fixing the consistency of thermodynamic quantities, 25 and the corresponding parameter is a kind of "Lagrange multiplier." Nevertheless, it is not the goal of the present paper to discuss this issue. Here, we want simply use this "tool" in order to formulate an extension of the theory to address the bosonic superfluid phase (see Sec. III). The knowledge of the correction to the order parameter and of the grand-potential Ω can then guide and facilitate a rigorous extension of SFA to deal with the superfluid phase. This is the goal of a future publication. 21 With the help of these operators, it is straightforward to write down a Hamiltonian which has the same energies and eigenvectors as the reference system, restricted within the subspace of single-particle and single-hole excitations from the ground statê with the (positive) excitation energies Since we are interested in zero temperature T = 0, the grand potential of the reference system is Ω ′ ≡ E ′ 0 . To proceed further, we need an expression forT , and, thus, of the original bosonic operators a i,R , in terms of the pseudoparticle operators. For simplicity, we drop the R index and concentrate on a given cluster. Within the pseudoparticle approximation the operators must coincide only within the constrained subspace. We thus approximate each a i by an operatorÔ i (b µ , b † µ ) which shares the same matrix elements ψ 0 | · |ψ 0 , ψ 0 | · |ψ ν , and ψ ν | · |ψ 0 . We expressÔ i by means of the ansatẑ where the first sum contains the n p indices associated with the single-particle excitations, and the second sum contains the n h indices for the single-hole excitations. The total number of excitations taken into account is n s = n p + n h . Here we have exploited particle-number conservation. Next, we use this expression to evaluate the following matrix elements where the coefficients γ i are zero so far, since the reference system conserves the particle number. We now introduce the compact notation i. e., the first part of the vector acts on particle states, and the second part on hole states. Notice that in this form B † (B) changes the number of particles by +1 (−1). We also introduce the Q matrix (which is the same as in Ref. 11) as The Q matrix can be used to express the original operators a and a † in terms of B operators [cf. Eq. (9)] in a compact form: Using the compact vector notation for B and B † , the reference Hamiltonian [Eq. (8)] can be written aŝ where we reintroduced the R dependence. Here we introduced the diagonal matrices Notice that S 2 = 1, while Λ contains the poles of the Green's function for the reference system. The constant with the function takes into account that some of the boson operators have been rearranged in order to obtain Eq. (12). The physical Hamiltonian introduced in Eq. (5) readŝ Using Eqs. (11) and (12) yields a quadratic expression in the B operators: We can now introduce a Fourier transform in the cluster vectors R leading toĤ Here, we have introduced the matrix is the Fourier transform oft(R−R ′ ). The non-Hermitian matrix M q is identically defined as in Ref. 11. Being quadratic in the B operators, H q can be quite generally put into diagonal form by a multimode Bogoliubov transformation. To achieve this, we look for "normal-mode" pseudoparticles described by the vector P with the same structure as B (in the following considerations we omit the q dependence for simplicity) where V is a nonsingular but in general nonunitary matrix. From a physical viewpoint the nonsingularity of V corresponds to a pseudoparticle conservation, meaning that there are as many pseudoparticles B as normalmode pseudoparticles P . The transformation V must satisfy two conditions. First it must be chosen such that P has appropriate bosonic commutation relations, i. e., This gives The second requirement on V is since after the transformation from B particles to P particles the Hamiltonian in Eq. (15) has to be diagonal. Multiplying Eq. (17) from the left by V S ′ and using Eq. (16a) yields the eigenvalue equation where D ≡ diag(d 1 , . . . , d ns ) = S ′ E contains the eigenvalues of the non-Hermitian matrix M . From Eq. (18) below, where we express the Hamiltonian in terms of the normal-mode pseudoparticles, it can be seen that the diagonal elements e i correspond to the excitation energies of the physical system. Since the energy of the physical system must be bounded from below, all e i have to be positive and real, leading to It will turn out that this stability condition is the only point, where the variables s ′ i of the auxiliary operators p s ′ i i show up. In App. A we show that, if M is completely diagonalizable with real eigenvalues and linear independent eigenvectors, which is of course not generally guaranteed for a non-Hermitian matrix M but necessary from the physical viewpoint, then V can be constructed so that both requirements of Eqs. (16) and (17) are fulfilled, and we can proceed with our analysis. If M is not completely diagonalizable or does not have real eigenvalues, the system is unstable, and it favors a different phase, which cannot be addressed by the reference system in this form. This instability toward a different phase, such as superfluidity, has to be cured by extending the reference system by proper additional variational parameters, as discussed in Sec. III.
In terms of the P operators we obtain for the Hamiltonian In the last line we have exploited the fact that in order for the system to be stable, i. e., the energy be bounded from below, all e ν must be positive. Inserting this expression in Eq. (14) yields the Hamiltonian in terms of diagonal normal modes. From this result one obtains immediately the grand-canonical groundstate energy per cluster As discussed, Λ and D q are diagonal matrices containing the poles of the reference Green's function and physical Green's function, respectively. Therefore, this expression being equivalent to Eq. (11) in Ref. 11 (see also Refs. 9, 12, and 26) is equivalent to the zero-temperature VCA grand potential. By using the expression for the Green's function of the noninteracting normal modes we readily obtained the Green's function for the physical system By simple algebra this expression can be rewritten such that it is independent of the auxiliary quantities S ′ and V and thus equivalent to Eq. (12) in Ref. 11 where we have used Eqs. (16) and (17). We, therefore, succeeded in proving that, for normal bosons, the pseudoparticle approach yields the same Green's function and grand potential as VCA. For fermions, this was shown in Ref. 14, see appendix therein. This result holds for T = 0, although extension to T > 0 is straightforward.

III. SUPERFLUID PHASE
When trying to apply VCA to bosonic lattice systems in regions of the phase diagram outside the Mott phase, one encounters instabilities which manifest in the form of noncausal Green's functions, i. e., in spectral functions with negative (positive) spectral weight for positive (negative) frequencies ω, or in complex poles. Within the pseudoparticle approach these instabilities show up as complex eigenvalues or negative diagonal elements of the matrix E. This kind of instability is well known in approaches based on the bosonic Bogoliubov approximation, such as the spin-wave approximation.
Quite generally, such an instability signals the occurrence of a phase transition toward a new phase. Quite often, as in the case of the BH model studied in Sec. IV, the new phase is the superfluid phase, which is accompanied by a Bose-Einstein condensation. Bose-Einstein condensation is described by a finite value of the order parameter a R . This suggests to include a source-anddrain term in the reference system, which breaks the U (1) symmetry of the reference system, leading to the "perturbation"T [see Eq.
where f R ≡ (f 1 , f 2 . . . f L ) T is a vector of size L and is identical for all clusters. The index R, however, will be kept for notational reasons. Due to these terms, the reference system Hamiltonian does not conserve particle number anymore. Its eigenstates will thus consist of superpositions of states with different particle numbers. Numerically, a cutoff in the maximum number of boson is necessary in order to solve the reference system on the cluster level exactly. We again introduce pseudoparticle operators b R connecting the ground state with excited states. Note that we cannot distinguish between particle or hole states anymore. The pseudoparticles are defined by Eq. (7) and are connected to the original boson operators a R by means of Eq. (9). Now, all matrix elements in Eq. (10) are nonzero in general. Therefore, the two sums over µ in Eq. (9) are extended to µ = 1, . . . , n s , where n s is the number of excited states considered in each cluster.
For the following considerations it is convenient to express the boson operators within a Nambu notation. For the particle operators we introduce in real space which after a Fourier transformation in the cluster vectors, see Eq. (13), becomes For pseudoparticle operators we have in real space Similarly to Sec. II, we have an approximate linear relation between the A operators and the B operators of the form

After the Fourier transformation in the cluster vectors it reads
Here, The constants γ i ≡ ψ 0 | a i |ψ 0 , will be nonzero as the reference system does not conserve the particle number.
In terms of pseudoparticle operators we can again write the reference Hamiltonian for a cluster R, similarly to Eq. (12) asĤ Here, the matrices S and Λ have a slightly different definition ) , To express the "perturbation"T of Eq. (20), we need to introduce a similar Nambu notation for the sourceand-drain terms, which, being R independent, become in q space After the Fourier transformation in the cluster vectors, we can rewritê whereT q = diag(t q ,t T −q ). Replacing the A operators in terms of the B operators with the help of Eq. (21), and combining Eq. (22) with the expression above forT , we finally obtain the complete Hamiltonian, defined in Eq. (5), in terms of pseudoparticleŝ The expression can be further simplified by using the fact that F and Γ are equal in all clusters, and thus have only q = 0 components. In addition we take advantage of since t(R − R ′ = 0) = 0 is a pure intercluster term. For notational convenience we introducẽ This giveŝ The term linear in B can be eliminated by a shift where clearly only the q = 0 term of X q is nonzero. Considering only the q = 0 part of Eq. (26), which we term Y 0 , and plugging in the shifted operators, we obtain The linear term is eliminated by setting yielding for the q = 0 term above which is the Green's function defined in (33) but evaluated for q = 0 and ω = 0. In total we havê with the constant terms In the last term of Eq. (29), we restrict the summation over half of the Brillouin zone, which we denote by q ∈ BZ/2, and thus removed the factor 1/2 in front of the sum. Due to Nambu representation, two summands with +q and −q are identical and therefore the restriction to half of the Brillouin zone is convenient. In our convention, the q = 0 term is included in the sum and retains the factor 1/2.

A. Condensate density
Before turning to the diagonalization of the Hamiltonian in Eq. (29), let us evaluate the condensate density.
Since there are no terms linear inB, its expectation value B vanishes. Therefore, we can immediately calculate the condensate density where we used Eqs. (23), (25), (27) and (28). We now exploit the fact that Γ = A ′ is the condensate density in the reference system. From the Dyson equation for the Green's function of the physical and the reference system we have 28 By multiplying (30) with G −1 (0) we obtain which corresponds to Eq. (2).

B. Diagonalization of the Hamiltonian
The Hamiltonian of Eq. (29) is finally quadratic and its diagonalization proceeds in the same way as in Sec. II. Again we introduce P operators B q = V q P q , and find the solution of the non-Hermitian eigenvalue equation where V q satisfies the relation The diagonal matrix S ′ , which is in principle q-dependent as well, consists of +1 or −1 terms. It is chosen according to the prescription derived in App. A. The stability condition is again that the pseudoparticle eigenenergies are all positive. The physical Hamiltonian in terms of P -particles now readŝ From that we readily obtain (see App. B) the grand potential per cluster of the physical system Ω, which is the ground state expectation value Ĥ /N c By considering the fact that Λ and D q contain the poles of G ′ and G, respectively, we conclude that, in the Thus, Eq. (32) is equivalent to Eq. (1) in the introduction in the T = 0 limit. An extension to T > 0 is straight forward.
The connected Green's function now contains anomalous contributions, but formally is obtained as in Eq. (19), where we have neglected the shifts Γ and X 0 since they only contribute to disconnected parts. Notice that Eq. (33) is a 2L × 2L matrix in Nambu and cluster-site space. The q vectors above refer to the reduced Brillouin zone originating from the cluster tiling, therefore G is expressed in a mixed representation. In translationinvariant systems, the Green's function is expected to be diagonal in the wave vectors k of the full Brillouin zone. This symmetry is notoriously broken in cluster methods such as VCA or C-DMFT. In order to obtain a k-diagonal 2 × 2 Nambu Green's function G(k, ω) we need to apply a periodization prescription. 29 This gives and r i is the position of site i within the cluster.
A nontrivial test for VCA is the noninteracting limit, for which this approximation becomes exact. In Appendix C we carry out this check for the noninteracting BH model, i. e., we set U = 0, and for a reference system consisting of single-site clusters. In this test case the grand potential Ω of the physical system can be evaluated analytically both using the VCA prescription as well as directly from the Hamiltonian of noninteracting lattice bosons.

C. Particle density and momentum distribution
The total particle density is defined as where N = N c L is the total number of lattice sites present in the physical system. The particle density can be easily expressed in Nambu formalism where the last term describes the contribution from the condensate, which can be deduced from Eq. (30). The term with the sum over q can be rewritten to obtain the known form of the particle density 11 The momentum distribution n(k) can be extracted by the Fourier transform within the cluster leading to where v † k is given by Eq. (34).

IV. APPLICATION TO THE BOSE-HUBBARD MODEL
In this section, we present the first nontrivial application of the extended VCA theory to the two-dimensional BH (BH) model and compare the results with unbiased quantum Monte Carlo (QMC) calculations. The BH Hamiltonian, 15,18 which describes strongly correlated lattice bosons, readŝ where a † i (a i ) creates (destroys) a bosonic particle and n i = a † i a i counts the number of particles at lattice site i. The parameter t is the hopping strength, which originates from the overlap of the localized wave functions belonging to lattice sites i and j, respectively. The first sum (indicated by angle brackets) is restricted to ordered pairs of nearest neighbor sites. The repulsive on-site interaction is termed U , and µ is the chemical potential, which controls the particle number. For increasing ratio t/U the system undergoes a quantum phase transition from the Mott to the superfluid phase. We evaluate static quantities, such as the particle density n and the condensate density n c as well as the dynamic single-particle spectral function A(k, ω). The phase boundary of the first three Mott lobes as obtained in VCA is shown in Fig. 1. The data for the first two lobes have been published in Ref. 11. Static quantities are evaluated for constant hopping strength t/U = 0.02 and distinct values of the chemical potential µ/U ranging from 0 to 3, scanning through various Mott lobes separated by the superfluid phase; see the dashed line in Fig. 1. The single-particle spectral function is evaluated for the parameter set marked by x in Fig. 1, which is located in the superfluid phase close to the tip of the first Mott lobe. For the numerical evaluation we used the chemical potential µ ′ and the strength of the source-and-drain coupling term F of the reference system as variational parameters. If not stated differently, the reference system consists of a cluster of size L = 2 × 2.   Fig. 1 (t/U = 0.02). The inset compares VCA and QMC results for the condensate fraction nc/n. VCA results are obtained for reference systems of size L = 1 × 1 and L = 2 × 2 and essentially infinitely large physical systems. QMC results are obtained for physical systems of size 32 × 32 inverse temperature U/T = 128.
The total particle density n evaluated using Eq. (35) is shown in Fig. 2 along with the condensate density n c = A † A /2L, and the density of the particles which are not condensed n − n c . From Fig. 3 it can be observed that the particle density n evaluated for reference systems of size L = 1 × 1 and of size L = 2 × 2 are almost identical. The same holds for the condensate fraction n c /n, which is shown in the inset of Fig. 3. In the same figure, we also compare our results with QMC calculations. The densities obtained from the two methods show an excellent agreement. The QMC data have been obtained for a system of size 32 × 32 and temperature U/T = 128 using the ALPS library 30 and the ALPS applications. 31 The single-particle spectral function A(k, ω) evaluated for the parameter set, marked by x in Fig. 1, i. e., in the superfluid phase close to the tip of the first Mott lobe, is depicted in Fig. 4. The colored density plot corresponds to VCA results and the dots with errorbars to latest QMC results of Ref. 32. The VCA spectral function A(k, ω) consists of four bands, which is in agreement with results obtained by means of a variational mean field calculation, 33 a strong coupling approach, 34 and random phase approximation (RPA) calculations. 35,36 The advantage of VCA in comparison to the above mentioned approaches is that the results can be systematically improved by increasing the cluster size of the reference system. For each wave vector k the weight is concentrated in one of the two bands present at positive and negative energy, respectively. We observe that the outer two modes exhibit a wide gap at k = 0, which is approximately of size U . The inner two, low-energy modes are also gapped at k = 0. However, the gap is tiny, and away from k = 0 the spectrum quickly develops a linear behavior, which is in agreement with the expected dispersion of Goldstone modes. The failure in obtaining a gapless long-wavelength excitation is a common problem of conserving approximations, i. e., of approximations for which macroscopic conservation laws are fulfilled. Similar aspects occur in dynamical mean-field theory calculations of two-component ultracold atoms as well. 37 In VCA there exists the additional possibility to systematically improve the obtained results by increasing the cluster size L of the reference system. Figure 5 compares the k = 0 gap of the inner modes for reference systems of size L = 1 × 1 and L = 2 × 2. The gap is evaluated along the dashed line shown in Fig. 1. The first observation is that the gap present in the condensed phase is almost an order of magnitude smaller than the gap in the Mott phase. It vanishes at the Mott-to-superfluid transition and, most importantly, shrinks with increasing cluster size L. This behavior signals convergence toward the correct result.
In Fig. 4 we also compare our VCA results for the  Fig. 1 for reference systems of size L = 1 × 1 and L = 2 × 2, respectively.
single-particle spectral function to latest QMC results obtained in Ref. 32. In this figure QMC results are indicated by dots with errorbars, which quantify the peak position of the spectral weight. Overall, we find good agreement in the low-energy spectrum. Only very close to k = 0 the two results differ slightly and the QMC dispersion possesses the correct gapless behavior. The QMC spectral function, exhibits only two instead of four bands. This is, however, not surprising since for the considered parameter set and at a specific wave vector k the weight of one positive (negative) energy band dominates drastically over the other one located at positive (negative) energy. Thus the four bands are extremely difficult to resolve by means of the maximum entropy method, which has been used to infer the spectra from QMC data; see Ref. 32 for details concerning the QMC results. This reference also contains a comparison between VCA data and QMC data for the spectral function evaluated in the Mott phase, where the results obtained from the two approaches coincide very well for all k values. We also evaluated the particle density n for the parameter set in the superfluid phase used in Fig. 4 and compared it to the QMC results. VCA yields n = 1.0321 in excellent agreement with the QMC result n QMC = 1.03068(2) obtained at U/T = 128 for a system of size 32 × 32.
In the following we provide some additional remarks on the fact that the single-particle spectral function, obtained within our approach, is gapped in the long wavelength limit, i. e., close to k = 0, for modes which ought to be identified with the Goldstone modes. This issue is a commonly known problem of conserving approximations. 38 One condition for an approximation to be conserving is, for example, to be Φ-derivable and self-consistent (see Ref. [39][40][41][42][43] for details). VCA is Φderivable but not self-consistent: The self-energy is the derivative of a functional of the Green's function, but the latter is not the Green's function obtained from Dyson's equation. Thus VCA is not completely conserving. However, many conservation laws are fulfilled at the stationary point of the self-energy functional, depending on which variational parameters are taken into account (see Ref. 44 for a more detailed discussion).
To obtain a gapless spectrum, a system of condensed bosons has to fulfill an independent condition, which is the Hugenholtz-Pines theorem. 42,45,46 There are only very few systematic approximation schemes which satisfy both conditions simultaneously. One notable exception occurs for interacting bosons composed of paired fermions. In this case, a consistent and gapless approximation can be developed provided the theory is expressed in terms of the constituent fermions. 47 In a different work 48 it was suggested to include an additional Lagrange multiplier in the form of a chemical potential, in order to explicitly enforce the Hugenholtz-Pines condition. Unfortunately, the Hugenholtz-Pines theorem is not fulfilled in VCA, and thus the low-energy modes of the single-particle spectral function are gapped in the long wavelength limit. Yet, the gap present in the VCA single-particle spectral function is small, and the spectrum quickly develops a linear behavior reminiscent of the gapless and linear Goldstone modes. Furthermore, in VCA there exists the possibility to systematically improve the results by increasing the cluster size of the reference system.
It is also interesting to mention that the related strong coupling approximation RPA, which yields a gapless sectrum, yet is not conserving 35,36 can be obtained within certain limits of the extended VCA formalism. Specifically, the limits to consider are (i) to use clusters of size L = 1 × 1, (ii) not to use the chemical potential µ as variational parameter and (iii) to determine the source-and-drain coupling strength F self-consistently within a mean-field approach, whereby intercluster hopping terms a † i a j are replaced with their mean-field value a † i a j +a † i a j in the reference Hamiltonian. This leads to the selfconsistency condition F = z t A , where A is given by Eq. (30) and z is the coordination number of the lattice. Our formalism provides a natural way to improve on RPA in a gapless, yet nonconserving, way by simply increasing the cluster size L and fixing F using the mean-field condition discussed above. However, it has to be emphasized that VCA yields much better results than RPA, even if RPA is extended to clusters of size L. Specifically, the particle density, the condensate density and the location of the phase boundary 11,49 can be determined much more accurately by means of VCA only because we allowed for a variation in the chemical potential µ ′ , i. e., allowed for macroscopic conservation laws to be fullfilled.

V. CONCLUSIONS
In the present paper, we introduced a pseudoparticle formalism for interacting bosonic systems, and showed that the results of the variational cluster approach can be derived within this formalism. We used it to extend the variational cluster approach to the superfluid phase of strongly correlated lattice bosons. We derived expressions for the grand potential and for other quantities, which are necessary to investigate the superfluid properties. Our results suggest that the pseudoparticle formalism is a quite versatile approach, which can be applied to a large variety of other problems.
As a first nontrivial application of the extended version of the variational cluster approach we choose the twodimensional Bose-Hubbard model and evaluated static quantities such as the total particle density and the condensate density, as well as the dynamic single-particle spectral function. We compared the single-particle spectral function with recent Quantum Monte-Carlo results 32 and found good agreement between the two approaches. It has to be pointed out that our extended variational cluster approach, while fulfilling many conservation laws, does not fulfill the Hugenholtz-Pines theorem. From this fact follows that the low-energy excitations of the spectrum have a small but nonzero gap in the long wavelength limit. This is a common aspect, which is already present in theories of the dilute Bose gas. 38,50,51 However, for wavevectors away from k = 0 the spectra obtained within this approach quite soon exhibit a correct linear behavior and agree very well with the Quantum Monte-Carlo results. Moreover, the gap shrinks with increasing cluster size, corroborating that the variational cluster approach becomes exact in the infinite cluster limit. Due to the fact this approach fulfills several conservation laws, the particle density, the condensate density as well as the phase boundary 11,49 delimiting the Mott from the superfluid phase can be evaluated very accurately. In the present paper we demonstrated, that our variational cluster approach results for the densities evaluated in both, the Mott and the superfluid phase, match perfectly with Quantum Monte-Carlo results.
The first condition of Eq. (16) requires that the Hermitian matrix Commuting Hermitian matrices have a common set of orthonormal eigenvectors. The matrix D is already diagonal. Hence for indices belonging to nondegenerate eigenvalues, X is also diagonal. Within the set of indices belonging to a degenerate eigenvalue, the corresponding Hermitian submatrix of X can be diagonalized by a unitary transformation U . In the following we term the diagonalized matrix as X ′ . The diagonalization also results in a new matrix V = V U of eigenvectors. We still have V −1 M V = D, but now . . , x ′ ns d ns ) . For the condition Eq. (16) we still need to ensure that x ′ α = ±1. Provided no x ′ α vanishes, which we will show below, this can easily be achieved by a suitable normalization of the column vector of V →Ṽ = V Z, with Z being a diagonal matrix, defined as Z αα ≡ 1/ |x ′ α |. We eventually havẽ . . . , e ns ) .
We are merely left with the proof that wherev α stands for the αth column of V . To this end we assume ad absurdum thatv α † Sv α = 0. In this case, v α would belong to the (n s − 1)-dimensional space S α orthogonal to the vector Sv α . According to Eq. (A1) the vectorsv 1 , . . .v α−1 ,v α+1 , . . .v ns also belong to S α and they are linear independent. Thus they span S α . Due to the fact that all vectorsv 1 , . . .v ns are linear independent,v α cannot belong to S α , which proves Eq. (A2).
Writing down the expressions in matrix form yields Using the expressions above and the relation A = QB +Γ we obtain for the pseudoparticle operators In order to evaluate part B we need the matrix M q , which is given by where we used thatT q = t q 0 0t T −q andt q =t T −q = µ ′ − µ − 2t cos q. Since M q is already diagonal we can readily evaluate part B as sum over the negative eigenvalues, which is µ + 2 t cos q, since µ < −2t. When summing over half of the q values the second term of the eigenvalue containing cos q is zero. For the calculation of part F we need the inverse of M 0 , which is simply andF , which reads Collecting all terms yields the grand potential Ω = 0, which is identical to the result obtained from the direct calculation.