3/2-Body Correlations and Coherence in Bose-Einstein Condensates

TL;DR

This paper introduces a variational wave function for weakly interacting Bose-Einstein condensates that accounts for 3/2-body processes, leading to lower energy states and a better understanding of coherence and excitation spectra.

Contribution

It presents a novel variational approach incorporating 3/2-body processes, improving energy estimates and revealing their role in coherence and excitation dynamics.

Findings

Lower energy ground state compared to mean-field theories.

Transformation of excitation spectrum into a bubbling mode.

Identification of 3/2-body processes as key to macroscopic coherence.

Abstract

We construct a variational wave function for the ground state of weakly interacting bosons that gives a lower energy than the mean-field Girardeau-Arnowitt (or Hartree-Fock-Bogoliubov) theory. This improvement is brought about by incorporating the dynamical 3/2-body processes where one of two colliding non-condensed particles drops into the condensate and vice versa. The processes are also shown to transform the one-particle excitation spectrum into a bubbling mode with a finite lifetime even in the long-wavelength limit. These 3/2-body processes, which give rise to dynamical exchange of particles between the non-condensate reservoir and condensate absent in ideal gases, are identified as a key mechanism for realizing and sustaining macroscopic coherence in Bose-Einstein condensates.

Equations

$$\hat{H}\equiv$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_1+{\hat{H}}_{3/2}+{\hat{H}}_2.$$

$${\hat{H}}_0\equiv$$

$${\hat{H}}_1\equiv$$

$$+\frac{1}{2\mathcal{V}}{\sum _{\mathbf{k}}}^{\prime }{U}_k\left({\hat{c}}_0^{†}{\hat{c}}_0^{†}{\hat{c}}_{\mathbf{k}}{\hat{c}}_{-\mathbf{k}}+{\hat{c}}_{-\mathbf{k}}^{†}{\hat{c}}_{\mathbf{k}}^{†}{\hat{c}}_0{\hat{c}}_0\right),$$

$${\hat{H}}_{3/2}\equiv$$

$${\hat{H}}_2\equiv$$

$$|n{⟩}_0\equiv \frac{({\hat{c}}_0^{†}{)}^n}{\sqrt{n!}}|0⟩(n=0,1,2,\cdots ,\mathcal{N}),$$

$$\displaystyle\hat{\beta}_{\bf 0}^{\dagger}|n\rangle_{\bf 0}\equiv|n+1\rangle_{\bf 0},\hskip 14.22636pt\hat{\beta}_{\bf 0}|n\rangle_{\bf 0}\equiv\left\{\begin{array}[]{ll}\vspace{1mm}|n-1\rangle_{\bf 0}&:n\geq 1\\ 0&:n=0\end{array}\right..$$

$${\hat{\tilde{c}}}_{\mathbf{k}}^{†}\equiv {\hat{c}}_{\mathbf{k}}^{†}{\hat{\beta }}_0,{\hat{\tilde{c}}}_{\mathbf{k}}\equiv {\hat{c}}_{\mathbf{k}}{\hat{\beta }}_0^{†}.$$

$$\displaystyle\hat{\beta}_{\bf 0}^{\nu}(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}|n\rangle_{\bf 0}=|n\rangle_{\bf 0},\hskip 14.22636pt(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}\hat{\beta}_{\bf 0}^{\nu}|n\rangle_{\bf 0}=\left\{\begin{array}[]{ll}|n\rangle_{\bf 0}&:\nu\leq n\\ 0&:\nu>n\end{array}\right.$$

$$\displaystyle\bigl{[}\hat{\tilde{c}}_{{\bf k}},\hat{\tilde{c}}_{{\bf k}^{\prime}}^{\dagger}\bigr{]}\approx\delta_{{\bf k}{\bf k}^{\prime}},\hskip 14.22636pt\bigl{[}\hat{\tilde{c}}_{{\bf k}},\hat{\tilde{c}}_{{\bf k}^{\prime}}\bigr{]}=0,$$

$${\hat{\pi }}^{†}\equiv \frac{1}{2}{\sum _{\mathbf{k}}}^{\prime }{\varphi }_{\mathbf{k}}{\hat{c}}_{\mathbf{k}}^{†}{\hat{c}}_{-\mathbf{k}}^{†},$$

$${\hat{\tilde{\pi }}}^{†}\equiv {\hat{\pi }}^{†}{\hat{\beta }}_0^2=\frac{1}{2}{\sum _{\mathbf{k}}}^{\prime }{\varphi }_{\mathbf{k}}{\hat{\tilde{c}}}_{\mathbf{k}}^{†}{\hat{\tilde{c}}}_{-\mathbf{k}}^{†},$$

$$\displaystyle\bigl{[}\hat{\tilde{c}}_{{\bf k}},\hat{\tilde{\pi}}^{\dagger}\bigr{]}=\phi_{{\bf k}}\hat{\tilde{c}}_{-{\bf k}}^{\dagger}.$$

$$|{\Phi }_{\mathrm{GA}}⟩\equiv$$

$$=$$

$${\hat{\tilde{\gamma }}}_{\mathbf{k}}|{\Phi }_{\mathrm{GA}}⟩=0(\mathbf{k}\ne 0),$$

$${\hat{\tilde{\gamma }}}_{\mathbf{k}}\equiv u_{\mathbf{k}}{\hat{\tilde{c}}}_{\mathbf{k}}-v_{\mathbf{k}}{\hat{\tilde{c}}}_{-\mathbf{k}}^{†},$$

$$u_{\mathbf{k}}\equiv \frac{1}{(1-|{\varphi }_{\mathbf{k}}{|}^2{)}^{1/2}},v_{\mathbf{k}}\equiv \frac{{\varphi }_{\mathbf{k}}}{(1-|{\varphi }_{\mathbf{k}}{|}^2{)}^{1/2}},$$

$${\hat{\tilde{c}}}_{\mathbf{k}}|{\Phi }_{\mathrm{GA}}⟩=$$

$$=$$

$$[{\hat{\tilde{\gamma }}}_{\mathbf{k}},{\hat{\tilde{\gamma }}}_{{\mathbf{k}}^{\prime }}^{†}]={\delta }_{\mathbf{k}{\mathbf{k}}^{\prime }},[{\hat{\tilde{\gamma }}}_{\mathbf{k}},{\hat{\tilde{\gamma }}}_{{\mathbf{k}}^{\prime }}]=0.$$

$${\hat{\tilde{c}}}_{\mathbf{k}}=u_{\mathbf{k}}{\hat{\tilde{\gamma }}}_{\mathbf{k}}+v_{\mathbf{k}}{\hat{\tilde{\gamma }}}_{-\mathbf{k}}^{†}.$$

$${\mathcal{E}}_{\mathrm{GA}}\equiv$$

$$=$$

$${\hat{\tilde{\pi }}}_3^{†}\equiv$$

$$|\Phi ⟩\equiv$$

$$\mathcal{E}\equiv ⟨\Phi |\hat{H}|\Phi ⟩=⟨\Phi |({\hat{H}}_0+{\hat{H}}_1+{\hat{H}}_{3/2}+{\hat{H}}_2)|\Phi ⟩$$

$${\hat{H}}_1\equiv$$

$$+\frac{1}{2\mathcal{V}}{\sum _{\mathbf{k}}}^{\prime }{U}_k\left({\hat{c}}_0^{†}{\hat{c}}_0^{†}{\hat{\beta }}_0{\hat{\beta }}_0{\hat{\tilde{c}}}_{\mathbf{k}}{\hat{\tilde{c}}}_{-\mathbf{k}}+{\hat{\tilde{c}}}_{-\mathbf{k}}^{†}{\hat{\tilde{c}}}_{\mathbf{k}}^{†}{\hat{\beta }}_0^{†}{\hat{\beta }}_0^{†}{\hat{c}}_0{\hat{c}}_0\right),$$

$${\hat{H}}_{3/2}\equiv$$

$$+{\hat{\tilde{c}}}_{{\mathbf{k}}_1}^{†}{\hat{\tilde{c}}}_{{\mathbf{k}}_2}^{†}{\hat{\tilde{c}}}_{-{\mathbf{k}}_3}{\hat{\beta }}_0^{†}{\hat{c}}_0),$$

$${\hat{H}}_2\equiv$$

$$({\hat{c}}_0^{†}{)}^n{\hat{c}}_0^m\approx {\mathcal{N}}_0^{(n+m)/2}({\hat{\beta }}_0^{†}{)}^n{\hat{\beta }}_0^m$$

$${\rho }_{\mathbf{k}}\equiv$$

Full text

I Introduction

Among the fundamental problems in the theory of Bose-Einstein condensation (BEC) are to clarify (i) how the interaction between particles changes the properties of the condensate and one-particle excitations from those of ideal gases,Huang87 ; GB68 ; ZUK77 ; Politzer96 ; TK16 and (ii) how the macroscopic coherence indispensable for superfluidity emerges. This paper makes a contribution to these issues by constructing a variational wave function for the ground state with a new ingredient, i.e., the $3/2$-body processes where a collision of two non-condensed particles throws one of them into the condensate and vice versa. These are dynamical processes beyond the scope of the mean-field treatment that have not been considered non-perturbatively. This wave function is given as a superposition in terms of the number of condensed particles within the fixed-number formalism, instead of the total number of particles in a subsystem as discussed by Anderson,Anderson66 where depleted particles serve as the particle reservoir for the condensate exchanging particles dynamically. Thus, the superposition, which is indispensable for bringing macroscopic coherence to the condensate,Anderson66 emerges naturally due to the interaction and is also maintained dynamically. The $3/2$-body processes are also shown to transform the free-particle spectrum of non-condensed particles in ideal gases into that of a bubbling mode with an intrinsic decay rate, as expected naturally in the presence of the dynamical exchange of particles between the non-condensate reservoir and condensate.

In 1959, Girardeau and ArnowittGA59 constructed a variational wave function for the ground state of homogeneous weakly interacting bosons so that it is the vacuum of Bogoliubov’s quasiparticle operators.Bogoliubov47 They used it to evaluate the ground-state energy incorporating two-body interactions of non-condensed particles, i.e., process (d) in Fig. 1, in addition to processes (a) and (b) of the Bogoliubov theory.Bogoliubov47 However, this apparent improvement brought about an unphysical energy gap in the one-particle excitation spectrum,GA59 ; Griffin96 unlike the Bogoliubov spectrum with a gapless linear dispersion and an infinite lifetime in the long-wavelength limit,Bogoliubov47 which is in contradiction to the Hugenholtz-Pines theoremHP59 or Goldstone’s theorem I.Goldstone61 ; GSW62 ; Weinberg96 This fact suggests that something crucial may be missing from the Girardeau-Arnowitt wave function, which still remains unidentified explicitly.

The key observation here is that process (c) in Fig. 1, which involves a smaller number of non-condensed particles than (d), makes no contribution to the energy in the Girardeau-Arnowitt theory.GA59 Hence, improving the variational state so as to make process (c) active is expected to lower the energy further and approximate the true ground state more closely. Such a state will be constructed below.

It is also interesting to see how process (c) affects properties of one-particle excitations. We investigate this using the moment method developed previously.TKK16 Widely accepted results on the excitations may be summarized as follows: (i) a finite repulsive interaction between particles turns the free one-particle spectrum $\propto k^{2}$ of ideal gases into the Bogoliubov spectrum $\propto k$ with a lifetime $\tau\propto k^{-5}$ that tends to infinity as the wavenumber $k$ approaches [math];Bogoliubov47 ; Beliaev58 ; GN64 ; SK74 ; WG74 ; Griffin93 ; PS02 ; PS03 ; OKSD05 ; Leggett06 (ii) the Bogoliubov mode is also identical to the density-fluctuation mode (phonons) in the two-particle channel;GN64 ; SK74 ; WG74 ; Griffin93 ; PS02 ; OKSD05 ; Leggett06 (iii) the Bogoliubov mode is the Nambu-Goldstone mode of broken $U(1)$ symmetry.Griffin93 ; PS03 ; WM13 On the other hand, an alternative picture has been presented recently based on a self-consistent perturbation expansion satisfying Goldstone’s theorem I and conservation laws simultaneously:Kita09 ; Kita10 ; Kita14 (i*′) the excitation in the one-particle channel is a bubbling mode with a finite lifetime $\tau<\infty$ even for $k\rightarrow 0$;Kita11 ; TK14 ; TKK16 (ii′) excitations in the one- and two-particle channels are different from each other;Kita10 ; Kita11 ; TKK16 (iii′*) the distinct modes in the two channels correspond to two different proofs of Goldstone’s theorem,GSW62 ; Weinberg96 the contents of which are not equivalent and should be distinguished as I and II. Indeed, it has been shownKita11 ; Kita14 that the first proof using the invariance of the effective action under a linear transformationGSW62 ; Weinberg96 is relevant to the poles of the one-particle Green’s function in the context of BEC, whereas the second (and more familiar) one based on the vacuum expectation of the commutator of the current and field is concerned with the poles of the two-particle Green’s function.

Here, it will be shown that including process (c) in Fig. 1 naturally produces the dynamical exchange of particles between the non-condensate reservoir and condensate, thereby giving rise to an intrinsic decay rate for non-condensed particles even for $k\rightarrow 0$. Thus, the result here also supports (i*′*), as does our previous study.TKK16 The finite lifetime of non-condensed particles may be regarded as a crucial element for realizing and sustaining temporal coherence in the condensate.

This paper is organized as follows. Section II constructs a variational wave function for the ground state with the $3/2$-body processes, obtains an expression for the ground-state energy, and derives the equations to determine the energy minimum. Section III presents numerical results for the ground-state energy, one-particle excitation spectrum, and superposition over the number of condensed particles obtained using our wave function. Section IV summarizes the paper. Appendix A outlines how to obtain an expression for the energy functional for the Girardeau-Arnowitt wave function. Appendix B gives a detailed derivation of the energy functional for our variational wave function with the $3/2$-body processes. Appendix C describes how to calculate the moments of the one-particle excitation spectrum.

II Formulation

II.1 System

We consider a system of ${\cal N}$ identical particles with mass $m$ and spin 0 in a box of volume ${\cal V}$ described by the HamiltonianAGD63 ; Kita-Text

[TABLE]

Here $\varepsilon_{k}\equiv\hbar^{2}k^{2}/2m$ is the kinetic energy, $(\hat{c}_{{\bf k}}^{\dagger},\hat{c}_{{\bf k}})$ are the field operators satisfying the Bose commutation relations, and $U_{q}$ is the interaction potential. We aim to describe the ground state of Eq. (1) involving BEC in the ${\bf k}={\bf 0}$ state. It is convenient for this purpose to classify $\hat{H}$ according to the number of non-condensed states involved as

[TABLE]

Each contribution on the right-hand side is given in terms of the primed sum $\displaystyle{\sum_{\bf k}}^{\prime}\equiv\sum_{\bf k}(1-\delta_{{\bf k}{\bf 0}})$ as

[TABLE]

The interactions in Eqs. (3a)-(3d) are expressible diagrammatically as (a)-(d) in Fig. 1, respectively, by symmetrizing the interaction potential as $U_{k_{1}}\rightarrow(U_{k_{1}}+U_{k_{2}})/2$ in Eq. (3c) and $U_{q}\rightarrow(U_{q}+U_{|{\bf k}+{\bf q}-{\bf k}^{\prime}|})/2$ in Eq. (3d).

II.2 Number-conserving operators

Following Girardeau and Arnowitt,GA59 ; Girardeau98 we introduce the number-conserving creation-annihilation operators as follows. First, orthonormal basis functions for the one-particle state ${\bf k}={\bf 0}$ are given by

[TABLE]

where $|0\rangle$ is the vacuum defined by (i) $\langle 0|0\rangle=1$ and (ii) $\hat{c}_{\bf k}|0\rangle=0$ for any ${\bf k}$.Kita-Text The ground state without interaction is $|{\cal N}\rangle_{\bf 0}$. Second, we introduce operators $(\hat{\beta}_{\bf 0}^{\dagger},\hat{\beta}_{\bf 0})$ byGA59 ; Girardeau98

[TABLE]

These operators are also expressible in terms of $(\hat{c}_{\bf 0}^{\dagger},\hat{c}_{\bf 0})$ as $\hat{\beta}_{\bf 0}^{\dagger}=\hat{c}_{\bf 0}^{\dagger}(1+\hat{c}_{\bf 0}^{\dagger}\hat{c}_{\bf 0})^{-1/2}$ and $\hat{\beta}_{\bf 0}=(1+\hat{c}_{\bf 0}^{\dagger}\hat{c}_{\bf 0})^{-1/2}\hat{c}_{\bf 0}$.Girardeau98 We then define the number-conserving creation-annihilation operators for ${\bf k}\neq{\bf 0}$ by

[TABLE]

Operator $\hat{\tilde{c}}_{{\bf k}}^{\dagger}$ has the physical meaning of exciting a particle from the condensate to the state ${\bf k}\neq{\bf 0}$.

It follows from Eq. (7) that

[TABLE]

holds for $\nu=1,2,\cdots$; thus, $\hat{\beta}_{\bf 0}^{\nu}(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}=1$ and $(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}\hat{\beta}_{\bf 0}^{\nu}\!\approx\!1$. For $\nu\ll{\cal N}$, the latter approximation becomes practically exact in the weak-coupling region where the ground state is composed of the kets $|n\rangle_{\bf 0}$ with $n=O({\cal N})$. This fact also implies that the operators in Eq. (8) satisfy the commutation relations of bosons almost exactly in the weak-coupling region as

[TABLE]

where $[\hat{A},\hat{B}]\equiv\hat{A}\hat{B}-\hat{B}\hat{A}$. Hereafter we replace “$\approx$” in Eq. (12) by “$=$”.

II.3 Girardeau-Arnowitt wave function

Next, we introduce the Girardeau-Arnowitt wave function in a form different from the original oneGA59 for our convenience. We define a pair operator $\hat{\pi}^{\dagger}$ with non-condensed states by

[TABLE]

satisfying

[TABLE]

Using them, we can express the Girardeau-Arnowitt wave function as

[TABLE]

where $[{\cal N}/2]$ denotes the largest integer that does not exceed ${\cal N}/2$, and ${\cal A}_{\rm GA}$ is the normalization constant determined by $\langle\Phi_{\rm GA}|\Phi_{\rm GA}\rangle=1$.

The ket of Eq. (15) is characterized by

[TABLE]

i.e., $|\Phi_{\rm GA}\rangle$ is the vacuum of the number-conserving quasiparticle operator

[TABLE]

where $u_{{\bf k}}$ and $v_{{\bf k}}$ are defined by

[TABLE]

satisfying $u_{-{\bf k}}=u_{{\bf k}}=u_{{\bf k}}^{*}$, $v_{-{\bf k}}=v_{{\bf k}}$, and $u_{{\bf k}}^{2}-|v_{{\bf k}}|^{2}=1$. To prove Eq. (16), let us operate $\hat{\tilde{c}}_{{\bf k}}$ on Eq. (15) and transform the resulting expression as

[TABLE]

where we used $\hat{\tilde{c}}_{{\bf k}}|{\cal N}\rangle_{\bf 0}=0$ and Eq. (14). Multiplying the equation in terms of $|\Phi_{\rm GA}\rangle$ by $u_{{\bf k}}$, we obtain Eq. (16).

It follows from Eqs. (12) and (18) that Eq. (17) obeys the Bose commutation relations

[TABLE]

The inverse of Eq. (17) is easily obtained as

[TABLE]

The ket presented by Girardeau and ArnowittGA59 is given by a unitary transformation on $|{\cal N}\rangle_{\bf 0}$ that appears to be different from Eq. (15). However, their equivalence can be confirmed by noting that both are (i) normalized and (ii) characterized as the vacuum of $\hat{\tilde{\gamma}}_{{\bf k}}$.

Evaluation of the ground-state energy using $|\Phi_{\rm GA}\rangle$ can be performed straightforwardly as outlined in Appendix A. Since the relevant expression is reproducible as a limit of the generalized version given below in Sect. II.4, we do not carry it out here. It may suffice to point out here that

[TABLE]

i.e., $\langle\Phi_{\rm GA}|\hat{H}_{3/2}|\Phi_{\rm GA}\rangle=0$, among the terms in Eq. (2). Neglecting the contribution of $\hat{H}_{2}$ in Eq. (21) corresponds to the Bogoliubov theory with a gapless excitation spectrum.Bogoliubov47 ; AGD63 Inclusion of the $\hat{H}_{2}$ contribution, which is supposed to improve the variational wavefunction, nevertheless gives rises to an unphysical energy gap in the excitation spectrumGA59 ; Griffin96 in contradiction to Goldstone’s theorem I.GSW62 ; Weinberg96

II.4 Including $3/2$-body correlations

Now, we improve $|\Phi_{\rm GA}\rangle$ so that $\hat{H}_{3/2}$ yields a finite contribution to lower the variational energy further below ${\cal E}_{\rm GA}$. First, we introduce an operator $\hat{\tilde{\pi}}_{3}^{\dagger}$ given by

[TABLE]

where $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}$ is a variational parameter that is symmetric in $({\bf k}_{1},{\bf k}_{2},{\bf k}_{3})$ by definition. Using Eqs. (15) and (22), we construct the following wave function:

[TABLE]

where ${\cal A}_{3}$ is determined by $\langle\Phi|\Phi\rangle=1$.

The variational ground-state energy

[TABLE]

can be estimated as follows. First, we insert either $\hat{\beta}_{\bf 0}^{\nu}(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}=1$ or $(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}\hat{\beta}_{\bf 0}^{\nu}=1$ with $\nu=1,2$ appropriately into Eqs. (3b)-(3d) to express them in terms of $(\hat{\tilde{c}}_{{\bf k}}^{\dagger},\hat{\tilde{c}}_{{\bf k}})$ as

[TABLE]

Then, we substitute the approximation

[TABLE]

in Eqs. (3a) and (25) with ${\cal N}_{\bf 0}$ denoting the number of condensed particles, and use $\hat{\beta}_{\bf 0}^{\nu}(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}=(\hat{\beta}_{\bf 0}^{\dagger})^{\nu}\hat{\beta}_{\bf 0}^{\nu}=1$ to eliminate $(\hat{\beta}_{\bf 0}^{\dagger},\hat{\beta}_{\bf 0})$ from the Hamiltonian in Eq. (24). The expectations of the remaining $(\hat{\tilde{c}}_{{\bf k}}^{\dagger},\hat{\tilde{c}}_{{\bf k}})$ operators can be calculated by performing the transformation of Eq. (20) and using Eq. (16), as detailed in Appendix B. Specifically, we obtain

[TABLE]

The density of condensed particles ${\bar{n}}_{{\bf 0}}\equiv{\cal N}_{{\bf 0}}/{\cal V}$ is expressible using the particle density ${\bar{n}}\equiv{\cal N}/{\cal V}$ and Eq. (27a) as

[TABLE]

Moreover, Eq. (24) is rewritten using Eqs. (27) and (28) as

[TABLE]

where the first term results from collecting all the contributions proportional to $U_{0}$. The Girardeau-Arnowitt functional of Eq. (21) is reproducible from Eq. (29) as

[TABLE]

II.5 Stationarity conditions

To derive the stationarity conditions of Eq. (29), we assume the symmetries

[TABLE]

in the variational parameters. Indeed, we will see that the symmetries are satisfied by the solutions. The conditions $\delta{\cal E}/\delta\phi_{{\bf k}}=0$ and $\delta{\cal E}/\delta w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=0$ for Eq. (29) can be calculated straightforwardly by performing the differentiations with the chain rule through the dependences in Eqs. (27) and (28), where

[TABLE]

also holds owing to Eq. (31). We thereby find that $\delta{\cal E}/\delta\phi_{{\bf k}}=0$ and $\delta{\cal E}/\delta w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=0$ yield

[TABLE]

[TABLE]

respectively. The quantities $\xi_{k}\equiv\delta{\cal E}/\delta\rho_{{\bf k}}$, $\Delta_{k}\equiv\delta{\cal E}/\delta F_{{\bf k}}$, and $\chi_{k}$ originating from the second line in Eq. (29) are given explicitly by

[TABLE]

By imposing $\phi_{{\bf k}}\rightarrow 0$ for $k\rightarrow\infty$, Eq. (32) can be transformed into

[TABLE]

Equations (33) and (35) with Eq. (34) form a set of self-consistent equations that can be used to determine $\phi_{{\bf k}}$ and $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}$. Equation (33) with Eqs. (34d) and (34e) indicates that $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=O({\cal N}^{-1/2})$; thus, it is more convenient for numerical calculations to rewrite the whole expressions above in terms of $\tilde{w}_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}\equiv w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}{\cal N}^{1/2}$. This procedure also enables us to confirm that the terms with $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}$ in Eq. (29) make finite contributions in the thermodynamic limit.

II.6 One-particle excitation spectrum

Now, we study one-particle excitations from Eq. (23) by calculating the first and second moments of the spectral function $A({\bf k},\varepsilon)$.TKK16 As shown in Ref. TKK16, , defining $A({\bf k},\varepsilon)$ using $(\hat{c}_{\bf k}^{\dagger},\hat{c}_{\bf k})$ necessarily shifts the excitation spectrum by the chemical potential $\mu$, as expected naturally whenever adding a particle to the system. To remove this undesirable shift, we here define the spectral function in terms of the number-conserving operators $(\hat{\tilde{c}}_{\bf k}^{\dagger},\hat{\tilde{c}}_{\bf k})$ instead of $(\hat{c}_{\bf k}^{\dagger},\hat{c}_{\bf k})$. The corresponding moments

[TABLE]

for $n=0,1,2$ can also be expressed asTKK16

[TABLE]

where $\rho_{{\bf k}}$ is defined by Eq. (27a). The mean value and width of the one-particle excitation spectrum are obtained from the moments asTKK16

[TABLE]

This $\overline{E}_{{\bf k}}$ in terms of $(\hat{\tilde{c}}_{\bf k}^{\dagger},\hat{\tilde{c}}_{\bf k})$ represents the true excitation spectrum without the chemical-potential shift, unlike the definition with $(\hat{c}_{\bf k}^{\dagger},\hat{c}_{\bf k})$. It is shown in Appendix C that Eq. (38) can be calculated straightforwardly but rather tediously. We thereby obtain the following expressions for the mean value and width of the one-particle spectrum:

[TABLE]

where we set ${\cal N}_{{\bf 0}}\approx{\cal N}$ and $A_{0}({\bf k})\approx u_{{\bf k}}^{2}$ in Eq. (39b) as justified in the weak-coupling region noting Eqs. (27a) and (28). Setting $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=0$ in Eq. (39a) reproduces the Girardeau-Arnowitt excitation spectrum,GA59 $E_{{\bf k}}^{\rm GA}=(\xi_{k}^{2}-\Delta_{k}^{2})^{1/2}$, as confirmed by using Eq. (35) with $\chi_{k}=0$ and Eqs. (18) and (27). Moreover, the Bogoliubov spectrumBogoliubov47 results from $E_{{\bf k}}^{\rm GA}$ by omitting the sums over ${\bf k}^{\prime}$ in Eqs. (34a) and (34b). Hence, our main interest in Eq. (39a) is how the presence of $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}$ changes the one-particle spectrum. On the other hand, the last term on the right-hand side of Eq. (39b) indicates that incorporating the $3/2$-body processes gives rise to a finite width $\overline{\Delta E}_{{\bf k}}>0$ even for the excitations from the mean-field wave function with $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=0$.

II.7 Superposition over the number of condensed particles

Finally, we derive an expression for the squared projection \bigl{|}_{0}\langle{\cal N}-n|\Phi\rangle\bigr{|}^{2} defined in terms of Eqs. (4) and (23), which enables us to study the superposition over the number of condensed particles in the wave function of Eq.  (23). For this purpose, we note that $\hat{\tilde{\pi}}$ and $\hat{\tilde{\pi}}_{3}$ in Eqs. (15) and (23) excite two and three particles from the condensate, respectively. With this observation, we expand the product ${\cal A}_{\rm GA}^{-2}{\cal A}_{3}^{-2}$ of Eqs. (51d) and (57) in a Taylor series and subsequently sort the terms according to the number of non-condensed particles. Multiplying the resulting expression by ${\cal A}_{\rm GA}^{2}{\cal A}_{3}^{2}$, we obtain the squared projection for $n$ excitations as

[TABLE]

where the summation is performed over all the distinct sets of $\{\ell_{1},\ell_{2},\cdots,\ell_{\nu}\}$, the quantities $I_{2\lambda}$ and $J_{3}$ are defined by

[TABLE]

and we have omitted the contribution of $J_{3\lambda^{\prime}}$ for $\lambda^{\prime}\geq 2$ as being negligible in the weak-coupling region. Equation (40) should obey the sum rule

[TABLE]

so as to be compatible with Eq. (28). Equation (42) can be used to check numerical results obtained with Eq. (40).

III Numerical Results

III.1 Model potential and numerical procedures

Numerical calculations were performed for the contact interaction potential $U_{k}=U$ used widely in the literatureLHY57 ; BS57 ; AGD63 to make a direct comparison possible. For convenience, we express this $U$ alternatively as $4\pi\hbar^{2}a_{U}/m$, i.e.,

[TABLE]

The ultraviolet divergence inherent in the potentialLHY57 ; BS57 ; AGD63 is removed by introducing a cutoff wavenumber $k_{\rm c}$ into every summation over ${\bf k}$ as

[TABLE]

The $s$-wave scattering length $a$ of this interaction potential is obtained byLL-Q

[TABLE]

which yields

[TABLE]

We choose $k_{\rm c}$ so that $k_{\rm c}a_{U}\ll 1$ is satisfied, i.e., $a\approx a_{U}$ up to the leading order.

The characteristic energy and wavenumber of this system are given by

[TABLE]

respectively. They are used to transform Eqs. (29) and (39) into the dimensionless forms ${\cal E}/{\cal N}\varepsilon_{U}$, $\bar{E}_{\bf k}/\varepsilon_{U}$, and $\overline{\Delta E}_{{\bf k}}/\varepsilon_{U}$ so that they are suitable for numerical calculations. Each sum over ${\bf k}$ in these quantities yields a factor of $(4/\pi)k_{U}a_{U}=8(2a_{U}^{3}\bar{n}/\pi)^{1/2}$, as seen by noting that $U/{\cal V}\varepsilon_{U}=1/{\cal N}$ and

[TABLE]

where $\tilde{k}\!\equiv\!k/k_{U}$. Hence, Eq. (29) in the weak-coupling region is given as a series expansion in terms of $(a_{U}^{3}\bar{n})^{1/2}\!\ll\!1$. In this context, its fourth term originating from $\hat{H}_{3/2}$ is of the same order as the last one originating from $\hat{H}_{2}$ due to the presence of $\delta_{{\bf k}_{1}+{\bf k}_{2}+{\bf k}_{3},{\bf 0}}$. Hence, process (c) of Fig. 1 yields as important a contribution as the mean-field estimation of process (d), meaning that it cannot be omitted even in the weak-coupling region.

Specifically, the energy per particle is expressible as

[TABLE]

The coefficient of $(a_{U}^{3}\bar{n})^{1/2}$ in the first expression were obtained by (i) substituting the result of the Bogoliubov approximation

[TABLE]

for Eq. (35) into the second and third terms of Eq. (29) with $\bar{n}_{{\bf 0}}\rightarrow\bar{n}$, (ii) carrying out the integrations analytically, and (iii) performing an expansion in $1/{\tilde{k}}_{\rm c}$ up to the next-to-the-leading order. On the other hand, we used Eq. (45) to derive the second expression for Eq. (48) involving $a$, which for $\tilde{k}_{\rm c}\rightarrow\infty$ reduces to the Lee-Huang-Yang expression for the ground-state energy.LHY57 ; BS57 Since Eq. (48) has a fairly strong $\tilde{k}_{\rm c}$ dependence, however, we use the first expression involving $a_{U}$ and focus on the coefficient $c_{2}$ representing corrections beyond the Bogoliubov theory, to which $\hat{H}_{3/2}$ as well as $\hat{H}_{2}$ contributes.

The sums with $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}\!\propto\!\delta_{{\bf k}_{1}+{\bf k}_{2}+{\bf k}_{3},{\bf 0}}$ in Eqs. (29) and (39) were calculated by using the transformation

[TABLE]

where we chose ${\bf k}$ along the $z$ direction, expressed ${\bf k}_{2}=(k_{2}\sin\theta_{2}\cos\varphi_{2},k_{2}\sin\theta_{2}\sin\varphi_{2},k_{2}\cos\theta_{2})$ in the polar coordinates, performed an integration over $0\!\leq\!\varphi_{2}\!\leq\!2\pi$, and made a change of variables from $\theta_{2}$ to $k_{3}=(k^{2}+k_{2}^{2}+2kk_{2}\cos\theta_{2})^{1/2}$. Integrals over $0\leq\tilde{k}\leq k_{\rm c}/k_{U}$ were calculated numerically by making a change of variables $\tilde{k}=x^{3}$ to evaluate the important region $\tilde{k}\lesssim 1$ efficiently using a small number of integration points, $N_{\rm int}\lesssim 100$.

We solved Eqs.  (33) and (35) with Eq. (34) iteratively from the initial values $\phi_{{\bf k}}=\phi_{{\bf k}}^{\rm B}$ and $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=0$ to obtain $\phi_{{\bf k}}$ and $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}$ self-consistently, where the condition $|\phi_{{\bf k}}|<1$ in Eq. (18) was incorporated by expressing $\phi_{{\bf k}}=-\cos\theta_{{\bf k}}$. The resulting solutions were substituted into Eqs. (29) and (39) to calculate ${\cal E}/{\cal N}\varepsilon_{U}$, $\bar{E}_{\bf k}/\varepsilon_{U}$, and $\overline{\Delta E}_{{\bf k}}/\varepsilon_{U}$.

Equation (40), which represents the superposition over the number of condensed particles in our wave function $|\Phi\rangle$, was calculated by omitting the contribution of $\lambda>\lambda_{\rm c}$. Choosing $\lambda_{\rm c}=10$ for ${\cal N}=20000$, $a_{U}^{3}\bar{n}=1.0\times 10^{-6}$, and $k_{\rm c}/k_{U}=10$, we verified that the sum rule of Eq. (42) was satisfied within $0.2\%$ by including $n\leq 70$.

III.2 Results

We now present our numerical results for $k_{\rm c}/k_{U}=10$, which corresponds to the cutoff energy $\varepsilon_{\rm c}=100\varepsilon_{U}$. We varied the key parameter $a_{U}^{3}\bar{n}$ between $10^{-10}$ and $10^{-5}$, where $k_{\rm c}a_{U}\ll 1$ is also satisfied so that $a\approx a_{U}$ holds in Eq. (45).

Figure 2 shows the $k$ dependence of the basic functions $\phi_{{\bf k}}$ and $\chi_{k}$, Eqs. (35) and (34c), in comparison with $\phi_{{\bf k}}^{\rm GA}$ and $\phi_{{\bf k}}^{\rm B}$ of the Girardeau-Arnowitt and Bogoliubov theories, respectively. We observe that the $3/2$-body correlations bring the basic function $\phi_{{\bf k}}$ much closer to the prediction $\phi_{{\bf k}}^{\rm B}$ of the Bogoliubov theory than $\phi_{{\bf k}}^{\rm GA}$ by making $\chi_{k}$ finite.

Figure 3 compares the coefficient $c_{2}$ in Eq. (48) with those of the Girardeau-Arnowitt theory ($\phi_{\bf k}\rightarrow\phi_{\bf k}^{\rm GA}$, $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=0$) and the Bogoliubov theory ($\phi_{\bf k}\rightarrow\phi_{\bf k}^{\rm B}$, $w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}=0$) as functions of $\log_{10}(a_{U}^{3}\bar{n})$. It shows that the ground-state energy is $20\%$ less than the estimation by the Girardeau-Arnowitt theory in the next-to-the-leading-order contribution. Although the reduction is not large, this fact clearly indicates that the $3/2$-body correlations yield the same order of a contribution as the mean-field interaction energy, meaning that they should be incorporated whenever the latter is included. In other words, the mean-field approximation for BEC is quantitatively not effective even in the weak-coupling region. Note that, in this context, the reduction of $c_{2}$ from $c_{2}^{\rm B}$ remains finite even for $a_{U}\rightarrow 0$, whereas the difference between $c_{2}^{\rm GA}$ and $c_{2}^{\rm B}$ vanishes in the weak-coupling limit.

The $3/2$-body correlations also bring a qualitative change to the one-particle excitation spectrum from the mean-field description. Figure 4 plots the mean value $\overline{E}_{\bf k}$ and the standard deviation $\overline{\Delta E}_{\bf k}$ of the one-particle excitation spectrum for $a_{U}^{3}\bar{n}=1.0\times 10^{-6}$ and $k_{\rm c}/k_{U}=10$ calculated by Eqs. (39a) and (39b), respectively. The one-particle excitation has a finite lifetime $\tau_{\bf k}=\hbar/\overline{\Delta E}_{\bf k}<\infty$ even for $k\rightarrow 0$, contrary to the predictions of the Girardeau-Arnowitt and Bogoliubov theories,Bogoliubov47 ; GA59 due to the $3/2$-body processes of Fig. 1(c). These processes also have the effect of reducing the mean value $\overline{E}_{\bf k}$, which roughly represents the peak of the excitation spectrum, from the Girardeau-Arnowitt spectrum $E_{{\bf k}}^{\rm GA}$ towards the Bogoliubov spectrum $E_{{\bf k}}^{\rm B}$. The reduction becomes larger for $k\rightarrow 0$, but $\overline{E}_{\bf k}$ finally approaches a finite value because ${\Delta E}_{\bf k}>0$ even for $k\rightarrow 0$. The finite width ${\Delta E}_{\bf k}>0$ can be traced to the dynamical exchange of particles between the non-condensate reservoir and condensate; its decrease for $k\rightarrow 0$ may be caused by the reduction of the available phase space. The inset in Fig. 4 shows the four curves over a wider range of $0\leq k\leq 3k_{U}$. The peak $\overline{E}_{{\bf k}}$, as well as $E_{{\bf k}}^{\rm GA}$ and $E_{{\bf k}}^{\rm B}$, approaches $\varepsilon_{k}$ for $k\gtrsim k_{U}$, but the width $\overline{\Delta E}_{\bf k}$ remains finite and decreases slowly for $k\gtrsim k_{U}$.

Finally, Fig. 5 shows the squared projection of Eq. (40) calculated for ${\cal N}=20000$, $a_{U}^{3}\bar{n}=1.0\times 10^{-6}$, and $k_{\rm c}/k_{U}=10$ by using Eq. (40). The superposition over the number of condensed particles has a peak at $n=24$ and becomes negligible for $n\gtrsim 60$. The profile with even integers is close to the one by the Bogoliubov approximation. On the other hand, there also is an extra contribution from odd integers due to the $3/2$-body correlations. It should be emphasized that the superposition is here realized physically and naturally due to the interaction, contrary to the case of photons with no interactions, for which SudarshanSudarshan63 and GlauberGlauber63 introduced the superposition mathematically to describe their coherence. Thus, the interaction plays a crucial role in establishing the superposition indispensable for the phase coherence within fixed-number Bose-Einstein condensates, which is also maintained here dynamically by the $3/2$-body processes.

IV Summary

We have constructed a variational wave function for the ground state of weakly interacting bosons given by Eq. (23). It incorporates the $3/2$-body processes of Fig. 1(c) to give a lower energy than the mean-field Girardeau-Arnowitt wave function, as shown in Fig. 3. This wave function is given as a superposition in terms of the number of condensed particles, as seen in Fig. 5, where non-condensed particles serve as a particle reservoir in the fixed-number formalism. Thus, the interaction naturally brings a superposition indispensable for coherenceAnderson66 to the condensate, which is sustained here temporarily by the dynamical $3/2$-body processes. The corresponding excitation spectrum is characterized by a finite lifetime even in the long-wavelength limit, as seen in Fig. 4, which reflects the dynamical exchange of particles between the non-condensate reservoir and condensate by the $3/2$-body processes. The unphysical energy gap appearing in the mean-field treatmentGA59 seems removed by the resulting broadening of the spectrum to give a finite spectral weight around $\varepsilon=0$ for $k\rightarrow 0$. However, it is still not clear whether or not Goldstone’s theorem I, which is given in terms of Green’s function for the non-condensate, is satisfied by the present treatment. This issue remains to be clarified in the future by developing a formalism to describe the one-particle excitations with the $3/2$-body correlations in terms of Green’s function.

It is widely accepted that the equilibrium in thermodynamics is realized and sustained by the exchange of momenta through collisions of particles. In contrast, little attention seems to have been paid to the origin of the coherence in Bose-Einstein condensates and superconductors. The present paper has clarified the key role played by the interaction in realizing the superposition over the number of condensed particles indispensable for coherence. An observation of the finite lifetime in the one-particle excitation spectrum, which has been predicted in the previousTKK16 and present papers, will provide a definite confirmation that the coherence is maintained dynamically.

Appendix A Calculations of ${\cal A}_{\rm GA}^{-2}$ and ${\cal E}_{\rm GA}$

The normalization constant ${\cal A}_{\rm GA}$ in Eq. (15) plays a key role in the evaluation of Eq. (21). Hence, we start by deriving its analytic expression. Imposing the condition $\langle\Phi_{\rm GA}|\Phi_{\rm GA}\rangle=1$ on Eq. (15) yields

[TABLE]

where we have replaced the upper limit $[{\cal N}/2]$ by $\infty$ to derive the second expression based on the observation that $Q_{\nu}$ for $\nu\!\sim\!{\cal N}/2$ can be set equal to zero in the weak-coupling region; see Fig. 5 regarding this point. The connected expectations $(I_{2},I_{4},I_{6},\cdots)$ have the common diagrammatic structure shown in Fig. 6(e) for $\lambda=4$ and are expressible generally as

[TABLE]

Here the factor $2^{2\lambda-1}\lambda!(\lambda-1)!$ originates from the number of combinations in connecting the $2\lambda$ pairs of field operators. Substituting Eq. (52) into Eq. (51d) yields

[TABLE]

where we used Eq. (18).

This quantity ${\cal A}_{\rm GA}^{-2}$ enables us to calculate various expectations with Eq. (15). First, $\langle\Phi_{\rm GA}|\hat{c}_{{\bf k}}^{\dagger}\hat{c}_{{\bf k}}|\Phi_{\rm GA}\rangle$ for ${\bf k}\neq 0$ can be transformed as

[TABLE]

Here we used $\hat{c}_{{\bf k}}|0\rangle=0$ and Eq. (14) for the first two equality signs, then expressed $\langle 0|\hat{\pi}^{\nu}\hat{c}_{{\bf k}}^{\dagger}\hat{c}_{-{\bf k}}^{\dagger}(\hat{\pi}^{\dagger})^{\nu-1}|0\rangle\nu$ in terms of a functional derivative of $Q_{\nu}$ in Eq. (51a) noting Eq. (13a) and that $\phi_{-{\bf k}}=\phi_{{\bf k}}$, and finally used Eqs. (53) and (18) for the last two equality signs. Note that Eq. (54) can be derived more easily by expressing $\hat{c}_{{\bf k}}^{\dagger}\hat{c}_{{\bf k}}=\hat{\tilde{c}}_{{\bf k}}^{\dagger}\hat{\tilde{c}}_{{\bf k}}$, performing the transformation of Eq. (20), and using Eqs. (16) and (19). Second, the number of condensed particles can be estimated as

[TABLE]

The second term in the final expression denotes the number of depleted particles, which can also be obtained from Eq. (54) by summing it over ${\bf k}$. Third, the expectation of $\hat{c}_{\bf 0}^{\dagger}\hat{c}_{\bf 0}^{\dagger}\hat{c}_{\bf 0}\hat{c}_{\bf 0}\approx(\hat{c}_{\bf 0}^{\dagger}\hat{c}_{\bf 0})^{2}$ is calculated as

[TABLE]

where “$\approx$” implies neglecting terms of $O({\cal N})$ compared with those of $O({\cal N}^{2})$. Equations (55) and (56) justify the procedure of replacing $\hat{c}_{\bf 0}$ by $\sqrt{{\cal N}_{\bf 0}}$ in the variational calculation using the Girardeau-Arnowitt wave function.

Appendix B Calculations of ${\cal A}_{3}^{-2}$ and ${\cal E}$

The transformation of ${\cal A}_{\rm GA}^{-2}$ in Eq. (51) is also applicable to that of ${\cal A}_{3}^{-2}$ from $\langle\Phi|\Phi\rangle=1$. Specifically, we only need to replace $(|0\rangle,\hat{\pi},[{\cal N}/2])$ in Eq. (51) by $(|\Phi_{\rm GA}\rangle,\hat{\tilde{\pi}}_{3},\infty)$. We thereby obtain

[TABLE]

This quantity is analogous to Eq. (51d) with the correspondence $(\hat{\tilde{\gamma}}_{{\bf k}},|\Phi_{\rm GA}\rangle,3\lambda)\leftrightarrow(\hat{c}_{{\bf k}},|0\rangle,2\lambda)$. Hence, we can also evaluate it analytically using Eq. (51b), the results of which can be classified diagrammatically as Fig. 7. In particular, the lowest-order contribution is obtained as

[TABLE]

It turns out that the terms of $\lambda\!\geq\!2$ in Eq. (57), which have increasing numbers of summations over ${\bf k}\!\neq\!{\bf 0}$, are negligible compared with Eq. (58) in the weak-coupling region.

Equation (57) enables us to calculate various expectations in terms of $|\Phi\rangle$ in Eq. (23). Among them, the expectations of $\hat{\tilde{c}}_{{\bf k}}^{\dagger}\hat{\tilde{c}}_{{\bf k}}$ and $\hat{\tilde{c}}_{{\bf k}}\hat{\tilde{c}}_{-{\bf k}}$ for ${\bf k}\!\neq\!{\bf 0}$ are transformed by substituting Eq. (20), using Eq. (19) to arrange the quasiparticle operators into the normal order, and noting that $\langle\Phi|\hat{\tilde{\gamma}}_{{\bf k}}\hat{\tilde{\gamma}}_{-{\bf k}}|\Phi\rangle=0$. We thereby obtain expressions for $\rho_{{\bf k}}\!\equiv\!\langle\Phi|\hat{\tilde{c}}_{{\bf k}}^{\dagger}\hat{\tilde{c}}_{{\bf k}}|\Phi\rangle$ and $F_{{\bf k}}\!\equiv\!\langle\Phi|\hat{\tilde{c}}_{{\bf k}}\hat{\tilde{c}}_{-{\bf k}}|\Phi\rangle$ as

[TABLE]

The expectation $\langle\Phi|\hat{\tilde{\gamma}}_{{\bf k}}^{\dagger}\hat{\tilde{\gamma}}_{{\bf k}}|\Phi\rangle$ in Eq. (59) can be transformed by using Eqs. (16), (19), (22), (23), (57), and (58) as

[TABLE]

Substituting Eq. (60) into Eq. (59), we obtain Eqs. (27a) and (27b).

We can also transform $W_{{\bf k}_{1}{\bf k}_{2};{\bf k}_{3}}\!\equiv\!\langle\Phi|\hat{\tilde{c}}_{-{\bf k}_{3}}^{\dagger}\hat{\tilde{c}}_{{\bf k}_{2}}\hat{\tilde{c}}_{{\bf k}_{1}}|\Phi\rangle$ for ${\bf k}_{1},{\bf k}_{2},{\bf k}_{3}\!\neq\!0$ by substituting Eq. (20), using Eq. (19) to arrange the quasiparticle operators into the normal order, and noting that $\langle\Phi|\hat{\tilde{\gamma}}_{{\bf k}_{1}}^{\dagger}\hat{\tilde{\gamma}}_{{\bf k}_{2}}\hat{\tilde{\gamma}}_{{\bf k}_{3}}|\Phi\rangle=0$ into

[TABLE]

Now, the last three lines of Eq. (60) indicate $\langle\Phi|\hat{\tilde{\gamma}}_{{\bf k}_{3}}\hat{\tilde{\gamma}}_{{\bf k}_{2}}\hat{\tilde{\gamma}}_{{\bf k}_{1}}|\Phi\rangle=w_{{\bf k}_{1}{\bf k}_{2}{\bf k}_{3}}$. Substituting this into Eq. (61), we obtain Eq. (27c).

Finally, the expectation of the operator product in Eq. (25c), which has the highest order among the terms on the right-hand side of Eq. (2), can be evaluated most easily by the Wick decomposition procedureAGD63 ; Kita-Text for $\hat{\tilde{c}}_{{\bf k}}$ within the order of the approximation we adopt. The result is expressible in terms of Eqs. (27a) and (27b) as

[TABLE]

Using Eqs. (26), (27), and (62) in the evaluation of Eq. (24) and collecting terms proportional to $U_{0}$, we obtain Eq. (29).

Appendix C Calculation of $A_{1}({\bf k})$ and $A_{2}({\bf k})$

Noting that $\hat{\tilde{c}}_{{\bf k}}=\hat{\beta}_{{\bf 0}}^{\dagger}\hat{c}_{{\bf k}}$, we express the commutator in Eqs. (37b) and (37c) as

[TABLE]

Subsequently, we substitute Eq. (2) into the right-hand side. The commutator $[\hat{\beta}_{{\bf 0}}^{\dagger},\hat{H}]$ can be evaluated by using

[TABLE]

which holds within the same order of approximation as Eq.  (26); this may be seen by replacing $|\Phi\rangle$ above by $|{\cal N}_{\bf 0}\rangle_{\bf 0}$ of Eq. (4) and calculating the commutator explicitly to the leading order in ${\cal N}_{\bf 0}$. On the other hand, $[\hat{c}_{{\bf k}},\hat{H}]$ in Eq. (63) can be calculated straightforwardly. After that, we can use the procedure for deriving Eq. (27) to evaluate Eqs. (37b) and (37c).

We first focus on Eq. (37b) and express it as a sum of the four contributions in Eq. (2) for convenience. The results for A_{1,\alpha}({\bf k})\equiv\langle\Phi|\bigl{[}\hat{\tilde{c}}_{{\bf k}},\hat{H}_{\alpha}\bigr{]}\hat{\tilde{c}}_{{\bf k}}^{\dagger}|\Phi\rangle ($\alpha=0,1,\frac{3}{2},2$) are summarized as follows:

[TABLE]

Substituting Eqs. (37a) and (65) into Eq. (38a), we obtain Eq. (39a) given in terms of Eqs. (34a) and (34b).

Calculations of A_{2,\alpha\alpha^{\prime}}({\bf k})\equiv\langle\Phi|\bigl{[}\hat{\tilde{c}}_{{\bf k}},\hat{H}_{\alpha}\bigr{]}\bigl{[}\hat{H}_{\alpha^{\prime}},\hat{\tilde{c}}_{{\bf k}}^{\dagger}]|\Phi\rangle ($\alpha,\alpha^{\prime}=0,1,\frac{3}{2},2$) can be performed similarly but rather tediously. Let us write it as

[TABLE]

It then follows that the $B_{2,\alpha\alpha^{\prime}}({\bf k})$ that contribute to Eq. (38b) are finite only for the combinations of $(\alpha,\alpha^{\prime})=(1,1),(1,\frac{3}{2}),(\frac{3}{2},\frac{3}{2})$ up to the leading order in the weak-coupling region with $B_{2,1\frac{3}{2}}({\bf k})=B_{2,\frac{3}{2}1}({\bf k})$. Moreover, it is the third term in Eq. (3b) that makes $B_{2,11}({\bf k})$ and $B_{2,1\frac{3}{2}}({\bf k})$ finite. Specifically, we obtain

[TABLE]

Substituting Eqs.  (66) and (67) into Eq. (38b) and using Eqs. (27) and (37a), we obtain Eq. (39b) up to the leading order.