Auxiliary field approach to dilute Bose gases with tunable interactions

Auxiliary field approach to dilute Bose gases with tunable interactions

Abstract

We rewrite the Lagrangian for a dilute Bose gas in terms of auxiliary fields related to the normal and anomalous condensate densities. We derive the loop expansion of the effective action in the composite-field propagators. The lowest-order auxiliary field (LOAF) theory is a conserving mean-field approximation consistent with the Goldstone theorem without some of the difficulties plaguing approximations such as the Hartree and Popov approximations. LOAF predicts a second-order phase transition. We give a set of Feynman rules for improving results to any order in the loop expansion in terms of composite-field propagators. We compare results of the LOAF approximation with those derived using the Popov approximation. LOAF allows us to explore the critical regime for all values of the coupling constant and we determine various parameters in the unitarity limit.

pacs:
03.75.Hh, 05.30.Jp, 67.85.Bc
1

I Introduction

In 1911, Kamerlingh Onnes found that liquid He, when cooled below K began to expand rather than contractOnnes (1911). ÊThe transition, later named the lambda-transition was recognized in 1938 as the onset of superfluidityKapitza (1938); Allen and Misener (1938). ÊThe connection with Bose-Einstein condensation (BEC), first argued by F. London on the basis of the near identical values of the lambda transition temperature and the critical temperature for BEC of noninteracting bosonsLondon (1938a, b) sparked a series of weakly interacting BEC studies when BogoliubovBogoliubov (1947) pointed out that the BEC elementary excitations satisfy the Landau criterion for superfluidityLandau (1941). ÊIn the weakly interacting limit, the interactions can be characterized by a single parameterLee et al. (1957) — the scattering length  — giving the results a powerful, general applicability. The hope of studying bosons with short-range inter-particle interactions of a strength that can be tuned all the way from weakly interacting () to universality (), appeared thwarted when it was found that the three-body loss-rate in cold atom traps scales as near a Feshbach resonanceFedichev et al. (1996); Esry et al. (1999). In cold atom traps, only fermions have been obtained in the strongly interacting, quantum degenerate regime in equilibriumShin et al. (2007), in which case three-body loss is reduced by virtue of the Pauli exclusion principle. Recently, however, it was pointed outDaley et al. (2009) that three-body losses can be strongly suppressed in an optical lattice when the average number of bosons per lattice site is two or less. The development of novel cold atom technologyHenderson et al. (2006, 2009) leads to the prospect of studying finite temperature properties, such as the BEC transition temperature, , superfluid to normal fluid ratio, depletion, and specific heat, at fixed particle density .

At finite temperature the description of BEC’s even in the weakly interacting regime remains a challenge. ÊStandard approximations such as the Hartree-Fock-Bogoliubov (HFB) and the Popov schemes, generally fall within the Hohenberg and Martin classificationHohenberg and Martin (1965) of conserving and gapless approximations which imply that they either violate Goldstone’s (or Hugenholz-Pines) theorem or general conservation lawsHohenberg and Martin (1965). Both these approximations predict the BEC-transition to be a first-order transition, whereas we expect the transition to be second-orderAndersen (2004). The calculation of , first undertaken by ToyodaToyoda (1982) to explain the difference between the lambda-transition temperature ( K) and the ( K) of the noninteracting BEC at the same density, exemplifies the difficulties of understanding the theory near : whereas Toyoda found a -decrease with increasing scattering length, K. Huang later pointed out that the calculation had a sign error, giving an increasing value of Huang (1999). ÊHowever, Baym and collaboratorsBaym et al. (1999, 2000) noted that the Toyoda expansion involves an expansion in a large parameter. ÊTheir calculation found a linear increase of with . ÊThe fact that the helium lambda transition temperature falls below may be explained by quantum Monte-Carlo calculationsBaym et al. (1999), which found that the critical temperature of a hard-sphere boson gas increases at low values of , then turns over and drops below near .

In this paper, we discuss in detail a theoretical description that we introduced recentlyCooper et al. (2010) to describe a large interval of values, satisfies Goldstone’s theorem, yields a Bose-Einstein transition that is second-order, gives the same critical temperature variation found in Refs. Baym et al., 1999, 2000 but at a lower order of the calculation, while also predicting reasonable values for the depletion. ÊThis method then resolves many of the main challenges in describing boson physics over a large temperature and regime and it’s predictions will be available for experimental testing in the near future. The approach we present here is different from other resummation schemes such as the large- expansion (which is a special case of this expansion), in that it treats the normal and anomalous densities on an equal footing.

In the following, we will discuss the general features that arise when rewriting the original theory in terms of composite fields. One aspect of this approach is that one can systematically calculate corrections to the mean-field results presented earlierCooper et al. (2010) in a loop expansion in the composite-field propagators. We derive the Feynman rules for such an expansion using the propagators and vertices of the mean-field approximation. At each level of this loop expansion one maintains the features that the results are both gapless and conserving. The broken symmetry Ward identities guarantee Goldstone’s theorem order-by-order in the loop expansion in terms of auxiliary-field propagatorsBender et al. (1977).

In our auxiliary field formalism, we introduce two auxiliary fields related to the normal and anomalous densities by means of the Hubbard-Stratonovitch transformationHubbard (1959); Stratonovich (1958), utilizing methods discussed in the quantum field theory community Bender et al. (1977); Coleman et al. (1974); Root (1974). This transformation has already been shown to be quite useful in discussing the properties of the BCS-BEC crossover in the analogous 4-fermi theory for the BCS phase Sá de Melo et al. (1993); Engelbrecht et al. (1997); Floerchinger et al. (2008). The path integral formulation of the grand canonical partition function can be found in Negele and Orland Negele and Orland (1988). The Hubbard Stratonovich transformation is used to replace the original quartic interaction with an interaction quadratic in the original fields. An excellent review of previous use of path integral methods to study dilute Bose gases is found in the review article of AndersenAndersen (2004). The use of path integral methods to study various topics in dilute gases began with the work of Braaten and Nieto Braaten and Nieto (1997). Path integral methods have recently been used to study static and dynamical properties of the dilute Bose gases Rey et al. (2004); Gasenzer et al. (2005); Temme and Gasenzer (2006); Berges and Gasenzer (2007); Friederich et al. (2010); Floerchinger and Wetterich (2008); Floerchinger et al. (2008). An excellent summary of this approach and its connection to the more traditional Hamiltonian approach is to be found in the recent book by Calzetta and Hu Calzetta and Hu (2008). We also point out that the 1/N expansion, which is a special case of the method being proposed here, has a long history of use in high-energy and condensed matter physics Brezin and Wada (1993); Moshe and Zinn-Justin (2003). It has been used to calculate the critical temperature by Baym, Blaizot and Zinn-Justin Baym et al. (2000). This calculation gives the same result for as the method we are describing here. However, our approach can be used at all temperatures. Corrections to the 1/N result to calculating were obtained by Arnold and Tomasik Arnold and Tomasik (2000).

The paper is organized as follows: In Sec. II we discuss the auxiliary-field formalism and rewrite the Lagrangian for weakly interacting Bosons in terms of two auxiliary fields. In Sec. III we derive the loop expansion by performing the path integral over the original fields and then performing the resulting path integral over the auxiliary fields by stationary phase. In Sec. IV we find the leading-order loop expansion in the auxiliary fields (LOAF) for the action. In Sec. V we set the auxiliary-field parameter and discuss the leading-order effective potential for both the ground state and at finite temperature. In Sec. VI we discuss related mean-field approximations. In Sec. VII we discuss numerical results for the theory at finite temperature and varying dimensionless coupling constant . We compare the LOAF approximation to the Popov approximation in detail. We conclude in Sec. VIII. Finally, in App. A we discuss the connection between regularization of the effective potential and renormalization of the parameters. In App. B we give the rules for determining all the Feynman graphs for the expansion using the mean-field propagators and vertices.

Ii The auxiliary-field formalism

The classical action is given by

 S[ϕ,ϕ∗]=∫[dx]L[ϕ,ϕ∗], (1)

where and where the Lagrangian density is

 L[ϕ,ϕ∗] =iℏ2[ϕ∗(x)(∂tϕ(x))−(∂tϕ∗(x))ϕ(x)] (2) −ϕ∗(x){−ℏ2∇22m−μ0}ϕ(x)−λ02|ϕ(x)|4.

Here and are the bare (unrenormalized) chemical potential and contact interaction strength respectively. We introduce two auxiliary fields, a real field, , and a complex field, , by means of the Hubbard-Stratonovitch transformationHubbard (1959); Stratonovich (1958), utilizing methods discussed in Refs. Bender et al., 1977; Coleman et al., 1974; Root, 1974. In our case, the auxiliary-field Lagrangian density takes the form

 Laux[ϕ,ϕ∗,χ,A,A∗]=12λ0[χ(x)−λ0coshθ|ϕ(x)|2]2 −12λ0∣∣A(x)−λ0sinhθϕ2(x)∣∣2, (3)

which we add to Eq. (2). Here is a parameter which provides a mixing between the normal and anomalous densities. In Sec. VI, we will see that choosing leads to the usual large- expansion which has only the auxiliary field Coleman et al. (1974); Root (1974). In lowest order, gives a gapless solution very similar to the free Bose gas in the condensed phase. If instead we choose such that , then in the weak coupling limit our results agree with the Bogoliubov theoryBogoliubov (1947); Andersen (2004), which represents the leading-order low-density expansion. Of course all values of lead to a complete resummation of the original theory in terms of different combinations of the composite fields.

For an arbitrary parameter , the action is given by

 S[Φ,J]= (4) −12∬[dx][dx′]ϕa(x)G−1ab[χ](x,x′)ϕb(x′) +∫dx{χi(x)χi(x)2λ0+Φα(x)Jα(x)}.

with

 (5) G−10ab=(h−μ000h∗−μ0),h=−ℏ2∇22m−iℏ∂∂t, Vab[χ](x)=(χ(x)coshθ−A(x)sinhθ−A∗(x)sinhθχ(x)coshθ).

Here we have introduced a two-component notation using Roman indices for the fields and and currents and ,

 ϕa(x) =(ϕ(x),ϕ∗(x)), ϕa(x) =(ϕ∗(x),ϕ(x)), (6a) ja(x) =(j(x),j∗(x)), ja(x) =(j∗(x),j(x)), (6b)

for , and a three-component notation using Roman indices for the fields , , and ,

 χi(x) =(χ(x),A(x)/√2,A∗(x)/√2), (7) Si(x) =(s(x),S(x)/√2,S∗(x)/√2),

and

 χi(x) =(χ(x),−A∗(x)/√2,−A(x)/√2), (8) Si(x) =(s(x),−S∗(x)/√2,−S(x)/√2),

for . For convenience, we also define five-component fields with Greek indices and currents . These definitions define a metric for raising and lowering indices. We use this notation throughout this paper.

The action is invariant under a global transformation, , , and . In components, the equations of motion are

 [ h−μ0+χ(x)coshθ]ϕ(x)−A(x)ϕ∗(x)sinhθ=j(x), χ(x)/λ0=|ϕ(x)|2coshθ−s(x), A(x)/λ0=ϕ2(x)sinhθ−S(x). (9)

We note that substituting and from the last two lines of Eqs. (9) (for zero currents) into the first line of Eqs. (9) gives the equation of motion for the field with no auxiliary fieldsAndersen (2004).

Parametrizing the Green function as

 G(x,x′)=(G(x,x′)K(x,x′)K∗(x,x′)G∗(x,x′)), (10)

and using

 ∫[dx′]G−1(x,x′)G(x′,x′′)=δ(x,x′′), (11)

we obtain the equations

 [h0−μ+χ(x)coshθ]G(x,x′) (12a) −A(x)K∗(x,x′)sinhθ=δ(x,x′), [h0−μ+χ(x)coshθ]K(x,x′) (12b) −A(x)G∗(x,x′)sinhθ=0,

and the complex conjugates. Here, and are the normal and anomalous correlation functions.

Iii Auxiliary-field loop expansion

The generating functional for connected graphs is

 Z[J]=eiW[J]/ℏ=N∫DΦeiS[Φ,J]/ℏ, (13)

with given by Eq. (4). Average values of the fields are given by

 ⟨Φα(x)⟩=ℏi1Z[J]δZ[J]δJα(x)∣∣J=0=δW[J]δJα(x)∣∣J=0. (14)

If we integrate out the auxiliary fields and , we obtain the path integral for the original Lagrangian of Eq. (2). The strategy we will use here is to reverse the order of integration and first do the path integral over the fields exactly and then perform the path integration over the auxiliary fields by stationary phase to obtain a loop expansion in the auxiliary fields. Performing the path integral over the fields , we obtain

 Z[J]=N′∫DχeiSeff[χ,J]/(ϵℏ), (15)

where the effective action is given by

 Seff[χ,J] =12∬[dx][dx′]ja(x)G[χ]ab(x,x′)ja(x) +∫[dx]{χi(x)χi(x)2λ0+χi(x)Si(x) −ℏ2iTr[ln[G−1(x,x)]]}. (16)

Here is defined in Eq. (7). As shown in Ref. Bender et al., 1977, the dimensionless parameter (which we eventually set equal to one) in Eq. (15) allows us to count loops for the auxiliary-field propagators in the effective theory in analogy with which counts loops for the -propagator for the original Lagrangian. The stationary point of the effective action are defined by , i.e

 χ0(x)λ0 ={|ϕ0(x)|2+ℏRe{G(x,x)}/i}coshθ−s(x) A0(x)λ0 ={ϕ20(x)+ℏK(x,x)/i}sinhθ−S(x), (17)

where is given by

 ϕa0[χ0,J](x)=∫[dx′]G[χ0]ab(x,x′)jb(x′). (18)

Both and include self consistent fluctuations and are functionals of all the currents . Expanding the effective action about the stationary point, we find

 Seff[χ,J] =Seff[χ0,J] (19) +12∬[dx][dx′]D−1ij[χ0](x,x′) ×(χi(x)−χi0(x))(χj(x′)−χj0(x′))+⋯

where is given by the second-order derivatives

 D−1ij[χ0](x,x′) =δ2Seff[χa]δχi(x)δχj(x′)∣∣∣χ0 (20) =1λ0ηijδ(x,x′)+Πij[χ0](x,x′),

evaluated at the stationary points. Here is the polarization and is calculated in App. B. We perform the remaining gaussian path integral over the fields by saddle point methods, obtaining

 ϵW[J] =S0+Seff[χ0,J] (21) −ϵℏ2i∫[dx]Tr[ln[D−1[χi,J](x,x)]],

where is a normalization constant. From this we calculate the order corrections to and . Schematically, these one-point functions are shown in Fig. 1.

The vertex function is constructed by a Legendre tranformation (see for example Ref. Itzykson and Zuber, 1980) by

 Γ[Φ]=∫[dx]Jα(x)Φα(x)−W[J]. (22)

Here is the generator of the one-particle-irreducible (1-PI) graphs of the theoryLuttinger and Ward (1960); Baym (1962); Cornwall et al. (1974), with

 δΓ[Φ]δΦα(x)=Jα(x). (23)

Keeping only the gaussian fluctuations in , we find

 ϵΓ[Φ]=12∬[dx][dx′]ϕa(x)G−1[χ]ab(x,x′)ϕb(x′) −∫[dx]{χi(x)χi(x)2λ0−ℏ2iTr[ln[G−1[χ](x,x)]] −ϵℏ2iTr[ln[D−1[Φ]](x,x)]}+ϵΓ0+⋯, (24)

which is the negative of the classical action plus self consistent one loop corrections in the and propagators. Here, is an adjustable constant used to set the minimum of the effective potential to have finite reference energy. The effective potential is defined for static fields by

 Veff[Φ]=ϵΓ[Φ]VT=V0+12ϕaV[χ]abϕb−χiχi2λ0 (25) −ℏ2iVTTr[ln[G−1[χ](x,x)]] −ϵℏ2iVTTr[ln[D−1[Φ](x,x)]]+O(ϵ2),

where

 V[χ]ab=(χcoshθ−μ−Asinhθ−A∗sinhθχcoshθ−μ). (26)

We will see below that for the static case, and are independent of .

For a system in equilibrium at temperature , we Wick rotate the time variable to Euclidian time according to the Matsubara prescription, . Then the effective potential becomes the grand potential per unit volume, . (Details of the Matsubara formalism can be found for example in Ref. Negele and Orland, 1988.) So to leading order in , the thermal effective potential is given by

 Veff[Φ] =V0+12ϕaV[χ]abϕb−χiχi2λ0 (27) −12βVTr[ln[G−1[χ](x,x)]],

and where is a normalization constant. At the next order we have the additional term

 V(1)eff[Φ]=−ϵ2βVTr[ln[D−1[Φ](x,x)]]. (28)

Here and throughout this section, we suppress the dependence of quantities on and the thermodynamic variables . The thermodynamic effective potential is obtained by evaluating the effective potential at zero currents. From (23), this is when the fields satisfy

 δVeff[Φ0]δΦα(x)=0,for α=1,⋯,5. (29)

We call these the “gap equations” in analogy with the corresponding equations in BCS theory.

The Green functions are periodic with Matsubara frequency with , and are expanded in a Fourier series,

 G[χ](x,x′)=1β∑k,n~G[χ](k,n)ei[k⋅(r−r′)−ωn(τ−τ′)]. (30)

Writing the Green function equation in - space as

 ~G−1[χ](k,n)~G[χ](k,n)=1, (31)

we find

 ~G−1[χ](k,n) (32) =(ξk+χcoshθ−iωn−Asinhθ−A∗sinhθξk+χcoshθ+iωn),

where . So

 det[~G−1[χ](k,n)]=ω2k+ω2n. (33)

where

 ω2k=[ξk+χcoshθ]2−|A|2sinh2θ. (34)

Stable solutions are possible for . The trace-log term then becomes

 12VβTr[ln[G−1[χ](x,x)]]=12Vβ∑k,nln[ω2k+ω2n] =1V∑k{ωk2+1βln[1−e−βωk]} =∫d3k(2π)3{ωk2+1βln[1−e−βωk]}. (35)

So from (27) the effective potential to leading order in the auxiliary field loop expansion (LOAF) is given by

 (36) −12[ϕ∗2A+ϕ2A∗]sinhθ −χ2−|A|22λ0+∫d3k(2π)3{ωk2+1βln[1−e−βωk]}.

It is useful to introduce new variables and as

 χ′=χcoshθ−μ0,andA′=Asinhθ. (37)

Then the effective potential (36) becomes

 Veff[Φ′]=V0+|ϕ|2χ′−12[ϕ∗2A′+ϕ2A′∗] −(χ′+μ0)22λ0cosh2θ+|A′|22λ0sinh2θ (38) +∫d3k(2π)3{ωk2+1βln[1−e−βωk]},

where now . The gap equations (29) are now written as

 (χ′0−A′0−A′∗0χ′0)(ϕ0ϕ∗0)=0, (39) χ′0+μ0λ0cosh2θ=|ϕ|2+∫d3k(2π)3ϵk+χ′02ωk[2n(βωk)+1], A′0λ0sinh2θ=ϕ2+A′0∫d3k(2π)3[2n(βωk)+1]2ωk,

where is the Bose-Einstein particle distribution. The solutions of Eqs. (39) substituted into Eq. (38) determine the effective potential.

To calculate the finite temperature effective potential to order we need to determine . For the static case in the imaginary time formalism, and are expanded in Fourier series’ analogous to Eq. (30). So from Eq. (20) we obtain

 ~Dij[Φ](k,n)=ηijλ0+~Πij[Φ](k,n). (40)

Iv The effective potential in the condensate phase to leading order

In the language of broken symmetry, the condensate phase is a phase where the U(1) symmetry of the theory is broken since then . From Eq. (38), the minimum of the effective potential is when

 δVeff[Φ]δϕ∗∣∣ϕ0=χ′ϕ0−A′ϕ∗0=0. (41)

Because of the gauge symmetry, we can choose to be real, which means that is also real. Hence, we have the broken symmetry condition , and the dispersion relation reads

 ω2k=ϵk(ϵk+2A′), (42)

The latter is a consequence of the the Hugenholz-Pines theorem which assures that the dispersion relation does not exhibit a gap. This is equivalent to the Goldstone theorem for a dilute Bose gas with a spontaneously-broken continuos symmetry. This connection is discussed in detail in Ref. Andersen, 2004. In the absence of quantum fluctuations in , one obtains the Bogoliubov dispersion, , by setting and .

In the spontaneously broken phase, the effective potential is

 Veff[χ′] =V0−(χ′+μ0)22λ0cosh2θ+χ′22λ0sinh2θ (43) +∫d3k(2π)3{ωk2+1βln[1−e−βωk]},

where is determined by the equation

 ∂Veff[χ′]∂χ′ =χ′λ0sinh2θ−χ′+μλ0cosh2θ (44) +∫d3k(2π)3ϵk2ωk[2n(βωk)+1]=0.

These equations for and contain infinite terms that need to be regulated. In order to regulate the effective potential, we first expand in a Laurent series in

 ωk=√(ϵk+χ′)2−|A′|2=ϵk+χ′−|A′|22ϵk+⋯, (45)

around . The first three terms in the series are responsible for the divergences in the integral in Eq. (38). To regularize the theory, we subtracting these three terms from in the integrand, and replace the constant the bare interaction strength and chemical potential by regulated ones. This procedure gives the regulated effective potential

 VReff[Φ′]=VR+|ϕ|2χ′−12[ϕ∗2A′+ϕ2A′∗] −(χ′+μ)22λcosh2θ+|A′|22λsinh2θ (46) +∫d3k(2π)3{12[ωk−χ′+|A′|22ϵk]+1βln[1−e−βωk]},

which is now finite. Similarly, the regulated gap equations for and are now give as

 χ′+μλcosh2θ A′λsinh2θ =ϕ2+A′∫d3k(2π)3{2n(βωk)+12ωk−12ϵk} (47)

which are also finite.

This regularization scheme is equivalent to dimensional regularization as done for example in Ref. Papenbrock and Bertsch, 1999, or to conventional renormalization of the coupling constant and chemical potential as described in the review article of Andersen and discussed in detail for the LOAF approximation in the App. A.

V Setting the parameter θ

Up to this point, we considered a one-parameter class of mean-field approximations governed by the parameter . The dispersion relation in the condensate phase for the leading-order auxiliary field (LOAF) approximation is given by Eq. (42)

 ωk=√ϵk(ϵk+2Asinhθ), (48)

where . Now, we will choose by demanding that in the weak coupling limit, when can be ignored, the dispersion relation agrees with the one-loop low-density result obtained by Bogoliubov. Using a Hamiltonian formalism, Bogoliubov assumed

 ϕ=ϕ0+ψ, (49)

subject to the constraint . Realizing that , he then wrote the theory in terms of the classical Hamiltonian plus a quadratic fluctuation Hamiltonian, which he diagonalized. Using Eq. (49) and limiting to at most quadratic fluctuations, one has

 [(ϕ∗0+ψ∗)(ϕ0+ψ)]2→(ϕ∗0ϕ0)2+4ψ∗ψ(ϕ∗0ϕ0) +ψψ(ϕ∗0ϕ0)+ψ∗0ψ∗0(ϕ∗0ϕ0). (50)

The minimum of the classical Hamiltonian defines . One can reformulateAndersen (2004) the Bogoliubov theory in path integral language as the classical approximation plus gaussian fluctuations. The inverse Green function in the gaussian fluctuation approximation now has

 Vab[ϕ](x)=λ(2ϕ∗0ϕ0ϕ0ϕ0ϕ∗0ϕ∗02ϕ∗0ϕ0) (51)

where is defined in Eq. (5). This leads to the dispersion relation at the minimum:

 ωk=√ϵk(ϵk+2λϕ20) (52)

We will choose such that our result for reduces to the Bogoliubov dispersion relation (52) when we ignore quantum fluctuations in the anomalous density. This sets and . With our choice of , the renormalized effective potential can be written as

 VReff[Φ]=VR+χ′|ϕ|2−12(A∗ϕ2+Aϕ∗2)−(χ′+μ)24λ +|A|22λ+∫d3k(2π)3{12[ωk−ϵk−χ′+|A|22ϵk] +1βln[1−e−βωk]}, (53)

where now and

 ω2k=(ϵk+χ′)2−|A|2. (54)

The equations for the auxiliary fields are obtained from , as

 Aλ =ϕ20+A∫d3k(2π)3{[2n(βωk)+1]2ωk−12ϵk}, (55a) χ′+μ2λ =|ϕ0|2+∫d3k(2π)3{ϵk+χ′2ωk[2n(βωk)+1]−12}. (55b)

From Eq. (41) we know that at the minimum of the effective potential we have , and we can replace by the physical density using

 ρ=−∂VReff[Φ0]∂μ=χ′+μ2λ. (56)

In the broken symmetry phase we have in which case Eqs. (55) become

 χ′λ =ρ0+χ′∫d3k(2π)3{[2n(βωk)+1]2ωk−12ϵk}, (57a) ρ =ρ0+∫d3k(2π)3{ϵk+χ