Nonperturbative renormalization-group approach to strongly-correlated lattice bosons

# Nonperturbative renormalization-group approach to strongly-correlated lattice bosons

A. Rançon and N. Dupuis Laboratoire de Physique Théorique de la Matière Condensée, CNRS UMR 7600,
Université Pierre et Marie Curie, 4 Place Jussieu, 75252 Paris Cedex 05, France
October 10, 2011
###### Abstract

We present a nonperturbative renormalization-group approach to the Bose-Hubbard model. By taking as initial condition of the renormalization-group flow the (local) limit of decoupled sites, we take into account both local and long-distance fluctuations in a nontrivial way. This approach yields a phase diagram in very good quantitative agreement with quantum Monte Carlo simulations, and reproduces the two universality classes of the superfluid–Mott-insulator transition. The critical behavior near the multicritical points, where the transition takes place at constant density, agrees with the original predictions of Fisher et al. [Phys. Rev. B 40, 546 (1989)] based on simple scaling arguments. At a generic transition point, the critical behavior is mean-field like with logarithmic corrections in two dimensions. In the weakly-correlated superfluid phase (far away from the Mott insulating phase), the renormalization-group flow is controlled by the Bogoliubov fixed point down to a characteristic (Ginzburg) momentum scale which is much smaller than the inverse healing length . In the vicinity of the multicritical points, when the density is commensurate, we identify a sharp crossover from a weakly- to a strongly-correlated superfluid phase where the condensate density and the superfluid stiffness are strongly suppressed and both and are of the order of the inverse lattice spacing.

###### pacs:
05.30.Jp, 05.10.Cc, 05.30.Rt

## I Introduction

In the last two decades, the nonperturbative renormalization group (NPRG) approach has been successfully applied to many areas of physics,Berges et al. (2002); Delamotte (2007) from high-energy physics to statistical and condensed-matter physics. It has proven to be a powerful tool to study not only the low-energy long-distance properties in the vicinity of second-order phase transitions but also non-universal quantities. In particular, the NPRG approach has been implemented in lattice models and used to compute the transition temperature and the magnetization in classical spin models (Ising, XY and Heisenberg models).Machado and Dupuis (2010) This implementation of the NPRG is referred to as the lattice NPRG.

The strategy of the NPRG is to build a family of models indexed by a momentum scale , such that fluctuations are smoothly taken into account as is lowered from a microscopic scale down to 0. In practice this is achieved by adding to the action of the system an infrared regulator term which vanishes for . For a scalar field theory, the regulator term is a mass-like term , where the cutoff function is chosen such that for and for , which effectively suppresses the low-energy modes . One can then define a scale-dependent partition function and a scale-dependent effective action defined as a slightly modified Legendre transform (see Sec. II.1 for the precise definition) of . Here is an external source which couples linearly to the field and . In the standard implementation of the NPRG, at the microscopic scale , all fluctuations are frozen by the term so that as in Landau’s (mean-field) theory of phase transitions. The effective action of the original model is obtained for () and can be determined by (approximately) solving the RG equation satisfied by .Berges et al. (2002); Delamotte (2007)

The lattice NPRG differs from the standard implementation in the initial condition.Machado and Dupuis (2010) The cutoff function is chosen such that at the microscopic scale the action corresponds to the local limit of decoupled sites. Local fluctuations are therefore included from the very beginning of the RG procedure. The intersite coupling is then gradually restored as decreases from down to 0. In the low-energy limit , acts as an infrared regulator suppressing fluctuations with momenta . The lattice NPRG is then equivalent to the standard NPRG and yields identical results for the critical properties. The hallmark of the lattice NPRG is thus to take into account both local and critical fluctuations in a nontrivial way.

In this paper, we present a NPRG study of the Bose-Hubbard modelFisher et al. (1989) at zero temperature and in dimension or . This model has been intensively studied in the last years following the experimental observation of the superfluid–Mott-insulator transition of an ultracold bosonic gas in an optical lattice.Jaksch et al. (1998); Greiner et al. (2002); Stöferle et al. (2004); Spielman et al. (2007) Phase diagram and thermodynamic quantities are known from the numerically exact lattice quantum Monte Carlo (QMC) simulations.Capogrosso-Sansone et al. (2007, 2008) On the other hand few studies have addressed the critical behavior at the superfluid–Mott-insulator transition,Krauth and Trivedi (1991) and most of our understanding goes back to the seminal work of Fisher et al.Fisher et al. (1989)

The standard NPRG scheme does not capture the superfluid–Mott-insulator transition in the Bose-Hubbard model. The reason is that near the transition the mean-field solution is too far away from the actual state of the system to provide a reliable initial condition for the NPRG procedure. The two-pole structure of the local (on-site) single-particle propagator is crucial for the very existence of the transition (see, e.g., Ref. Fisher et al., 1989). It is however impossible to reproduce this structure from a RG approach starting from the mean-field (Bogoliubov) theory within standard approximations of the RG equation satisfied by the effective action . This prevents a straightforward generalization of recent NPRG studiesDupuis and Sengupta (2007); Dupuis (2009a, b); Wetterich (2008); Floerchinger and Wetterich (2008); Sinner et al. (2009, 2010) of interacting bosons to the Bose-Hubbard model.

By contrast the lattice NPRG, which takes into account local fluctuations, is able to describe the superfluid–Mott-insulator transition.Rançon and Dupuis (2011a) Since the starting action is purely local, this approach is to some extent reminiscent of various expansions of the Bose-Hubbard model.Freericks and Monien (1994); *Freericks96; *Freericks09; Buonsante and Vezzani (2005); dos Santos and Pelster (2009); Teichmann et al. (2009a); *Teichmann09b; Koller and Dupuis (2006); Knap et al. (2010, 2011); Arrigoni et al. (2011) Moreover, the lattice NPRG is not restricted to the computation of thermodynamic quantities and allows us to study the critical behavior at the superfluid–Mott-insulator transition and compare with the predictions of Fisher et al.Fisher et al. (1989) based on scaling arguments.

In addition to the phase diagram and the critical behavior at the superfluid–Mott-insulator transition, the NPRG approach can also address the superfluid phase. Deep in the superfluid phase, localization effects are negligible and we expect the Bogoliubov theory to provide a good description of the system. However, even in this weak correlation limit, it is known that the Bogoliubov approximation breaks down below a characteristic (Ginzburg) momentum scale . In perturbation theory about the Bogoliubov approximation, the Ginzburg scale manifests itself by the appearance of infrared divergences below three dimensions (). Although these divergences cancel out in local gauge invariant quantities (condensate density, sound mode velocity, etc.),Beliaev (1958a, b); Hugenholtz and Pines (1959); Gavoret and Nozières (1964) they do have a physical origin: they result from the coupling between longitudinal and transverse (phase) fluctuations and reflect the divergence of the longitudinal susceptibilityNepomnyashchii and Nepomnyashchii (1978); Nepomnyashchii (1983) – a general phenomenon in systems with a continuous broken symmetry.Patasinskij and Pokrovskij (1973); Fisher et al. (1973); Anishetty et al. (1999); Sachdev (1999a); Zwerger (2004); Dupuis (2011) The normal and anomalous self-energies, and , are non-analytic functions of and when and ( denotes the velocity of the sound mode), while vanishes,Nepomnyashchii and Nepomnyashchii (1975) in marked contrast with the Bogoliubov approximation where the linear spectrum and the superfluidity rely on a finite value of the anomalous self-energy. A weakly-correlated superfluid is defined by the condition where the healing scale is the inverse of the healing length .Pistolesi et al. (2004); Dupuis (2009b, 2011); Capogrosso-Sansone et al. (2010) In this case, the Bogoliubov theory applies to a large part of the spectrum where the dispersion is linear () and breaks down only at very low momenta . The Goldstone regime , dominated by phase fluctuations, is conveniently described by Popov’s hydrodynamic theory (free of infrared divergences) based on a density-phase representation of the boson field .Popov (1972); Popov and Seredniakov (1979); Popov (1987) The NPRG approach yields a unified description of superfluidity which includes both Bogoliubov theory (valid for ) and Popov’s hydrodynamic approach (valid for ).Dupuis (2009a, b, 2011) The Bose-Hubbard model gives us the opportunity to understand the fate of the weakly-correlated superfluid phase as we increase the strength of the interactions and move closer to the Mott insulating phase in the phase diagram.

The paper is organized as follows. In Sec. II, we derive the lattice NPRG formalism for the Bose-Hubbard model. We introduce the scale-dependent effective action and compute its initial value by solving the single-site Bose-Hubbard model. We show that reproduces the result of the strong-coupling random-phase approximation (RPA).Sheshadri et al. (1993); van Oosten et al. (2001); Sengupta and Dupuis (2005); Ohashi et al. (2006); Menotti and Trivedi (2008) We also discuss the approximations used to solve the flow equation satisfied by . The phase diagram obtained from the NPRG equations is in very good quantitative agreement with the QMC results (Sec. III). Furthermore, the critical behavior derived from the NPRG analysis is in complete agreement with the predictions of Fisher et al. based on scaling arguments (Sec. IV).Fisher et al. (1989) We find multicritical points in the universality class of the -dimensional XY model where the transition takes place at constant density. The XY critical behavior is observed in the Mott gap, the condensate density, the compressibility and the superfluid stiffness when a multicritical point is approached at constant chemical potential by varying the ratio between the hopping amplitude and the local repulsion between particles. At a generic transition point, we observe mean-field behavior, with logarithmic corrections in dimension (corresponding to the upper critical dimension). The superfluid phase is discussed in Sec. V. In the dilute limit, the renormalization-group flow is controlled by the Bogoliubov fixed point down to a characteristic (Ginzburg) momentum scale which is much smaller than the inverse healing length . The Goldstone regime , dominated by phase fluctuations, is characterized by a (relativistic) Lorentz invariance of the effective action .Wetterich (2008); Dupuis and Sengupta (2007) In the vicinity of the multicritical points, when the density is commensurate, we identify a sharp crossover from a weakly- to a strongly-correlated superfluid phase where the condensate density and the superfluid stiffness are strongly suppressed and both and are of the order of the inverse lattice spacing. The main results are summarized in Sec. VI.

## Ii Lattice NPRG

The Bose-Hubbard model on a -dimensional hypercubic lattice is defined by the (Euclidean) action

 S=∫β0dτ{ ∑r[ψ∗r(∂τ−μ)ψr+U2(ψ∗rψr)2] −t∑⟨r,r′⟩(ψ∗rψr′+c.c.)}, (1)

where is a complex field and an imaginary time with the inverse temperature. denotes the sites of the lattice. is the on-site repulsion, the hopping amplitude between nearest-neighbor sites and the chemical potential. In the following we will sometimes write the boson field

 ψr=1√2(ψ1r+iψ2r) (2)

in terms of two real fields and .

We set and take the lattice spacing as the unit length throughout the paper.

### ii.1 Scale-dependent effective action

Following the general strategy of the NPRG, we consider a family of models with action indexed by a momentum scale varying from a microscopic scale down to 0. The regulator term is defined by

 ΔSk=∫β0dτ∑qψ∗qRk(q)ψq, (3)

where is the Fourier transform of and the sum over runs over the first Brillouin zone of the reciprocal lattice. The cutoff function modifies the bare dispersion of the bosons. is chosen such that the effective (bare) dispersion vanishes.Machado and Dupuis (2010) The action then corresponds to the local limit of decoupled sites (vanishing hopping amplitude). By choosing , rather than as in Ref. Machado and Dupuis, 2010, we ensure that does not modify the chemical potential but only the kinetic energy.

In practice, we choose the cutoff function

 Rk(q)=−ZA,kϵksgn(tq)(1−yq)Θ(1−yq), (4)

with , , and the step function (see Fig. 1). The -dependent constant is defined below (). Since , the action coincides with the action (1) of the original model. For small , the function gives a mass to the low-energy modes and acts as an infrared regulator as in the standard NPRG scheme.Berges et al. (2002); Delamotte (2007)

The scale-dependent effective action

 Γk[ϕ∗,ϕ]= −lnZk[J∗,J]+∫β0dτ∑r(J∗rϕr+c.c.) −ΔSk[ϕ∗,ϕ] (5)

is defined as a (slightly modified) Legendre transform which includes the explicit subtraction of . Here is the partition function obtained from the action , a complex external source which couples linearly to the bosonic field , and

 ϕr(τ)=δlnZk[J∗,J]δJ∗r(τ),ϕ∗r(τ)=δlnZk[J∗,J]δJr(τ) (6)

the superfluid order parameter. The variation of the effective action with is governed by Wetterich’s equation,Wetterich (1993)

 ∂kΓk[ϕ∗,ϕ]=12Tr{∂kRk(Γ(2)k[ϕ∗,ϕ]+Rk)−1}, (7)

where is the second-order functional derivative of . In Fourier space, the trace in (7) involves a sum over momenta and frequencies as well as the two components of the complex field .

We are primarily interested in two quantities. The first one is the effective potential defined by

 Vk(n)=1βNΓk[ϕ∗,ϕ]∣∣∣ϕconst (8)

where is a constant (uniform and time-independent) field. The U(1) symmetry of the action implies that is a function of . Its minimum determines the condensate density and the thermodynamic potential (per site) in the equilibrium state.

The second quantity of interest is the two-point vertex

 Γ(2)k,ij(r−r′,τ−τ′;ϕ)=δ(2)Γ[ϕ]δϕir(τ)δϕjr′(τ′)∣∣∣ϕconst (9)

which determines the one-particle propagator . Here the indices refer to the real and imaginary parts of [see Eq. (2)]. Because of the U(1) symmetry of the action (1), the two-point vertex in a constant field takes the formDupuis (2009b)

 Γ(2)k,ij(q;ϕ)=δi,jΓA,k(q;n)+ϕiϕjΓB,k(q;n)+ϵijΓC,k(q;n) (10)

in Fourier space, where , is a Matsubara frequency and the antisymmetric tensor. For , we can relate to the derivative of the effective potential,

 Γ(2)k,ij(q=0;ϕ)=∂2Vk(n)∂ϕi∂ϕj=δi,jV′k(n)+ϕiϕjV′′k(n), (11)

so that

 ΓA,k(q=0;n)=V′k(n),ΓB,k(q=0;n)=V′′k(n),ΓC,k(q=0;n)=0. (12)

Parity and time-reversal invariance implyDupuis (2009b)

 ΓA,k(q;n)=ΓA,k(−q;n)=ΓA,k(q,−iω;n),ΓB,k(q;n)=ΓB,k(−q;n)=ΓB,k(q,−iω;n),ΓC,k(q;n)=−ΓC,k(−q;n)=−ΓC,k(q,−iω;n). (13)

The one-particle propagator can be written in a form analogous to (10) or in terms of its longitudinal and transverse components,

 Gk,ij(q;ϕ)= ϕiϕj2nGk,ll(q;n)+(δi,j−ϕiϕj2n)Gk,tt(q;n) +ϵijGk,lt(q;n), (14)

where

 Gk,ll(q;n)=−ΓA,k(q;n)Dk(q;n),Gk,tt(q;n)=−ΓA,k(q;n)+2nΓB,k(q;n)Dk(q;n),Gk,lt(q;n)=ΓC,k(q;n)Dk(q;n), (15)

with . Note that the single-particle propagator entering the flow equation (7) is defined by , which is the propagator associated with the true Legendre transform, rather than .

### ii.2 Initial conditions

Since the action corresponds to the local limit, the initial value of the effective action reads

 ΓΛ[ϕ∗,ϕ]=Γloc[ϕ∗,ϕ]+∫β0dτ∑qϕ∗(q)tqϕ(q), (16)

where

 Γloc[ϕ∗,ϕ]=−lnZloc[J∗,J]+∫β0dτ∑r(J∗rϕr+c.c.) (17)

is the Legendre transform of the thermodynamic potential in the local limit. In Eq. (17), is related to by the relation and is the partition function obtained from .

It is not possible to compute the functional for arbitrary time-dependent fields. One can however easily obtain the effective potential and the two-point vertex in a time-independent field . These quantities are sufficient to specify the initial conditions of the flow within the approximations that we consider in Sec. II.5.

To obtain and in a time-independent field, it is sufficient to consider a single site with time-independent complex external source . The corresponding Hamiltonian reads

 ^H=−μ^n+U2^n(^n−1)−J∗^b−J^b†, (18)

where () is a creation (annihilation) operator and . In the basis [ with integer], the Hamiltonian is represented by a tridiagonal matrix,

 ⟨m|^H|m′⟩= δm,m′[−μm+U2m(m−1)] −δm+1,m′J∗√m+1−δm−1,m′J√m, (19)

which can be numerically diagonalized in the truncated Hilbert space . The low-energy eigenstates are independent of if the latter is large enough. If we denote by the source-dependent eigenstates and eigenvalues – with the ground state – we obtain the superfluid order parameter

 ϕ=−∂E0∂J∗,ϕ∗=−∂E0∂J, (20)

and the effective potential

 Vloc(n)=E0+J∗ϕ+Jϕ (21)

() in the zero-temperature limit . Figure 2 shows the superfluid order parameter as a function of the external source , and the local effective potential obtained by numerically inverting (20). The special case where the ground state in the local limit is degenerate for ( integer) is discussed in Appendix A.

To determine the two-point vertex , we start from the (source-dependent) normal and anomalous local Green functions

 Gn(τ)=−⟨Tτ^b(τ)^b†(0)⟩+|⟨^b⟩|2,Gan(τ)=−⟨Tτ^b(τ)^b(0)⟩+⟨^b⟩2, (22)

where and is a time-ordering operator. The Fourier transforms and are easily expressed in terms of the eigenstates of the Hamiltonian,

 Gn(iω)=−∑α≠0[|⟨α|^b|0⟩|2iω+Eα−E0−|⟨0|^b|α⟩|2iω+E0−Eα],Gan(iω)=−∑α≠0⟨α|^b|0⟩⟨0|^b|α⟩2(Eα−E0)ω2+(Eα−E0)2. (23)

From the relation , we obtain

 Γloc,A(iω;n)=−12D[Gn(iω)+Gn(−iω)+2Gan(iω)],Γloc,B(iω;n)=Gan(iω)nD,Γloc,C(iω;n)=i2D[Gn(iω)−Gn(−iω)], (24)

where . is expressed in terms of the condensate density (rather than the external source ) by inverting (20).

The large frequency limit of the two-point vertex is given bynot (a)

 lim|ω|→∞Γloc,A(iω;n)=−μ−U⟨ψ2r⟩+2U⟨ψ∗rψr⟩,lim|ω|→∞Γloc,B(iω;n)=U⟨ψ2r⟩n,lim|ω|→∞Γloc,C(iω;n)=ω, (25)

where and are assumed real (which corresponds to a real external source ),

 ⟨ψ∗rψr⟩ =−Gn(τ=0−)+|ϕ|2 =−∫∞−∞dω2πGn(iω)eiω0++|ϕ|2, (26)

and

 ⟨ψ2r⟩ =−Gan(τ=0−)+ϕ2 =−∫∞−∞dω2πGan(iω)eiω0++ϕ2. (27)

The asymptotic forms (25) are reached for . and are shown in Figs. 4 and 4.

#### ii.2.1 Large-field limit

To obtain the large-field limit of the effective potential , we must compute the local partition function for . Using a loop expansion about the saddle-point approximation (see Appendix B), we find

 Vloc(n)=−¯μn+U2n2+O(n0), (28)

where .

#### ii.2.2 Strong-coupling RPA

The initial effective action [Eq. (16)] treats the local fluctuations exactly but includes the intersite hopping term at the mean-field level, thus reproducing the strong-coupling RPA.Sheshadri et al. (1993); van Oosten et al. (2001); Sengupta and Dupuis (2005); Ohashi et al. (2006); Menotti and Trivedi (2008) The effective potential reads

 VΛ(n)=Vloc(n)−2dtn, (29)

while the two-point vertex takes the RPA-like form

 Γ(2)Λ,ij(q;n)=Γ(2)loc,ij(iω;n)+δi,jtq. (30)

 VΛ(n)=Vloc(0)+[Γ(2)loc,11(iω=0;n=0)−2dt]n+O(n2), (31)

where

 Γ(2)loc,ii(iω=0;n=0)=−Gn(iω=0;n=0)−1 (32)

is determined by the local Green function

 Gn(iω;n=0)=¯nloc+1iω+μ−U¯nloc−¯nlociω+μ−U(¯nloc−1) (33)

for vanishing source (). Here is the number of bosons per site in the local limit: if and if . The ground state is a Mott insulator as long as . Thus the transition to the superfluid state is determined by the criterion , i.e.

 Gn(iω=0;n=0)−1+2dt=0, (34)

which reproduces the mean-field (or strong-coupling RPA) phase diagram.Fisher et al. (1989) Equation (34) can also be obtained from the condition , which signals the appearance of a pole at zero momentum and frequency in the one-particle propagator .

In the strong-coupling RPA, the condensate density in the superfluid phase is determined by

 V′Λ(n0)=V′loc(n0)−2dt=0. (35)

The hopping amplitude acts as a source term for the local potential . For , i.e. deep in the superfluid phase, the source term is large and we are effectively in the large field limit discussed in Sec. II.2.1. From Eqs. (16), and the fact that when , we then obtain (with ), which is nothing but the result of the Bogoliubov approximation.Dupuis (2009b); not (b) The strong-coupling RPA reduces to the Bogoliubov theory in the limit .Menotti and Trivedi (2008) Table 1 compares the initial conditions given by the Bogoliubov theory and the strong-coupling RPA.

It should be noted that the true Legendre transform is for since the action is local. The lattice NPRG is an expansion about the local limit. The scale-dependent effective action is however the right quantity to consider to analyze the physical properties of the system with action (without the regulator term). In , the regulator term is compensated, in a mean-field manner, by subtracting from the true Legendre transform. It follows that the physical quantities at scale , such as the condensate density , are obtained from rather than from the true Legendre transform.

### ii.3 Gauge invariance and Ward identities

The invariance of the action in the local (time-dependent) gauge transformation , and imposes important constraints on the effective action . In the superfluid phase, this implies that the two-point vertex satisfies the Ward identitiesDupuis (2009b)

 ∂∂ωΓC,k(q;n0,k)∣∣∣q=0=−∂2Vk∂n∂μ∣∣∣n0,k,∂2∂ω2ΓA,k(q;n0,k)∣∣∣q=0=−12n0,k∂2Vk∂μ2∣∣∣n0,k, (36)

where the effective potential is considered as a function of both and , the condensate density being then defined by

 ∂Vk(n,μ)∂n∣∣∣n0,k=0. (37)

Since Eq. (37) is valid for any , we deduce the relation

 0=ddμ∂Vk∂n∣∣∣n0,k=∂2Vk∂n∂μ∣∣∣n0,k+∂2Vk∂n2∣∣∣n0,kdn0,kdμ, (38)

which will be used below together with (36).

In the Mott insulator (), the Ward identities (36) become

 ∂∂ωΓC,k(q;n=0)∣∣∣q=0=−∂2Vk∂n∂μ∣∣∣n=0,∂2V0,k∂μ2=d2V0,kdμ2=0, (39)

which implies that the compressibilitynot (c)

 κk=d¯nkdμ=−d2V0,kdμ2 (40)

vanishes. denotes the boson density (mean boson number per site).

### ii.4 Derivative expansion and infrared behavior

The low-energy behavior of the system is best understood from a derivative expansion of the two-point vertex. Since the cutoff function (4) acts as an infrared regulator, is a regular function of for . In the infrared limit, we can therefore use the derivative expansion

 ΓA,k(q;n)=ZA,k(n)tq2+VA,k(n)ω2+V′k(n),ΓB,k(q;n)=V′′k(n),ΓC,k(q;n)=ZC,k(n)ω, (41)

(, etc.) in agreement with the symmetry properties (13). For the following discussion, it is convenient to introduce

 δk=∂Vk∂n∣∣∣n0,k,λk=∂2Vk∂n2∣∣∣n0,k, (42)

with vanishing in the superfluid phase. If the spectrum is gapped, Eqs. (41) will always be valid for energy scales below the gap. Otherwise their validity requires and where is the lowest excitation energy for (see Sec. II.5).

#### ii.4.1 Superfluid phase

Using (37) and (38), we can rewrite the Ward identities (36) as

 ZC,k(n0,k)=λkdn0,kdμ,VA,k(n0,k)=−12n0,k∂2Vk∂μ2∣∣∣n0,k, (43)

while the compressibility (40) is expressed as

 κk =−∂2V0,k∂μ2∣∣∣n0,k−∂2V0,k∂n∂μ∣∣∣n0,kdn0,kdμ =2n0,kVA,k+Z2C,kλk. (44)

The superfluid stiffness , defined as the rigidity with respect to a twist of the phase of the order parameter, can be obtained from the transverse part of the two-point vertex,Dupuis (2009b)

 ΓA,k(q,ω=0;n0,k)=ρs,k2n0,kq2(q→0), (45)

 ρs,k=2tZA,k(n0,k)n0,k. (46)

The excitation spectrum is given by the zeros of the determinant of the matrix (after analytical continuation ),

 detΓ(2)k(q) =ΓA,k(q)[ΓA,k(q)+2n0,kΓB,k(q)]+ΓC,k(q)2 ≃2λkn0,k(ZA,ktq2+VA,kω2)+(ZC,kω)2 (47)

(all quantities are evaluated for ) for . This equation yields a gapless (Goldstone) mode with a velocity

 ck =(ZA,k(n0,k)tVA,k(n0,k)+ZC,k(n0,k)2/(2λkn0,k))1/2 =(ρs,kκk)1/2 (48)

which can be expressed in terms of the compressibility and superfluid stiffness.Fisher et al. (1989) The existence of a gapless mode is a consequence of the Hugenholtz-Pines theoremHugenholtz and Pines (1959) which, in our formalism, readsDupuis (2009b)

 ΓA,k(q=0;n0,k)=V′k(n0,k)=0. (49)

Equations (46,48) are identical to those obtained in continuum models if we identify with the (effective) mass of the bosons in the lattice potential. From , we also obtain a gapped mode, with a gap which is however larger than , and therefore beyond the domain of validity of the derivative expansion (). The existence of two modes in the superfluid phase follows from being of order . Pushing the derivative expansion to higher order in would yield additional modes. These modes are not in the domain of validity of the derivative expansion and do not show up in the spectral function.Dupuis (2009b); Sinner et al. (2010)

#### ii.4.2 Mott insulator

In the Mott insulator (), the Ward identities (39) yield

 ZC,k=−dδkdμ. (50)

Since , the excitation spectrum is obtained from after analytical continuation . This gives two gapped modes

 ω±(q) =−ZC,k2VA,k±12VA,k[Z2C,k+4VA,k(ZA,ktq2+δk)]1/2 =Δk±±ZA,ktq2(Z2C,k+4VA,kδk)1/2+O(|q|4), (51)

where

 Δk±=−ZC,k2VA,k±12VA,k(Z