Finite Temperature Description of an Interacting Bose Gas

We derive the equation of state of a Bose gas with contact interactions using relativistic quantum field theory. The calculation accounts for both thermal and quantum corrections up to 1-loop order. We work in the Hartree-Fock-Bogoliubov approximation and follow Yukalov's prescription of introducing two chemical potentials, one for the condensed phase and another one for the excited phase, to circumvent the well-known Hohenberg-Martin dilemma. As a check on the formalism, we take the non-relativistic limit and reproduce known non-relativistic results. Finally, we translate our results to the hydrodynamical, two-fluid model for finite-temperature superfluids. Our results are relevant for the phenomenology of Bose-Einstein Condensate and superfluid dark matter candidates, as well as the color-flavor locking phase of quark matter in neutron stars.

The goal of this paper is to revisit this calculation using a relativistic quantum field theory (QFT) framework. As usual, the upshot of an effective field theory approach is that it offers a systematic treatment of the relevant degrees of freedom and the symmetries governing their dynamics. Although a complete relativistic calculation of a BEC/superfluid might not be of immediate relevance to dark matter, it is nevertheless important for various other phenomena, such as the color-flavor locking superfluid phase of quark matter [23] conjectured to reside in the core of neu-tron stars. The associated observable phenomena, like pulsar glitches [24], would be sensitive to the relativistic corrections.
To properly account for the depletion of a BEC with increasing temperature, one must perform a self-consistent calculation. In QFT, this requires working with the Cornwall-Tomboulis-Jackiw (CJT) or the two particle-irreducible (2PI) framework [25]. In this approach, the effective action is computed in terms of the background field as well as the dressed propagator. The CJT formalism has also been extensively applied to thermal field theory [26][27][28]. The 2PI effective action is expressed in terms of an infinite set of diagrams having partially resummed propagators. Thus, for practical purposes, one must truncate the 2PI effective action at finite order in the loop diagram expansion. This truncation, however, gives rise to residual violations of global symmetries, which prevent the Goldstone boson of the spontaneously broken phase from being gapless. The Euler-Lagrange equations of motion of the truncated effective action are not consistent with the Ward identities. On the other hand, the complete theory must be invariant under the global U (1) symmetry, and as such must give rise to a massless Goldstone mode [29]. The inability of the CJT effective action to simultaneously satisfy the mean-field equation and display a gapless mode is similar to the Hohenberg-Martin dilemma [30] and has been noted in literature [31][32][33].
Various working solutions have been proposed to this problem depending on the particular application. It has been found that in O(2N ) theories, a Goldstone mode can be recovered in the large-N limit [32][33][34]. Other solutions add a phenomenological or an ad hoc term to satisfy both constraints [35][36][37][38]. Another way to obtain a massless Goldstone boson is to define a constrained version of the CJT formalism with the so-called external propagators [39,40] by giving up on a second-order phase transition [41]. Lastly, the symmetry-improved CJT formalism introduces new constraints on the loop-truncated CJT effective action to enforce the Ward identities [42]. While offering different insights into the nature of the problem, each of these approaches either carries a pathology, introduces ad hoc terms, or enforces constraints by hand.
For concreteness, in this paper, we follow the framework proposed by Yukalov [43]. This approach relies on introducing two different chemical potentials -one for the condensed phase, and a second one for the normal phase. The two chemical potentials are distinct below the critical temperature, and allow us to simultaneously satisfy the self-consistency condition for the meanfield while having a gapless Goldstone mode. The two chemical potentials become equal at the critical temperature. Above the critical temperature, there is of course a single chemical potential, associated with the conserved particle number. Like the other aforementioned approaches to the Hohenberg-Martin dilemma, the framework of [43] introduces an additional variable to account for the different constraints. One might argue, however, that the usage of two chemical potentials can be physically motivated, and can be naturally extended to a QFT exhibiting spontaneous symmetry breaking.
The paper is organized as follows. We first set the stage for the calculation by reviewing the imaginary time formalism and the Hartree-Fock-Bogoliubov (HFB) approximation [44] in Sec. 1.
We then introduce Yukalov's framework in Sec. 2, justifying the use of two chemical potentials. We will then work out the Matsubara summation in Sec. 3 to compute the equation of state in Sec. 4, before computing the necessary correlators and renormalizing our theory in Sec. 5. As a check on our results, we work out the non-relativistic limit of our theory in Sec. 6 and compare it to the results of [22]. As another application, in Sec. 7, we translate our results to the hydrodynamical two-fluid model of superfluidity [45][46][47], in which the system is treated as a mixture of the superfluid component and a normal fluid, made up of quasi-particles.

One-loop effective action
Similar to the grand partition function, the functional integral for a complex scalar field Φ, in the imaginary time formalism, is written as [48] where Π is the conjugate momentum, and the field Φ is periodic in imaginary time, Φ( x, 0) = Φ( x, β), with β = 1/T . 1 The free energy density is given by where V is the volume. For a complex scalar field, the Hamiltonian density H is where the potential U (Φ † , Φ) is assumed to arise from self-interactions rather than coupling to an external field. Furthermore, we assume that the predominant interactions are particle-conserving contact interactions, in which case U (Φ † , Φ) = U (|Φ| 2 ). Thus the potential U (Φ † , Φ) has a U (1) symmetry. The outline of the subsequent calculation will be valid for any such U (|Φ| 2 ) [49]. For concreteness, however, we specialize to the simple potential It is convenient to split the field into the condensate (zero-momentum modes) and the excitations: where ρ 1 , ρ 2 represent the condensate, and φ 1 , φ 2 are the excitations. Owing to the U (1) symmetry, we have freedom in choosing the condensate such that ρ 2 1 + ρ 2 2 is a constant. For simplicity, let us fix Expanding the potential in powers of the excitations, we have, at zeroth-order, The first-order terms can be ignored as usual since they do not contribute to the effective action in the one-loop approximation. The second and the fourth-order terms are respectively given by We will not be needing the third-order piece U (3) , since we will apply the HFB approximation.
The HFB method is a self-consistent mean-field approximation scheme (see, e.g., [44]). In this approach, the fourth-order terms are approximated by where the expectation values account for both thermal and quantum corrections. These HFBapproximated terms can be added to the zeroth and second-order contributions to the potential, giving us Since this HFB potential is not U (1) invariant, the symmetry of the original theory is lost. This means that the would-be Goldstone boson has become massive. This conundrum, which is of course an artifact of the field decomposition and HFB approximation, has been pointed out previously [31][32][33] and is one version of the so-called Hohenberg-Martin dilemma [30].

Two chemical potentials
We have seen that the HFB approximation gives an effective action that fails to simultaneously satisfy the self-consistent mean-field equation while having a gapless Goldstone mode. As discussed in the Introduction, several solutions have been proposed in the literature to deal with this issue. In this paper, our strategy is to introduce a second chemical potential, generalizing the prescription of Yukalov et al. in a non-relativistic setting [43]. The approach is similar to the one used by [42], in which they modify the equations of motion to ensure both constraints are satisfied.
This framework entails that we have one chemical potential µ 0 for the condensate and a separate chemical potential µ 1 for the excitations. These encode the fact that the fraction of particles in the condensate and the excitations are fixed at a given T , V and N . Mathematically, let us minimize the free energy F (N 0 , N 1 ) [50], The derivative of the free energy with respect to the number of particles is the chemical potential: µ 0 = ∂F ∂N 0 and µ 1 = ∂F ∂N 1 . Furthermore, since the total number of particles is conserved, δN 0 = −δN 1 . Thus (2.1) becomes If, in a phase transition at a particular critical temperature, the fraction of particles in either of the phases is not fixed, δN 0 = 0, then µ 0 = µ 1 . Instead, in our case the fraction of particles in the condensate and in excited states are each fixed, δN 0 = 0, therefore the two chemical potentials are no longer required to be equal. In fact, we will see that they are indeed different for T < T c and become equal at T = T c , when the condensate vanishes.
This prescription of two chemical potentials will allow us to circumvent the Hohenberg-Martin dilemma discussed previously, with each chemical potential satisfying a different constraint. The condensate chemical potential, µ 0 , enforces the mean-field equation of motion, and the one associated with the excitations, µ 1 , ensures that Goldstone's theorem is satisfied.
To incorporate two chemical potentials, the generating functional (1.1) is generalized to 2 where π i and χ i are the momenta conjugate to φ i and ρ i , respectively. The conserved charge densities for the two phases are The Hamiltonian H is given by In the above expressions, we have used that ρ i correspond to zero-momentum modes. 2 Since ρi correspond to the ground state, their Fourier transform gives the energy of the ground state which is zero. Hence, terms havingρi are ignored. Terms having ∇ρi are subsequently dropped for the same reason.
Performing the functional integrals over the conjugate momenta, fixing ρ 1 = ρ 2 = ρ as before, and applying the HFB approximation, we obtain where N is an irrelevant multiplicative constant, and the action S is given by where the HFB corrected potential terms, U HFB and U HFB , are once again given by (1.10). The Euler-Lagrange equation of motion for ρ, given by δS δρ = 0, fixes the condensate chemical potential: Notice that the effective action S and the condensate chemical potential µ 0 both depend on the expectation values φ 2 1 , φ 2 2 and φ 1 φ 2 . These expectation values are set by the excitations, and therefore vanish at T = 0 (ignoring quantum effects). To proceed, therefore, we need to compute the one-loop corrections to the free energy, which in turn will allow us to evaluate the relevant expectation values. Moving forward, we can eliminate ρ from our calculation by using (2.8). However, this substitution will make the equations more complicated. Hence, we will retain ρ in our calculation while keeping in mind that it has already been fixed using the Euler-Lagrange equation of motion.

One-loop free energy
In this Section we calculate the free energy density (1.2), including finite temperature and quantum corrections at one loop. The zeroth-order contribution can be read off from (2.7): To compute one-loop corrections, we focus on quadratic terms, written in the form The matrix elements M ij = δ 2 S δφ i (x)φ j (y) are given by Our choice of setting ρ 1 = ρ 2 in the field decomposition also leads us to expect that φ 2 1 = φ 2 2 . As a consistency check, notice that the diagonal elements of M become equal when φ 2 1 = φ 2 2 , as expected. We will henceforth make this choice.
At this stage, we decompose the fields in Fourier modes. One option is to perform this decomposition in terms of the usual ladder operators. However, working with creation and annihilation operators becomes complicated at finite temperature and requires the formalism of thermofield dynamics [51,52]. In this formalism, one introduces a dual field space, conjugate to the φ 1 , φ 2 space. Furthermore, one must go beyond the imaginary time formalism by charting a path in the complex time plane, such that dynamical effects can be captured as evolution along the real-time axis, while thermal effects correspond to evolution along the imaginary axis.
Since we are not interested in the dynamical nature of the relevant expectation values, we can avoid using thermofield dynamics altogether and instead Fourier decompose the fields as with ω n = 2πn/β. Substituting into (3.2), the quadratic action becomes where M is the momentum-space representation of M: Here we have defined It is convenient to transform to a new basis ψ ± in which M is diagonal: The transformation is unitary, and its determinant is unity. In the diagonal basis, S (2) can be written as The coefficients, ω 2 n + p k ± q 2 − 4µ 2 1 ω 2 n , are of course identified as the two eigenvalues of M. Our goal in this Section is to compute the free energy, which only depends on the eigenvalues, and not on the exact form of the transformation. We will use the explicit transformation (3.8) in Sec. 5 to compute the expectation values φ 2 1 and φ 1 φ 2 . Letting Z (1−loop) denote the contribution of the one-loop corrections to the generating functional, we obtain where the dispersion relations are given by For each term in (3.10) we can perform the Matsubara summation using the identity Ignoring the divergent constant, we obtain the one-loop corrected free energy, The first term is the zeroth-order contribution F (0) given in (3.1), the second term is the 1-loop finite-temperature correction, and the last term is the zero-point energy. This form of the free energy (3.13) is well known in prior literature, e.g., [53]. The difference with our result is that the zeroth-order contribution features a chemical potential that is different from the chemical potential entering the contribution from thermal excitations.

Thermodynamic relations
In this Section we derive various thermodynamic relations. As a first step, we can fix the excitation chemical potential µ 1 by demanding the existence of a gapless mode: lim k→0 E − (k) = 0. Using (3.11), we obtain where we have reintroduced φ 2 2 by allowing for φ 2 1 = φ 2 2 . Comparing with the condensate chemical potential µ 0 in (2.8), we see that the only difference between the two chemical potentials is due to φ 1 φ 2 being non-zero. However, setting this term to zero arbitrarily would sacrifice self-consistency. Had we started with just one chemical potential, the Euler-Lagrange equation for ρ would not be consistent with the existence of a gapless mode. This shows that the chemical potentials must in general be different in our approach.
Using the usual thermodynamic relation n = −∂F/∂µ, we obtain the number density for the condensate and excitations, As temperature increases, particles will transition from the condensed phase to the excited phase, thereby reducing n 0 while increasing n 1 , such that the total number density n = n 0 + n 1 (4.4) is conserved.
As a quick check on this result, consider the simplest example of an ideal Bose gas (λ = 0). In this case, the dispersion relations (3.11) reduce to Substituting into (4.2b) gives We recognize the difference in the number of particles and anti-particles, which is the usual result for the net charge density in a relativistic field theory. From this result, we can also see that the non-thermal part in (4.2b) arises from contact interactions between particles. Since contact interactions are an approximation to a potential that falls off with distance, we expect that our results will be UV-divergent. We will confirm later that the non-thermal part indeed diverges and must be suitably regularized.

Expectation values and renormalization
The expression for the free energy (3.13), as well as those for the chemical potentials and number densities, all depend on the expectation values φ 2 1 and φ 1 φ 2 . The goal of this Section is to calculate these expectation values. This will require the use of a renormalization scheme. Throughout the calculation we set φ 2 1 = φ 2 2 , as justified below (3.3).

Correlation functions
Our starting point is the transformation (3.8) from the original basis φ 1,2 to the diagonal basis ψ +,− . It implies φ 1,n φ 1,−n = φ 2,n φ 2,−n = 1 2 ψ +,n ψ +,−n + ψ −,n ψ −,−n ; (5.1a) where q was defined in (3.7), and we have suppressed the k dependence in the fields with the understanding that k has the same sign as n. Furthermore, note that we have dropped the cross terms ψ +,n ψ −,−n , since we are interested in calculating expectation values, and ψ + ψ − = 0 in the diagonal basis.
To compute φ 2 1 and φ 1 φ 2 , we first need to evaluate ψ 2 + and ψ 2 − , the expectation values for the diagonal basis fields. These can be expressed as a Matsubara sum over the Fourier correlators ψ 1,n ψ 1,−n and ψ 2,n ψ 2,−n , which in turn are easily determined by performing functional integrals with the generating functional using the action (3.9). This gives Let us first study φ 2 1 . From (5.1a) we see that Notice that the square root present in the individual correlators ψ 2 + and ψ 2 − , disappears in their sum. This is important because a square root in the contour integral would give rise to two branch cuts, which would prevent us from using a contour that encompasses all relevant poles. Hence, even though computing the expectation values for ψ 2 ± is rather difficult, computing their sum is straightforward since all the singularities are poles.
The Matsubara sum (5.3) can be performed through standard manipulations [48]. We first write the sums as a contour integral encompassing the poles of coth βω 2 : As shown in the left panel of Fig. 1, the contour C is oriented counter-clockwise and runs parallel on both sides of the imaginary ω-axis with infinitesimal segments crossing the imaginary axis and closing the contour at infinity on both sides. This contour encompasses all the poles of coth βω 2 , which lie at ω = iω n = 2πinT . at ω = iω n = 2πinT . Since the infinitesimal parts at the very the top and bottom of C are zero, we get the Right Panel by closing the contour with two semi-circular arcs C 1 and C 2 of infinite radius. These enclose the poles that correspond to ω = ±E ± .
Since the two infinitesimal segments at the very top and bottom of C vanish, we can transform the contour C into semi-circular arcs C 1 and C 2 of infinite radius (right panel of Fig. 1). Contour C 1 encompasses the poles at E ± , while contour C 2 encompasses −E ± . The line integrals corresponding to the infinite semi-circular arcs are zero since the integrand scales as 1/R 2 as R → ∞, where R is the radius of the arc. Thus the Matsubara sum is equal to the integral over the contours C 1 and C 2 computed using the residue theorem. This gives The expectation value φ 1 φ 2 is evaluated following similar steps. In this case we obtain The part proportional to ω n sums to zero because it is odd. The remaining sum is evaluated using the same methodology as before, with the result The two equations (5.6) and (5.8) are implicit and coupled and must be solved together at a given temperature.

Renormalization
As alluded to earlier, the correlators (5.6) and (5.8), as well as the zero-point energy term in (3.13), are all UV-divergent. We take care of these divergences using the renormalization scheme of [38].
The divergence in φ 2 1 and φ 1 φ 2 arises from coth(βE/2) = 1 + 2f B (E). The Bose factor approaches zero exponentially as k → ∞ and therefore gives a finite contribution, but the constant term is problematic. To cure this divergence, we introduce a hard momentum cut-off Λ. We first separate out the temperature-dependent term, and, from the remaining expression, we then separate out the term with µ 1 = 0. In other words, for the divergent part of the integrals for φ 2 1 and φ 1 φ 2 , denoted respectively by I φ 2 1 (µ 1 , Λ) and I φ 1 φ 2 (µ 1 , Λ), we separate out the Λ-dependent term as: I(µ 1 , Λ) = I(0, Λ) + I(µ 1 ) . (5.9) The first part of these integrals is given by where Thus the renormalized correlators are Looking at the free energy (3.13), we can similarly subtract from the zero-point energy the modified dispersion relation (ν ± ) obtained for µ 1 = 0. The result is The zeroth-order term F (0) and finite-temperature terms are both finite, thanks to the renormalized correlators. For the zero-point energy term, the ν ± subtraction removes the leading k divergence, but still leaves behind a divergent answer. One could add further counter-terms to make the result finite, as done in [38]. Instead, to parallel the analysis in [22], we will deal with the zero-point integral in the non-relativistic regime using dimensional regularization in Sec. 6.
To summarize, for a given choice of m, λ, T and n, we can now solve numerically the implicit expressions (5.12) for the renormalized correlators, (2.8) and (4.1) for the chemical potentials, as well as the requirement of charge conservation (4.4). Substituting the solution to these equations in (5.13) gives the renormalized free energy.

Non-Relativistic Limit
As a check, we will work out the non-relativistic limit of our relativistic calculation to verify that the result is consistent with the non-relativistic analysis of [22].
The non-relativistic chemical potentials µ NR 0,1 are related to their relativistic counterparts via µ NR 0,1 = µ 0,1 − m. Using (2.8) and (4.1) we obtain where we have assumed that m 2 λρ 2 , λ φ 2 1 and λ φ 1 φ 2 . Expanding the dispersion relations (3.11) for small k gives Thus we recognize ψ + as the massive mode, and ψ − as the gapless mode. Furthermore, we can read off that the gapless mode has a linear dispersion relation for low momentum, with sound speed This massless mode dominates the contribution to the excitations in the non-relativistic limit, so we henceforth ignore the contribution of the massive mode by setting ψ 2 + 0.
The condensate number density (4.2a) is approximately Meanwhile, the excitation number density (4.2b) is given by n 1 = −∂F/∂µ 1 = − δS/δµ 1 . Ignoring the massive mode, each time derivative in the action (2.7) can be replaced by iE − . Thus we obtain Therefore the non-relativistic chemical potentials (6.1) reduce to is the so-called anomalous average. Equations (6.6) agree precisely with Eqs. (42) and (44) in [22]. Figure 2: Comparison of the condensate fraction n 0 n = 1 − η between the non-relativistic case with m = 0.5 eV and n = 10 −9 eV 3 (Left Panel), and the mildly relativistic case with m = 0.5 eV and n = 4×10 −6 eV 3 (Right Panel). The dimensionless temperatures t 1 and t 2 are normalized with respect to the critical temperature of an ideal Bose gas in the non-relativistic and mildly relativistic regimes, respectively [54]. See the main text for details. Implicit expressions for n 1 and σ can be obtained by substituting the renormalized correlators (5.12). The latter also simplify in the non-relativistic limit, with the result which is consistent with Eq. (47) of [22].
These implicit equations can be solved numerically, once we specify m, λ, T and n. To do so, it is convenient to define dimensionless excitation and anomalous fractions: The condensate fraction is just n 0 n = 1 − η. Furthermore, instead of working with T and λ, we can define dimensionless temperature and interaction strength respectively as is the critical temperature for a non-relativistic ideal Bose gas. The non-relativistic condensate fraction n 0 n is plotted in the left panel of Fig. 2 as a function of t 1 and γ 1 , for m = 0.5 eV and n = 10 −9 eV 3 . Over the range 0 ≤ t 1 ≤ 1 shown in the Figure, T /m ranges from 0 to 10 −6 , which confirms the validity of the non-relativistic approximation.
For comparison, the right panel of Fig. 2 shows the condensate fraction in the mildly relativistic regime as a function of t 2 = T m 3n and γ 2 = λn 2/3 m 2 , with m = 0.5 eV and n = 4 × 10 −6 eV 3 . Note that t 1 = T /T R c , where T R c is the critical temperature for a relativistic ideal Bose gas. Over the range 0 ≤ t 2 ≤ 1 shown in the Figure, T /m ranges from 0 to 0.1, corresponding to a mildly relativistic regime. Our 1-loop effective description breaks down for T mc 2 s [59], such that we cannot reliably describe the ultra-relativistic regime. In the relativistic case, our framework breaks down for small m close to the critical temperature. Interestingly, we see from the Figure that increasing the interaction strength increases the condensate fraction significantly in the nonrelativistic case, but has no noticeable effect on the condensate fraction in the mildly relativistic case.
The expression for the renormalized free energy (5.13) also simplifies in the non-relativistic limit. In terms of the normal and the anomalous fractions (6.10), the zeroth-order contribution F (0) becomes Meanwhile, the zero-point energy contribution is given by As mentioned below (5.13), this contribution remains divergent. To parallel the analysis in [22], we evaluate the momentum integral using dimensional regularization [55]. The dominant contribution in dimensional regularization comes from the gapless mode, which yields: (1 − η + ξ) 5/2 . (6.14) Equations (6.12) and (6.14) agree with Eqs. (51) and (23) of [22], respectively.
It is instructive to consider our results at T = 0 and in the limit of a dilute Bose gas, a 3 s n 1, where a s = λ 8πm is the s-wave scattering length. As argued in [22], in the limit T → 0 the normal density n 1 goes to zero, but the anomalous average σ remains finite. In the dilute limit the integral in (6.8) can be performed explicitly, with the result Substituting into (6.12) and (6.14), we obtain the free energy at T = 0 This differs from the result of Lee and Yang [56], which ignores the contribution from the fourth-order terms. In our case, we have a non-zero anomalous average due to quantum corrections, which result from the fourth-order terms.

Hydrodynamics of a Superfluid
The existence of a BEC is related to the phenomenon of superfluidity, though there are some technical differences between the two [57]. In Landau's phenomenological model, a superfluid at finite (sub-critical) temperature behaves as a mixture of two fluids [45]: an inviscid superfluid component, and a "normal" component, which is viscous and carries entropy. In this Section we will use the results above to split the field into superfluid and normal fluid components, and derive an explicit dictionary to the hydrodynamical description.
For this purpose it is helpful to generalize the field decomposition (1.5) to Allowing the phases ψ 0 and ψ 1 to have spatial gradients enables the condensate and the excitations, respectively, to have finite velocity with respect to the frame of interest. 3 The gradient of the phases is proportional to the velocity of the superfluid in a particular frame, as we will see, and thus vanish in the rest frame of the superfluid.
Instead of implementing the chemical potentials as Lagrange multipliers, it is convenient to include them as part of the phases. Concretely, the results of the previous Sections are recovered by setting ψ α (t) = −µ α t, α = 0, 1. The Lagrangian can once again be evaluated order by order in powers of the excitations. Ignoring odd-order terms, since they do not contribute in the HFB approximation, we only concern ourselves with even-order terms: This reproduces the Hamiltonian of Sec. 1 once we set ψ α (t) = −µ α t, α = 0, 1.
In the case of a single chemical potential, there is a single conserved current [59] j µ = n s ∂ µ ψ χ + n n u µ , with χ = −∂ µ ψ∂ µ ψ, and where n s and n n are the number density for the superfluid and normal components, respectively. Meanwhile, u µ is the four-velocity of the normal component. This current satisfies the usual continuity equation ∂ µ j µ = 0.
In our approach with two chemical potentials, there are two conserved currents, given by [60] j µ α = n sα with χ α = −∂ µ ψ α ∂ µ ψ α . These reflect our demand that the charge in the condensate and excited states are individually conserved at fixed temperature. The total superfluid and normal component densities are given by n s = n s 0 + n s 1 ; n n = n n 0 + n n 1 .
Thus, in general, n s and n n receive contributions from both the condensate and the excitations.
In the normal fluid rest frame, where u µ = (1, 0, 0, 0), the currents become  α = n sα ∇ψα χα , which implies On the other hand, the conserved currents derive from the free energy via Noether's theorem, In the limit of small superflow, which is the regime of interest, we can expand the free energy as Substituting this into (7.6), we obtain where we have used χ α µ α at this order. This expression differs from the result obtained with one chemical potential [59], but agrees with it once we set µ 0 = µ 1 and ψ 0 = ψ 1 .
To compute the superfluid densities n sα explicitly, we must generalize the free energy to include spatial gradients of the phases. The dependence on ∇ψ 0 2 comes solely from the zeroth-order term F (0) . It is easy to see that (3.1) generalizes to This is recognized as the condensate density n 0 obtained in (4.2a), which tells us that the condensate only contributes to the superfluid component (i.e., n n 0 = 0).
Meanwhile, the dependence on ∇ψ 1 2 comes from the 1-loop corrections. These take the same form as in (3.13), but with modified E ± (k) to account for non-zero superflow. Specifically, we obtain (7.13) Thus it remains to calculate ∂E ± (k) . To do so, we go back to the mass matrix M. Allowing for spatial gradients of the phases, (3.6) generalizes to with k µ ≡ (iω n , k), and where p k and q are defined in (3.7). Similar to the calculation of the previous Section, the vanishing of the determinant gives the dispersion relations: where we have used ω = iω n . The solution gives the desired dispersion relations: This is an implicit relation, however, because k 0 = E ± (k) in the above. Thus a closed form solution is difficult to obtain. Fortunately, the relevant quantities, i.e, the derivatives of E ± with respect to | ∇ψ 1 |, are easy to extract: where A k ≡ 4µ 2 1 p k + µ 2 1 + q 2 , and k = k· ∇ψ 1 ∇ψ 1 . All quantities on the right-hand side are evaluated at vanishing spatial gradients, e.g., with E ± (k) given by (3.11). Substituting into (7.13) gives the excitation contribution to the superfluid density.
The total superfluid density, to leading order in the superflow, is given by the sum of (7.12) and (7.13): (7.18) This expression greatly simplifies in the non-relativistic regime, where the massive excitations can be neglected and hence f B (E + ) 0. Furthermore, in this regime it is easy to show from (7.17) that ∂ 2 E + (k) gives a suppressed contribution. Thus (7.18) becomes The different terms are to be interpreted as follows: • As already mentioned, the first term is recognized as the condensate density (4.2a), n 0 = 2µ 0 ρ 2 .
• The second term (on the first line) is independent of T and represents the contribution due to quantum corrections in the form of contact interactions at T = 0. It is easy to show that it matches the T = 0 part of the excitation density n 1 given by (6.8). 4 • Lastly the second line in (7.19) vanishes exponentially as T → 0 thanks to the Bose factors.
It follows that, in the limit T → 0, the superfluid density is equal to the sum of the condensate density plus the excitation density: where we have used (4.4). Thus the superfluid fraction is equal to unity at zero temperature, and there are no particles in the normal phase. This also agrees with the experimental observation that liquid helium, which can be modeled as having strong interactions, has a superfluid fraction close to 1, while the condensate fraction is O(10%) [61] since the excitation density increases with interaction strength. Thus, while the condensate depletes as the interaction strength between particles increases, there is no corresponding depletion of the superfluid.
Along these lines, it also follows from (7.19) that the only way for the superfluid to deplete is through the temperature-dependent terms in the second line. These terms grow with increasing temperature, resulting in a depletion of the superfluid. In particular, for sufficiently large temperature (or sufficiently weak interaction strength), the second term in (7.19) arising from quantum corrections can be neglected compared to the thermal corrections. In this case, we obtain (7.21) This matches the known result in the non-relativistic limit and for weak coupling, as shown in Eq. (66) of [22].

Outlook
The problem of describing a BEC through a scalar field exhibiting spontaneous symmetry breaking has been well-known for over five decades in the condensed matter community. After the development of the CJT formalism for studying self-consistent QFTs, this problem was again noted in terms of the inability to simultaneously satisfy the Euler-Lagrange equation of motion and Goldstone's theorem. A number of different approaches to solve this problem have been proposed over the years. Each offers different insights into the problem, but also usually suffers from a pathology or carries some undesirable baggage in the form of additional ad hoc terms or constraints.
Our own motivation for revisiting this problem is the recent interest in BEC and superfluid candidates for dark matter. This paper is a natural follow-up to our earlier work [22], where we derived the non-relativistic, finite-temperature equation of state for dark matter superfluids, using a self-consistent mean-field approximation. In this paper we extended the calculation using tells us that − µ1 ∂ 2 L ∂| ∇ψ1| 2 ∇ψ 1 =0 = 2µ0ρ 2 + µ1 φ 2 1 + φ 2 2 + terms that vanish as T → 0 . a relativistic QFT framework. As in [22], we followed Yukalov's proposal of using two chemical potentials to describe a BEC -one for the condensed phase, a second one for the normal phase. The two chemical potentials allow us to simultaneously satisfy the self-consistency condition for the mean-field while having a gapless Goldstone mode.
Our main results can be summarized as follows. We applied this proposal in the context of an imaginary time formalism QFT to describe thermal effects. We worked out the free energy of the system, incorporating the renormalization scheme of [38]. Since our calculation was done self-consistently, the resulting expressions for the condensate (excitation) and anomalous densities are implicit and can be evaluated numerically. We then worked out the non-relativistic limit and showed its consistency with the earlier results in [22]. Finally, we sought to clarify the relationship between superfluidity and BEC by translating our results to the hydrodynamical language and working out the superfluid fraction.
Though we performed an explicit calculation for a |Φ| 4 theory, our analysis can be easily generalized to any theory with a potential having |Φ| 2n terms. It would be illuminating to repeat the analysis for the more realistic superfluid effective theory proposed in [12][13][14], with hexic potential. It would be interesting, in particular, to study various observable consequences of dark matter superfluidity in our language, such as the effect of core fragmentation [16].
Even though our solution is naturally framed in a way that makes it easier to map it to the physics of a BEC, it can be easily checked against other results in the literature. Comparing with the results of [38], for instance, we find that our calculation yields the same results as the usual CJT calculation, with the only differences arising from our choice of the two chemical potentials. This choice allows us to avoid the ad hoc method used in that particular calculation, as well as others, by introducing a physically well-motivated scheme of two chemical potentials. A similar method is also used in [42], wherein a Lagrange multiplier is introduced to define a new, truncated 2PI effective action, which essentially serves the same purpose as our second chemical potential. Understanding the origin of the various approaches to this problem, as well as their similarities/differences, can help us provide deeper insights into its resolution.