Strong connection between single-particle and density excitations in Bose–Einstein condensates

Strong connection between the single-particle excitation and the collective excitation stands out as one of the features of Bose–Einstein condensates (BECs). We discuss theoretically these single-particle and density excitations of BECs focusing on the exact properties of the one-body and two-body Green’s functions developed by Gavoret and Nozières. We also investigate these excitations by using the many-body approximation theory at nonzero temperatures. First, we revisited the earlier study presented by Gavoret and Nozières, involving the subsequent results given by Nepomnyashchii and Nepomnyashchii, in terms of the matrix formalism representation. This matrix formalism is an extension of the Nambu representation for the single-particle Green’s function of BECs to discuss the density and current response functions efficiently. We describe the exact low-energy properties of the correlation functions and the vertex functions, and discuss the correspondence of the spectra between the single-particle excitation and the density excitation in the low-energy and low-momentum limits at T = 0. After deriving the exact low-energy structures of the one-body and two-body Green’s functions, we develop a many-body approximation theory of BECs with making the use of the matrix formalism for describing the single-particle Green’s function and the density response function at nonzero temperatures. We show how the peaks of the single-particle spectral function and the density response function behave with an increasing temperature. Many-body effect on the single-particle spectral function and the density response function is included within a random phase approximation, where satellite structures emerge because of beyond-mean-field effects. Criticisms are also made on recent theories casting doubt upon the conventional wisdom of the BEC: the equivalence of the dispersion relations between the single-particle excitation and the collective excitation in the low-energy and low-momentum regime.


Introduction
One of the motives for the study of the condensed matter physics is to know excitations in a quantum many-body system, which provides deep understandings of physics behind the system [1][2][3][4]. Various kinds of response functions, such as the single-particle spectral function, density response function, pair-correlation function and spin response function, are useful to understand excitations including the single-particle excitation, density excitation, pair-breaking, and spin excitation. Generally, even if the two-body correlation function is constructed from the one-body correlation function, the peak structure of the single-particle excitation does not directly clearly emerge as the exact same peak structure in the density response function generated from the two-body correlation function, where effect of the single-particle excitation may emerge as the broad continuum [1][2][3][4]. In contrast to this wisdom, Bose-Einstein condensates (BECs) is of particular interest, since the single-particle property strongly relates to the collective property [5]. The Josephson sum-rule concludes that the outcome of the coherent flow of particles is related to the single-particle spectral function, which explicitly gives the relation among the superfluid density, the condensate density, and the single-particle Green's function [6][7][8]. Gavoret and the Nepomnyashchii-Nepomnyashchii identity [38]. Section 6 reviews the earlier experimental and theoretical studies focused on the sound mode in the superfluid helium as well as ultracold atomic BECs, where variant sound modes in the superfluid, such as the second sound, are important but beyond the scope of the present paper. This section also serves as criticisms of recent theories casting doubt upon the paradigm about the BEC: the equivalence of the dispersion relations between the single-particle excitation and the collective excitation in the low-energy and low-momentum regime. Section 7 develops the formulation of the random phase approximation (RPA) in terms of the matrix formalism. Section 8 discusses the single-particle spectral function and the density response function at nonzero temperatures within the many-body approximation developed in the previous section 7. This section also addresses the correspondence between peaks of the single-particle spectral function and the density response function, and studies the sound speed estimated from the compressibility zero-frequency sum-rule by using the density response function obtained in the RPA. We end with the summary and conclusions in section 9.
Throughout this paper, we set = k B = 1, and take the system volume V to be unity. The terms, n-particle irreducible (nPI) and n-particle reducible (nPR), are applied to represent diagrams that cannot and can be separated into two pieces by cutting n single-particle lines, respectively. The regular part called in this paper means the proper part, which represents a diagram that cannot be separated into two pieces by cutting a single interaction line.

Correlation functions
We consider the Hamiltonian of an interacting Bose system with the atomic mass m, given bŷ H = dr 1 2m ∇Ψ † (r)∇Ψ(r) whereΨ(r) andΨ † (r) are bosonic annihilator and creator, respectively. In the BEC ordered phase, the field operator may be treated by the so-called Bogoliubov prescription: where Φ 0 (r) represents the order parameter of the condensate wave function, andφ(r) the non-condensate part. In the uniform system with the condensate density n 0 , we may supposê where the condensate wave function is taken to be real, i.e., Φ 0 (r) = √ n 0 . We consider a contact interaction U int (r − r ) = Uδ(r − r ), where the interaction strength U is related to the s-wave scattering length a s through the relation 4πa s /m = U/[1 + U p c p 1/(2ε p )], where p c is the cutoff momentum, and ε p the kinetic energy of the bosonic particle ε p = p 2 /2m.
An average of an operatorÔ at temperature T is given by Ô = Tr[exp(−Ĥ/T)Ô]/Ξ. Here, the HamiltonianĤ with the chemical potential μ is given byĤ =Ĥ − μ drφ † (r)φ(r), where the Bogoliubov prescription is applied toĤ. The partition function is given by Ξ = Tr[exp(−Ĥ/T)], which may be regarded as the quasi-grand partition function because the term −μn 0 is omitted from the HamiltonianĤ. It is sufficient to evaluate an average by using exp(−Ĥ/T) with the Bogoliubov prescription, because the term −μn 0 is the c-number, which is reduced in the form of the average.
The Dyson equations for the Green's functions are given by G(p) = G 0 (p) + G 0 (p)Σ(p)G(p), with p = (iω n , p). Each matrix equation provides equivalent equations for the Green's function G ij with i, j = 1, 2. Here, G −1 0 (p) = iω n σ 3 − ε p + μ is the (2 × 2)-matrix Green's function for non-interacting bosons, where σ j=0, 1,2,3 are the Pauli matrices. The identity matrix of the size 2 is given by σ 0 . The non-interacting parts G 0 and G † 0 are given by where G (0) (p) = 1/(iω n − ε p + μ). Interaction effects are included into the (2 × 2)-matrix self-energy Σ, and we may introduce the (4 × 1)-matrix self-energy Σ, and the (1 × 4)-matrix self-energy Σ † . Diagrammatic representations of matrix elements of (G, G, G † ) as well as (Σ, Σ, Σ † ) are summarized in section 3. Matrix elements G ij are not independent of each other: G 11 (p) = G 22 (−p) and G 12 (p) = G 21 (p) = G 12 (−p) = G 21 (−p), where the self-energies Σ ij satisfy the same relations [1,9]. The Green's function G provides the non-condensate density where n is the total particle density. In the following, we omit the convergence factor exp(σ 3 iω n δ) for simplicity. The formalism at T = 0 is introduced by applying the analytic continuation (iω n → ω + iδ) as well as the following replacement −T n → i dω/(2π) [1,33]. The (4 × 4)-matrix two-particle Green's function K(p, p ; q) is composed of the one-particle reducible (1PR) and one-particle irreducible (1PI) parts, i.e., K 1PR and K 1PI , where the 1PR part is specific to the condensed Bose system. (In reference [9], these are called the singular and regular parts, respectively.) The two-particle Green's function K is given by (see figure 1) where For the 1PI part, K 0 is a bare part of the (4 × 4)-matrix two-particle Green's function and Γ is the (4 × 4)-matrix four point vertex, given by (see figure 2(a))  Here, I(p, p ; q) is a two-particle irreducible (2PI) part of the (4 × 4)-matrix four point vertex. The matrixX is given byX which exchanges upper and lower ends of a two-particle Green's function (see section 3). For the 1PR part, the (4 × 2)-and (2 × 4)-matrix vertices Q and Q † are given by Here, the (4 × 2)-and (2 × 4)-matrix three point vertices P and P † are given by (see figures 2(b) and (c)) where J and J † are 2PI parts of P and P † . The condensate contributions here are included by the (4 × 2) and (2 × 4)-matrix condensate Green's functions where G 1/2 and G † 1/2 are the condensate Green's functions In the case at T = 0, the factor √ −1 is replaced with √ −i [40]. The density and current correlation functions χ μν (q) are defined as where the density-density and current-current correlation functions are χ 00 (q) and χ ij (q) for i, j = 1, 2, 3, respectively. Here, i, j = 1, 2, 3 are the index of the Cartesian coordinate. The density-current correlation functions are χ 0i (q) and χ i0 (q) for i = 1, 2, 3. The density and current vertex vector |λ μ (p; q) is given by where and The density vertex vector is simply given by |λ 0 (p; q) = |f 0 . The correlation functions (24) are constructed from the two-particle Green's function K, which are decomposed into the 1PI and 1PR parts, giving the form The 1PI and 1PR parts are of the form (see figure 3) where we used relationsX 2 = 1,X|λ ν (−p − q; q) = |λ ν (p; q) andXK 0 (−p − q; q)X = K 0 (p; q), and introduced the density and current vertices

Diagrammatic representations and matrix forms of correlation and vertex functions
Correlation and vertex functions are presented in the matrix form in this paper. It is very convenient to explicitly provide these representations in terms of diagrams. We apply the following rules to satisfy the conservation law. The point with the filled circle (•) connects to an outgoing external particle line ( figure 4(a)). The point with the open circle (•) connects to an incoming external particle line ( figure 4(b)). The point (•) can also connect to the point (•) and vice versa. The (2 × 1)-matrix vertex functions include the condensate Green's function G 1/2 , and the two point vertex Υ ν that connects to an external particle line and an external interaction line (figure 5). The (1 × 2)-matrix vertex functions G † 1/2 and Υ † μ are their counterparts (figure 6). The (2 × 2)-matrix correlation and vertex functions include the single-particle Green's functions G 0 and G, as well as the self-energy Σ (figure 7). These functions are also given in the (4 × 1)-or (1 × 4)-matrices. The (4 × 1)-matrix correlation and vertex functions include the Green's functions G 0 and G, the self-energy Σ, as well as the three-point vertex γ that connects to two external particle lines and an external interaction     The (4 × 2)-matrix vertex functions include the condensate Green's functions G 1/2 (= σ 0 ⊗ G 1/2 ) and XG 1/2 (= G 1/2 ⊗ σ 0 ), as well as the three point vertex P that connects three external particle lines (figure 10). The (2 × 4)-matrix vertex functions are their counterparts such as G † 1/2 (= σ 0 ⊗ G † 1/2 ), G † 1/2X (= G † 1/2 ⊗ σ 0 ) as well as P † (figure 11). The (4 × 4)-matrix correlation and vertex functions include the four-point vertex Γ, as well as the two-particle Green's functions, such as K, K 0 ,XK 0 and K 0X ( figure 12 and 13).
The matrixX provides exchange contributions of the two-particle Green's function. In the diagram of XK 0 (K 0X ), upper and lower left (right) ends of the two-particle Green's function K 0 are exchanged Three point vertex P that connects to three external particle lines. The matrix Q as well as the 2PI part J of the three point vertex P are also given in the same matrix form.
(c) Three point vertex P † that connects to three external particle lines. The matrix Q † as well as the 2PI part J † of the three point vertex P † are also given in the same matrix form.

Relations between vertex functions
Vertex functions in the static and zero-momentum limits can be systematically generated from all the possible linked diagrams that construct the thermodynamic potential Ω = −T ln Ξ. An exact many-line vertex M(n out , n in , n U ) is given by [36] M(n out , n in , n U ) = n Here, n in(out) is the number of incoming (outgoing) particle lines that can connect to the vertex function M, and n U is the number of external interaction lines U that can also connect to the vertex function M. The operator n n in(out) /2 0 (∂/∂n 0 ) n in(out) works as the elimination of the n in(out) condensate lines Φ ( * ) 0 (= √ n 0 ) from the linked diagrams. The operator (−∂/∂μ) n U affects on the Green's function G 0 in the linked diagrams, which provides the n U vertex points for the interaction line due to the relation −∂G 0 /∂μ = G 2 0 . This prescription was originally invented for the ground state energy at T = 0 [36]. Since the liked diagrammatic structures for the thermodynamic potential at nonzero temperature are formally the same as those of the ground state energy [41], this prescription is also applied to the nonzero temperature case [42].
The equation (33) generates the self-energy matrix Σ, the three point vertex matrix P, and the density vertex Υ 0 in the zero-energy and zero-momentum limits, respectively given by where |± = (1, ±1) T , ±| = (1, ±1) and Vertices Σ and P are related with each other, giving the form where Â = diag(1, 0, 0, −1). The thermodynamic potential Ω is related to the grand potential Ω by introducing the chemical potential of the condensate μ 0 . We have the relation Ω = Ω − μ 0 n 0 with the condition μ 0 = μ [36]. Since the condensate density is determined from the condition ∂Ω/∂n 0 = 0, we have Given this relation, we may derive the Hugenholtz-Pines relation Σ 11 (0) − Σ 12 (0) = μ [43] (or its matrix form Σ(0)|− = μ|− ). The Nepomnyashchii-Nepomnyashchii identity [35,36], giving the form reduces the Hugenholtz-Pines relation to the following form The derivation of the Nepomnyashchii-Nepomnyashchii identity with the use of the matrix formalism is summarized in reference [38], and physics of this identity can be found in references [5,35,36,38,39,[44][45][46][47][48][49][50]. The Nepomnyashchii-Nepomnyashchii identity Σ 12 (0) = 0, which is strongly related to the weak infrared divergence of the longitudinal susceptibility caused by the convolution of the phase-phase correlation function [38,44], can be obtained from the relation between vertex functions and the nature of the infrared divergence in the self-energy diagrams [35,36,38]. As a result, this identity is also valid at nonzero temperature [38,42]. The identity (40) also provides a relation P(0; 0)|− = 0. The Nepomnyashchii-Nepomnyashchii identity also provides the zero-frequency density vertex identity, i.e., the vanishing density vertex in the limit p = 0, giving the form This is valid in the isothermal condition, and the derivation is summarized in appendix A.
In the remaining part of this section, we summarize low-energy behaviors of vertex functions in our matrix representation. We first consider a relation between Σ and P as well as a relation between G and L, where the (4 × 2)-matrix L (diagrammatically described in figure 14) is given by With respect to (Σ, P) or (G, L), relations at small but finite q = ( , q) are given by [9] P(p, +q)|0 where |0 = (1, 0) T , |1 = (0, 1) T andD(q) = D 1 (q) +BD 2 (q) withB = diag(0, −1, +1, 0) and Here, two types of derivatives are introduced: δ/δx ν = (∂/∂μ, δ/δp) and ∂ ν = (∂/∂ω, ∂/∂p). The partial derivative ∂/∂p i and the total derivative δ/δp i for i = 1, 2, 3 are respectively defined as [G (0) (±(p + dp i e i )) − G (0) (±p)]/dp i , and where e i is the unit vector in the Cartesian coordinate [9]. The total derivative δ/δp i is related to an observation of the system from a reference frame with a speed −δp/m. (Details can be found in section 5 in reference [9].) In the limit q = 0, (44) and (45) are reduced into The three point vertex P (or L) is created from the two-point vertex Σ (or the Green's function G) by eliminating a condensate line from Σ (or G) that provides an extra vertex point at q = 0 [9].
We also have relations [9] The upper equality in (49) can be derived as follows [9]; the self-energy Σ can be constructed from two parts. One is the three point vertex J, where one of the three vertex points is blocked by a condensate line √ n 0 . The other is the four point vertex I, where two of the four vertex points are blocked by a Green's function G 12 or G 21 . It gives the following relation (see appendix B in reference [9]): By comparing equations (20) with (50) with the use of a mathematical identity ÂG(p) = K 0 (p; 0)ÂΣ(p), we obtain the first equality of (49), which is consistent with (38) at p = 0. The second equality in (49) with respect to (L, G) is also obtained by following the similar way. According to symmetries, the density and current vertices can be given by the matrix element of which in the low-energy regime behaves as (see appendix C)
The first equality indicates that the differential ∂ ω is related to ∂ μ , since the self-energy is constructed from the non-interacting Green's function G (0) (p) = 1/(ω − ε p + μ) and then the infinitesimally small increase of the energy ω + δω in the self-energy Σ can be regarded as the infinitesimally small increase of the chemical potential μ + δω in a Green's function G (0) that constructs Σ [9]. With respect to the second order of p or ω, similar relations are obtained, giving the forms [9] where c T is the isothermal sound speed (see appendix A). Note that although the relations between the vertex functions as shown in section 4 hold at nonzero temperatures because the diagrammatic structure is the same as in the case at T = 0, the relations between the thermodynamic quantities and the differentiations of vertex functions with respect to p shown here may not hold at nonzero temperatures. The single-particle Green's function in the low-energy regime is reduced into [36] (see appendix B) The first term in (57) is the leading term of G, which provides the phonon spectrum of the single-particle excitation, whose sound speed c T corresponds to the isothermal sound speed related to the compressibility (A7). This first term is important to the low-energy behavior of the density and current correlation functions, and essential for the transverse susceptibility G ⊥ (p) = − −|G(p)|− /4 with respect to the BEC order parameter [38,47]. The second term in (57) also provides the infrared divergence, because of the Nepomnyashchii-Nepomnyashchii identity, where the infrared divergence of the second term is much weaker than that of the first term in (57). This weak infrared divergent second term never plays an essential role in the density and current correlation functions. However, it plays an important role in the longitudinal susceptibility G (p) = − +|G(p)|+ /4 [38, 45-47, 49, 52]. The transverse and longitudinal fluctuation operators are not commutable [46,47]. Since the transverse fluctuation is regarded as the phase fluctuation, the longitudinal fluctuation might be expected to represent the amplitude (Higgs) mode. However, the longitudinal susceptibility does not describe the gapped amplitude mode, and shows the weak infrared divergence in the low-energy limit, because it is provided from the convolution of the phase-phase correlation function. The response function that can capture the Higgs mode is the scalar susceptibility [53].
In the low-energy regime, the 1PI part χ 1PI μν (q) can be given by [9] (see appendix C) The density-density correlation function remains nonzero in the low-energy and low-momentum limit. On the other hand, the density-current and current-current correlation functions vanish in the same limit.
Using (52), we have the following simple expression of the density and current vertex in the low-energy and low-momentum limits: The density and current vertex vanishes in the static and low-momentum limits, i.e., Υ ν (0) = 0, which is consistent with the exact identity (42). To obtain (62), we employed the identity (A2) in the isothermal condition for ν = 0, which is the consequence of the Nepomnyashchii-Nepomnyashchii identity, and employed the relation The 1PR parts of correlation functions are then of the forms which provide where c T is the isothermal sound speed that can be found in the single-particle Green's function (57). Here, we used a relation which is conveniently obtained from (57). From equations (57) and (69), the effect of the Nepomnyashchii-Nepomnyashchii identity and the infrared divergence of the longitudinal susceptibility G = − +|G|+ /4 = 1/4Σ 12 (0), which comes from the second term of (57), are clearly found to be irrelevant to the density and current correlation functions. The correlation functions can be obtained from the 1PI parts (58)-(61) and the 1PR parts (66)-(68), given by [9] χ 00 (p) n m The density-density correlation functions (70) satisfies the compressibility zero-frequency sum-rule, giving the form [5,40] The compressibility zero-frequency sum-rule is exhausted by the 1PI part. In the low-energy limit, the 1PI and 1PR parts of the density-density correlation function in (58) and (66) The 1PR part vanishes in the static and low-momentum limits, and does not contribute to the compressibility zero-frequency sum-rule.
The leading term of the single-particle Green's function (57) and the density and current correlation functions (70)- (72) share the pole, which provides the phonon excitations. Because of the presence of the BEC, the two-particle Green's function involves the single-particle Green's function as in the 1PR part (14). This contribution directly involves the single-particle property to the density correlation function. Since the self-energy in the single-particle Green's function can be related to thermodynamic quantities as discussed in this section, which can provides the phonon dispersion relation with the isothermal sound speed, the density correlation function can consistently describe the sound mode, the speed of which is defined in terms of the macroscopic compressibility. The paper by Huang and Klein [51] also provides a useful discussion about the phonon mode in BEC.
The single-particle excitation is also related to the superfluidity. An interesting relation between them owes to the Josephson sum-rule [6,7], given by where ρ s is the superfluid mass density. By using equation (57) as well as the relation ΔΣ(p) p 2 in the small momentum regime, we find the relation which indicates that the superfluid mass density is exactly the total mass density at T = 0. The current-current response function can be decomposed into the longitudinal and transverse response functions, given by [8,10,54] These longitudinal and transverse response functions are extracted from the relations [8] The longitudinal response function satisfies the f-sum rule χ L (p → 0, 0) = −n/m and the transverse response function provides the normal fluid density n n , given by χ T (p → 0, 0) = −n n /m [3,8,10,54]. As a result, the superfluid mass density can be given by ρ s = m 2 lim p→0 [χ T (p, 0) − χ L (p, 0)]. Using the results (61) and (65), we obtain which is consistent with the f-sum rule as well as with the fact that at T = 0, the normal fluid density is absent and the superfluid mass density is equal to the total mass density as in equation (76).

Experimental and theoretical studies of sound modes in superfluid
This section presents an overview of the experimental and theoretical studies of excitations in superfluid 4 He and in BECs of ultracold atomic gases, which will be helpful to bridge both fields and to push further the study of the single-particle and collective excitations in BECs in ultracold atoms. The static structure and dynamic structure have been intensively and extensively studied on the superfluid liquid 4 He experimentally [55][56][57][58][59][60][61]. The dynamic structure factor S(q, ω) consists of a sharp peak superimposed on a broad background in the superfluid 4 He [42,55,[58][59][60][62][63][64][65][66][67]. The sharp peak in S(q, ω) is interpreted as a collective density mode as well as a (single)-quasiparticle excitation arising from the 1PR part of the density response function χ 1PR 00 [59,60,65,66], where the density and single-particle responses have the same pole [60,66,68]. The broad component is interpreted as the multi-particle excitation originated from the 1PI part of the density response function χ 1PI 00 [60,63,67]. The temperature dependence of S(q, ω) is quite different above and below T c [59], and abruptly changes at T c [60]. As the temperature increased, the sharp peak broadens [60,65,66], which is well described by quasiparticle-quasiparticle scattering [60], and it loses intensity [58-60, 65, 66], since the condensate density decreases, which includes the single-particle Green's function to the density response function.
At low momentum regime, the superfluid has a single phonon mode [60,65,69], whose peak is very sharp at low T, where the phase-space for the decay of a single phonon into two is limited [59]. The sharp peak at the maxon momentum region is also interpreted as a contribution from a quasiparticle excitation [65]. In the high momentum regime, the superfluid 4 He does not support a collective density mode [60], and the density response function in this momentum regime broadens in the normal and superfluid phases [60,65,70].
The broad component is considered as multi-quasiparticle excitations with the high-energy tail, which originates from roton-roton, maxon-maxon, and maxon-roton contributions [70,71]. This broad continuum does not contain a collective mode in the superfluid phase [60], which starts from a finite positive energy [62]. The broad multiphonon component and high-frequency tail are largely temperature independent [59,60,72].
Above the critical temperature, the sharp peak phonon-maxon-roton excitation disappears from S(q, ω) [58-60, 65, 66], where the single-particle Green's function does not contribute to the density response function [59]. In more detail, the sharp component disappears in the maxon and roton momentum regions, but the peak remains well defined in the low-momentum phonon region, which indicates the existence of a collective density mode [60,65].
In the theoretical framework, it can be clearly seen that the condensate plays an essential role in coupling the density excitation and the quasiparticle excitation [60,65,66,[83][84][85], where this hybridization disappears above the critical temperature [83]. In the low momentum phonon regime, the single-particle Green's function and density response function share the pole [9,84]. Above the critical temperature, where the hybridization is absent, the maxon-roton peak vanishes in S(q, ω), which suggests that the sharp maxon-roton intensity originates from the single-particle excitation and the BEC in the superfluid 4 He [82].
For the hybridization, the dielectric formalism [66,82,86,87] is an approach that fulfills the Ward identities related to the conservation of particle number and the breaking of the gauge symmetry, i.e., a conserving and gapless approach by using the continuity equation [87]. It gives the same pole in the single-particle Green's function and the density correlation function in the superfluid phase [66,87], and the density fluctuation is coupled into the single-particle excitation though the condensate [86].
The collective description is a theory described by the canonical collective variables, i.e., the density fluctuation and velocity operators [90-92, 103, 106-110], which is a divergent free approach [90][91][92]. The collective description is employed to study energy spectrum [90][91][92]109], focusing on effects of the phonon-phonon interaction [90], and phonon-roton interaction [92], which play an important role in the phonon velocity and roton minimum, and is employed to study the temperature dependence of the sound velocity and the absorption coefficient including the thermal roton effect [103].
The kinetic equation is also applied to study the sound velocity and absorption [98,99,102,111]. In this approach, collisions between excitations are assumed to be not frequent in the superfluid helium at low temperatures, and thus the kinetic equation in the collisionless regime is employed. The sound velocity in the sufficiently low temperatures increases as T 4 ln(const./T) [99], where the constant is very small [102], and the absorption is reported to follow the T 6 -law [99].
Since successful creation of the BEC in alkali atom gases [112,113], the condensate excitation in ultracold gases has been intensively and extensively studied [15,16,114]. Sudden modification of the trapping potential can create the local density fluctuation, and the dynamical propagation of the density fluctuation has been measured by using the phase-contrast images, where the propagation speed is consistent with the Bogoliubov theory [20]. Two-photon Bragg scattering is a useful tool to study the excitation in the BEC of ultracold gases [21]. The Bragg spectroscopy has been applied to measure the structure factor of the BEC in the phonon regime, the line shift and line strength of which are consistent with the results of the local density approximation [22]. The Bragg pulses have also been applied to observe the Bogoliubov transformation for a BEC [115,116], and to reveal the wide range of the excitation spectrum from the phonon regime to the single-particle regime, which is also consistent with the Bogoliubov theory with the local density approximation [23]. By using the Bragg spectroscopy, experiments have probed the excitation in a strongly interacting BEC [24], as well as the roton-type excitation in BECs with cavity-mediated long-range interactions [117], with spin-orbit couplings [118], in shaken optical lattices [119], and with dipole interactions [120]. Recently, the sound propagation of the BEC trapped in a box trap has been intensively and extensively studied, including a uniform two-dimensional Bose gas [121], and a cylindrical box trap with tuning the atomic density [122], which are free from the conventional restriction of the harmonic trap potential.
Through the development of the study on BECs in ultracold atoms, theories have been proposed [27][28][29][30] that cast doubt upon the conventional wisdom about the BEC, where those recent theories claim that the dispersion relation of the single-particle excitation is not phonon and not equal to that of the collective excitation in the low-energy and low-momentum regime, which contradicts the earlier result given by Gavoret and Nozières [9]. It is concluded from two different approaches: the Luttinger-Ward thermodynamic functional approach (Φ-derivable approximation) [27][28][29]123] and a functional renormalization group approach [30,124].
The Luttinger-Ward thermodynamic functional approach [27][28][29]123] is useful for considering the theory satisfying the Noether's theorem and the Goldstone's theorem, which may cure the so-called conserving-gapless dilemma [40,[125][126][127]. The papers [28,29] are concluded that the self-energy contribution should be one-particle reducible (1PR), because the 1PR contribution cures the conserving-gapless dilemma. As a result, the two-particles Green's function has the pole showing the collective sound mode; on the other hand, the single-particle Green's function provides a bubbling mode with a considerable decay rate rather than the sound mode, which results in no well-defined quasiparticle in BECs [28]. However, in general, in the case where the self-energy contribution is included to the Green's function through the Dyson-Beliaev equation, the one-particle irreducible part should be employed. Otherwise, multi-counting of diagrammatic contribution emerges in the full Green's function. In this regard, even if the 1PR approximation may avoid the conserving-gapless dilemma, it provides a problem, namely, the trilemma among conserving, gapless and 1PR approximation in the BEC theory.
Other approach, the extension of Bijl-Feynman formula, also provides the result that the energy spectrum of the single-particle excitation is distinct from that of the collective excitation, and the lifetime of the quasiparticles remains finite even in the long-wavelength limit [31]. In this respect, the Josephson sum-rule (75) and the Bogoliubov operator inequality [54,126,[130][131][132] could be useful criteria for the result contradictory to the conventional wisdom about the single-particle excitation and the collective excitation in BECs.

Density response function in random phase approximation
The matrix formalism is a useful tool to develop many-body theories, such as the RPA, for studying many-body effects as well as the density-density correlation function in BECs. The same idea of the matrix formalism for the BEC may be found in the study of an effective roton-maxon interaction in liquid He II [81]. In the BEC phase, the density-density correlation function is constructed from the sum of the 1PI and 1PR parts as in (28), i.e., χ 00 = χ 1PI 00 + χ 1PR 00 . In the following, we omit the subscript describing the density vertex μ = ν = 0 in the polarization function χ as well as the density vertices Υ and Υ † , for simplicity.
We consider the 2PI parts I(p, p ; q), J(p; q) and J † (p; q) introduced in section 2 as the simplest contributions, given by where the first and second terms inÛ provide the Hatree and Fock contributions in the present matrix formalism, respectively. By assuming the momentum and frequency-dependence of the four and three point vertices as Γ(p, p ; q) = Γ(q), P(p; q) = P(q), and P † (p; q) = P † (q), we construct the RPA by using equations (16), (20) and (21), giving the forms where Π(q) = −T p K 0 (p; q). The 1PI part of the density correlation functions and the density vertices are also given by By using the relations such as G 22 (p) = G 11 (−p) as well as G 12 (p) = G 12 (−p), the polarization function can be reduced into The four point vertex in this approximation can be conveniently decomposed into the T-matrix T (q) given by the ladder type diagrams and the effective interaction U eff (q) including the density fluctuation, given by where γ(q) = T (q)Π(q)|f 0 , γ † (q) = f 0 |Π(q)T (q) and Here, χ R (q) is the regular part of the density-density correlation function including the vertex correction, giving the form The 1PI part of the density correlation function (88), and the density vertices (89) and (90) are also reduced into .
where A(q) = [1 + U eff (q)χ R (q)]. It can be clearly seen from the present formalism that the condensate plays an essential role in coupling the density excitation and the quasiparticle excitation as in references [60,65,66,[83][84][85], where this hybridization disappears above the critical temperature [83]. Note that because of the relation (91), five elements T 11,12,14,22,23 are needed to construct the T-matrix T , which is given by We take the following bare part of the two-particle Green's function where g(p) is the single-particle Green's function, given by where ξ p = ε p + Un 0 . At T T c , we employed the Hartree-Fock-Bogoliubov-Popov (Shohno) approximation [48,125,[133][134][135][136]. At T T c , the effective chemical potential is taken to be μ − Σ 11 (0), since the Green's function g has a pole of a gapless dispersion law at the critical temperature T c , with satisfying the Hugenholtz-Pines relation μ = Σ 11 (0). The detailed expressions of the polarization function Π 11,12,14,22 in this approximation are summarized in appendix D. Above the critical temperature, the 1PI part of the density-density correlation function is given by because g 12 (p) = Π 12,14 (p) = 0 at T T c . At the same temperature regime, the regular part is given by The T-matrix T at T T c has a diagonal form, whose matrix elements are given by For the single-particle Green's function G(p), we include many-body effects to the self-energy by using the RPA for focusing on density fluctuations, which is given by where ñ = −T p g 11 (p). The density vertex Υ in the RPA given in (97) does not satisfy the zero-frequency density vertex identity Υ(0) = 0. In the static and low-momentum limit, the density vertices given in (97) and (98) are reduced to with Π (q) = Π 11 (q) + Π 22 (q) + 2Π 14 (q) + 4Π 12 (q). (109) Each polarization function Π ij exhibits an infrared divergence. For example, in the three dimensional system at T = 0, the polarization functions exhibit the infrared divergence, giving the form Π 11,12,22,14 (p, 0) ∝ 1/|p| at small p [38]. Because of a relation g 11 (p) = −g 12 (p) in the low-energy limit, the following exact relation holds: All the infrared divergences are thus canceled out each other in Π , and then the function Π (0) converges at T < T c . By using (D1)-(D4), we have its explicit form given by where n p is the Bose-distribution function n p = 1/[exp(E p /T) − 1] with E p = ε p (ε p + 2Un 0 ). At T < T c , therefore, the density vertex parts in (107) provide Υ(0) = 0. This problem can be avoided by adopting the simplified regular part of the density-density correlation function that does not include the vertex correction, giving the form Using this simplified version, we may take a variant of the density vertices Υ s (q) and Υ s † (q), which are given by replacing A(q) in equations (97) and (98) In the low-energy limit, the simplified density vertex Υ s is reduced to Since the simplified regular part is given by χ s R (0) = χ 22 (0) + χ 14 (0), which shows the infrared divergence, the simplified density vertex Υ s satisfies the identity Υ s (0) = 0.
According to the same reason, the off-diagonal self-energy (106) does not satisfy the Nepomnyashchii-Nepomnyashchii identity Σ 12 (0) = 0. This problem is also avoided by replacing the effective interaction U eff with U s eff in the off-diagonal self-energy (106), because the infrared divergence of χ s R provides U s eff (0) = 0. As a result, the off-diagonal self-energy Σ 12 = n 0 U s eff (p) is one of the candidates to satisfy the Nepomnyashchii-Nepomnyashchii identity [36,37]. Other approaches that satisfies the Nepomnyashchii-Nepomnyashchii identity have been also discussed, including the description in terms of the hydrodynamic variables [38, 44, 48-50, 52, 137], the renormalization group approach [49,50,129,138,139], the large-N expansion [49,140] and the division approach into singular and nonsingular self-energies [38].

Density and single-particle spectral function
This section serves as the study of the density response function and the single-particle spectral function in the BEC by using the formalism developed in the previous section. The condensate density is calculated as a function of temperature, by solving the particle number equation with the non-condensate density (11), where below T c , the chemical potential satisfies the Hugenholtz-Pines relation, and the self-energies are given in (105) and (106). We performed the analytic continuation based on references [141,142].
At the low temperature regime (0.2T 0 c ), where T 0 c is the critical temperature of an ideal Bose gas, the sharp peak emerges with the satellite structure in the density response function χ ( figure 15(a)). Since the 1PI part χ 1PI is found to be negligibly small compared with the 1PR part χ 1PR , the satellite peak is mainly originated from χ 1PR part at the low temperature. This is stark contrast to the case of the multi-particle excitation in the superfluid 4 He, which provides the significant broad peak. The multi-particle excitation in  (105) and (106), and the density vertex in (97) and (98), with the effective interaction (94) including the vertex correction. We also used the 1PI part in (96) also including the vertex correction. The critical temperature is given by T c /T 0 c 1 + 1.9an 1/3 at the gas parameter an 1/3 = 10 −2 . We take the momentum q = 0.05q 0 , where T 0 the superfluid 4 He is originated from the roton-roton, maxon-maxon, and roton-maxon scattering and their bound states. Since the dispersion relation of the quasiparticle has extremum at the roton and maxon region, those provides the very large density of states owing to the van Hove singularity. This effect leads the pronounced contribution of the 1PI part to the density response function. In the present case without roton and maxon excitations, however, the satellite peak is originated from the 1PR part. For increasing temperature, the contribution from the 1PI part is enhanced ( figure 15(b)), and the density response function is mainly organized by the 1PI part close to T c ( figure 15(c)). Although the main structure of χ in figure 15(c) is the broad peak with a tail in the high frequency side, one can see the small sharp peak structure at the low frequency side, which is originated from the 1PR part. Very close to the critical temperature, the 1PR part does not show the main contribution to the density response function, because the density vertex Υ proportional to √ n 0 is small. The temperature dependences of each contribution to χ are summarized in figure 16. The density response function gives the striking sharp peak with the satellite structure, but at the intermediate temperature, the peak strength becomes weak and the satellite peak structure changes into the tail structure ( figure 16(a)). The intensity of the density response function at the critical temperature is quite small compared with the case at the low temperature. The 1PR part shows the similar behavior to the total density response function χ; however, the 1PR part vanishes at the critical temperature ( figure 16(b)). The 1PI part exhibits the striking sharp structure with a broad tail at 0.1T 0 c ; on the other hand, as the temperature increase, this sharpness vanishes with the growth of the intensity ( figure 16(c)). The spectral function of the single-particle excitation is also shown in figure 16(d). The structure of G 11 in the low temperature regime provides the sharp peak with a small satellite peak, which is the same behavior as the density response functions χ and χ 1PR . However, at high temperature such as T 0 c and T c , the satellite peak disappears, where the intensity of −ImG 11 remains the same order as those in the low-temperature case, which is in contrast to the case for the density response function. In the density response function χ, the peak of the 1PR part emergent from the single-particle excitation is suppressed by the density vertex Υ proportional to √ n 0 . In figures 15 and 16, we have discussed the structure of the density response function and the single-particle spectral function by using the self-energies (105) and (106) and the density vertex (97) and (98) both including the vertex correction. These qualitative features do not change in the case where the vertex correction is eliminated. Figure 17 shows the results with the density vertices Υ s and Υ s † satisfying the identity Υ s (0) = 0, where the self-energy contribution is still given by (105) and (106). Figure 18 shows the results with the density vertex Υ s and Υ s † as well as the self-energy contribution with the use of the effective interaction (113), which satisfies the Nepomnyashchii-Nepomnyashchii identity Σ 12 (0) = 0 [34]. Although the satellite peaks without the vertex correction in figures 17 and 18 are very slightly enhanced compared with the result in figure 15, the qualitative features remain the same.
Above the critical temperature, the density response function is exhausted by the 1PI part, where the 1PR part is absent since Υ = 0. By using the RPA, we found that the density response function at T > T c has qualitatively the same structure at T = T c , where a broad structure emerges and very long-lived collective excitations are absent. This is because the RPA describes collisionless modes, and does not describe the hydrodynamic mode. In this sense, this result indicates that there is no long-lived collisionless sound modes in a normal Bose gas. The hydrodynamic analysis in the superfluid phase can be found in reference [40]. The total density response function χ, (b) the 1PR part χ 1PR , (c) the 1PI part χ 1PI , and (d) the single-particle Green's function G 11 . We used the same Feynman diagrams as used in figure 15. Figure 17. Structure of the response functions at an 1/3 = 10 −2 and q = 0.05q 0 . (a) The total density response function χ, (b) the 1PR part χ 1PR . We used the same structure of the self-energies and the 1PI part as in figure 15. For the density vertex, we employed Υ s and Υ s † , which satisfies the zero-frequency density vertex identity Υ s (0) = 0. Figure 18. Structure of the response functions at an 1/3 = 10 −2 and q = 0.05q 0 . (a) The total density response function χ, (b) the 1PR part χ 1PR , and (c) the single-particle Green's function G 11 . We used the self-energy satisfying the Nepomnyashchii-Nepomnyashchii identity Σ 12 (0) = 0, where U eff in (105) and (106) is replaced with U s eff in (113). For the density vertex, we employed Υ s and Υ s † . We used the 1PI part in (96) including the vertex correction.
The origin of the satellite peak of the density response function can be discussed as follows: as shown in figures 15-17, and 18, the satellite peak of the density response function χ is dominantly originated from the 1PR part χ 1PR that includes the single-particle Green's function through the density vertex Υ. We thus separately treat the self-energy contribution in the single-particle Green's function to discuss the origin of the satellite peaks [34]. The self-energy contribution in the BEC involves two-parts: diagonal and off-diagonal self-energies Σ 11 (12) , which are also consists of two parts: condensate part Σ 11(12),c and non-condensate part Σ 11 (12),n . For Σ 12,n , we consider the form Σ 12,n (p) = − q U eff (p)g 12 (p + q). In order to separately analyze each contribution, we first consider the Hartree-Fock-Bogoliubov type self-energies, which can include all the contributions Σ 11 (12),c , and Σ 11(12),n , diagrammatically described in figure 19(b). In this approximation, the satellite peak can be seen ( figure 19(a)), which is consistent with the case of the Hartree-Fock-Bogoliubov-Popov type self-energies that does not include Σ 12,n .
The emergent satellite peak is possibly originated from (i) the off-diagonal self-energy Σ 12 , (ii) non-condensate part Σ 11 (12),n , or (iii) condensate part Σ 11 (12),c . Figure 19(a) shows the result of these contributions, where the self-energy contribution is selectively eliminated. The satellite peak still survives even if we eliminate the off-diagonal self-energy Σ 12 , and the non-condensate part Σ 11 (12),n . On the contrary, the satellite peak vanishes when the condensate part of the self-energy Σ 11 (12),c is absent. In the Bogoliubov approximation, where we replace the effective interaction U eff (p) with the bare interaction U, the satellite and the broadening of the sharp peak structure never emerge, which gives the Bogoliubov excitation showing the sharp peak of the quasiparticle with infinite life-time. The origin of the satellite peak is thus concluded as the many-body BEC effect, namely, the condensate part of the self-energy, which gives the interaction between the condensate and the quasiparticle in the background of the many-body density-fluctuated medium. The non-condensate part of the self-energies, showing the quasiparticle-quasiparticle interaction effect, is not important for the satellite peak, where the many-body effect of the density fluctuation is smeared out by quasiparticles with various momenta.
One of the feature of the excitation of a BEC at T = 0 is that the correspondence of the spectrum between the single-particle excitation and the collective excitation in the low-energy regime [9]. We study the temperature dependence of these two excitations and discuss this correspondence by using the effective interaction including the vertex correction ( figure 20). Except close to the critical temperature, the density spectral function is dominated by the 1PR part, and thus the peak position of χ traces that of the 1PR part ( figure 20(a)). The intensity of the 1PI part is weak and its structure is very broad compared with the 1PR part (figures 20(b) and (c)). The peak of the 1PI part is not monotonic function of the temperature, and close to T c , the peak of χ traces that of the 1PI part instead of the 1PR part, because the density vertex in the 1PR part becomes small. At the very low-temperature regime, we corroborated that the correspondence between the single-particle excitation peak and the collective excitation peak within the resolution of the numerical calculation. On the other hand, as the temperature increases, the peak of the single-particle excitation and that of the collective density excitation have a slight difference. This is due to the diminishing density vertex and the relatively increasing 1PI part as a function of the temperature.
We discuss the approximation dependence on the result of the correspondence of the spectrum between the single-particle excitation and the density collective excitation ( figure 21). In contrast to the case of The total density response function χ, (b) the 1PR part χ 1PR , (c) the 1PI part χ 1PI , and (d) the single-particle Green's function G 11 . Red, yellow, white, and black points represent the maximum peak positions of −Imχ, −Imχ 1PR , −Imχ 1PI , and −ImG 11 , respectively. We used the same self-energies, density vertex, effective interaction and the 1PI part as shown in figure 15. figure 20, we employ the density vertex satisfying the identity Υ s (0) = 0, where the vertex correction is omitted. In this case, the temperature dependence of the peak position of χ as well as χ 1PR are quite different from the case in figure 20. As a result, the temperature dependence of the peak position may change, depending on approximations, such as the absence/existence of the vertex correction. However, in the very low temperature regime, we can still find the correspondence between the peak positions between the density response function and the single-particle Green's function.
The density response function and the single-particle spectral function are shown in the ω-q plane in figure 22. The peak of the single-particle excitation traces that of the density response function at low temperature (0.1T 0 c ). This correspondence cannot be seen at T c , because of the absence of the BEC. At moderate temperature (0.5T 0 c ), although two peak positions are slightly different at high-momentum and high-energy regime, the correspondence may survives in the low-momentum and low-energy regime. At very low temperature, the peak position is well described by the Bogoliubov approximation, although the satellite peak emerges which is not reproduced by the mean-field type Bogoliubov approximation [34]. As temperature increases, the width of the single-particle spectral function becomes broad, and the phonon structure disappears at T c . The density response function at higher temperature also becomes quite broad. From these results, we can reasonably expect that it is an essential feature in BECs that the peak position of the density response function overlaps with that of the single-particle Green's function not only at the zero temperature but also in the very low but nonzero temperature regime, which is irrespective of the approximation that we take. The linear dispersion at T = 0 is analytically discussed to be originated from the identity ∂ ω Σ 11 (0) = 1 [36] for the theory satisfying Σ 12 (0) = 0. In many-body approximations at nonzero temperatures, the numerical analytic continuation makes it difficult to analyze the origin of the Momentum and frequency dependence of the total density response function χ, and the single-particle Green's function G 11 . White points in panels (a), (c), and (e) represent the maximum peak position of −ImG 11 . We used the self-energies in (105) and (106), with the effective interaction (94) including the vertex correction. For the density vertex, we employed Υ s and Υ s † . We used the 1PI part in (96) including the vertex correction. structure of the excitation spectrum. Although the linear dispersion can be originated from Σ 12 (0) in the approximation Σ 12 (0) = 0 as discussed in reference [1], this problem is important all the more in the many-body approximation at nonzero temperatures satisfying the identity Σ 12 (0) = 0 [34].
The sound speed can be estimated by inversely solving the compressibility zero-frequency sum-rule c = −n/[mχ(0)]. Since this sum-rule is exhausted by the 1PI part because of χ 1PR (0) = 0, the sound speed is exactly given by If we employ Υ s and Υ s † , we can reproduce the exact identity χ 1PR (0) = 0, because of Υ s (0) = 0. The 1PI part (96) in the static and low-momentum limits is given by and the sound speed at T T c can be estimated as where c 0 ≡ Un/m. This sound speed (117) is found to be a positive real number if we are considering the repulsive interaction U > 0, because Π (0) given in (111) is a real negative number according to the relation ∂n p /∂E p = −βn p (1 + n p ) < 0. At T T c , the sound speed is given by where the 1PI part is given in (102). This sound speed (118) is also safely a positive real number for U > 0, because of the relation One may employ the simpler regular part (112) for the 1PI part (96). In this bubble diagram case not including the vertex correction, however, we obtain an unphysical temperature-independent sound speed c = c 0 = Un/m for all temperatures below T c . Since χ s R (0) exhibits the infrared divergence at T T c , the 1PI part (96) in the static and low-momentum limits is temperature-independent, given by χ 1PI s (0) = −1/U. We discuss the temperature dependence of the sound speed c using the RPA (117) and (118) (see figure 23). In this formalism, the sound speed is temperature dependent, and the sound speeds in (117) (117) and (118) below and above the critical temperature, respectively. this temperature. The sound speed is given by c = √ 2c 0 at T = T c within the RPA including the vertex correction, where the factor 2 comes from the many-body effect in this approximation. In the Bogoliubov-Popov mean-field calculation, the sound speed is given by c = Un 0 /m, and it drops to zero at the critical temperature. At absolute zero temperature case, the sound speed (117) is approximately given by c = √ 3c 0 for p c a 1. The sound speed is overestimated in the RPA with the vertex correction, although it reproduces the same order of the sound speed in the Bogoliubov approximation at T = 0. For the consistency, further improvements may be necessary for the calculation of the sound speed derived from the RPA with the zero-frequency compressibility sum-rule.
There has been a debate whether the sound speeds below and above the critical temperature converge to the same value or show the discontinuity at the critical temperature in superfluid 4 He. The measurement of the sound velocity very close to the λ-point has the fundamental difficulty [144,148,156]. No detectable discontinuity of the sound velocity was discussed at the λ-transition [143,144,156]. The specific heat shows the jump, which suggests the second order phase transition according to the Ehrenfest relations [156], and the isothermal compressibility κ T shows not the divergence but a discontinuity at the transition point [158]. On the other hand, the logarithmic singularity of κ T is also discussed at the λ-transition [150,159]. The sound velocity near the λ-point is also theoretically investigated [96,97], and ultrasonic attenuation is also studied based on the Pippard-Buckingham-Fairbanks relations [96]. Within the present formalism, the sound speeds converge to the same value from above and below the critical temperature, where it should be noted that thoughtful treatments are needed in fluctuation regions [160].
In the formalism used in this paper, we take the Hartree-Fock-Bogoliubov-Popov approximation g(p) for constructing the building blocks and self-energies. One of the directions for the future study is to develop the self-consistent approximation, such as the self-consistent T-matrix approximation [161]. In contrast to the Fermi gas, the BEC provides the infrared divergence in the single-particle Green's function with a relation G 11 (0) = −G 12 (0), which also provides a strong constraint for the infrared divergent polarization functions, given by Π 11,22,14 (0) = −Π 12 (0) [37]. Since the exact infrared property is important for studying the low-energy properties of the BEC [38], this constraint will be important in development of the self-consistent approximation for the BEC.
The matrix formalism for BECs presented in this paper will have potential to efficiently study the exact low-energy properties of the single-particle Green's function and the density response function at nonzero temperatures as an extension of the theory at T = 0 by Gavoret and Nozières [9]. In this nonzero temperature case, we will need the forth order expansion of the self-energy with respect to ω and p i , in order to study the second sound contribution [40]. Even in this case, the Bogoliubov operator inequality and the Josephson sum-rule are still important criteria for checking the validity of the results. The present matrix formalism will also have potential to extend theories for the spinor BEC [162], the dipolar BEC [163], the collisionless sound [164], the deep inelastic scattering [165], the Bose polaron problem [166], and the renormalization-group method [50,128,138].
Ultracold atomic gases may serve as a platform for directly addressing the strong connection between the single-particle and density excitations in BECs by employing useful tools, such as the Feshbach resonance, the uniform box trap, and the spectroscopy. Theoretical concepts of BECs that should be interesting to confirm experimentally are the Josephson sum-rule, as well as the equivalence of the dispersion relations between the single-particle and collective excitations. It is also interesting to experimentally study the phonon-maxon-roton excitation in dipolar BECs not only in the collective excitation [120], but also in the single-particle excitation below and above the critical temperature by controlling the relative strength of the dipolar to the contact interactions [120]. Since the sharp maxon-roton intensity has been considered to originate from the single-particle excitation and the BEC in the superfluid 4 He [82], it will provide deeper understanding of the maxon-roton excitations as well as the connection between the single-particle and density excitations in BECs, with extending the context of superfluid 4 He.

Conclusions
We investigated the single-particle excitation and the collective density excitation in BECs by using the single-particle Green's function and the density response function. First, we revisited the earlier study presented by Gavoret and Nozières [9], with including the subsequent results given by Nepomnyashchii and Nepomnyashchii [35,36]. We extended the Nambu representation of the single-particle Green's function for BECs to correlation functions and vertex functions by making the use of the matrix formalism, which reproduces the exact properties efficiently. By following the discussion given by Gavoret and Nozières [9] with the matrix formalism, we revisited the low-energy properties of the correlation functions and the vertex functions, and the correspondence of the spectrum between the single-particle excitation and the collective excitation in the low-energy and low-momentum regime. We also present an overview of the earlier experimental and theoretical studies on the collective excitations in superfluid 4 He as well as in ultracold atomic gases. We also gave criticisms on theories casting doubt upon the conventional wisdom of the BEC: the equivalence of the dispersion relations between the single-particle excitation and the collective excitation in the low-energy and low-momentum regime. The consistency with the Bogoliubov operator inequality and the Josephson sum-rule is an important criterion for the theory contradict to the conventional wisdom.
By applying the matrix formalism, we developed a RPA for BECs to describe a single-particle Green's function and the density response function at nonzero temperatures. Depending on the presence or absence of the vertex correction, approximations provide the quantitatively different temperature dependence of the density response function and the single-particle spectral function. However, the peak positions in both functions are consistent in the very low-temperature regime, which supports the correspondence of the spectrum between the single-particle excitation and the collective excitation. Many-body effect can be seen in the satellite structure of the single-particle spectral function, which comes from the interaction between the condensate and the quasiparticles in the medium with the density fluctuation. By using the compressibility zero-frequency sum-rule, the temperature dependence of the sound speed was evaluated, where the result within the RPA including the vertex correction shows no discontinuity at the critical temperature, although careful treatments are necessary in the fluctuation region.