Dynamics of Bose-Einstein condensate with account of pair correlations

The system of dynamic equations for Bose-Einstein condensate at zero temperature with account of pair correlations is obtained. The spectrum of small oscillations of the condensate in a spatially homogeneous state is explored. It is shown that this spectrum has two branches: the sound wave branch and branch with an energy gap.


I. INTRODUCTION
Bose-Einstein condensate of a low-density system of weakly interacting Bose particles at zero temperature is usually described by the Gross-Pitaevskii equation [1,2], which is nowadays widely used for study of the condensates created in magnetic and laser traps [3,4]. The Gross-Pitaevskii equation is obtained in the selfconsistent field approximation, where the short-range correlations between particles are neglected. In this case, the Bose system is described in terms of the coherent state vector [5]. Meanwhile, the consideration of the pair correlations, being essential at short distances, turns out to be important even for systems with low densities since it leads to some qualitatively new results. For example, in a dilute classical gas the account for pair correlations allows us to obtain the integral of collisions in the kinetic equation and, therefore, all effects which are described by the Boltzmann equation [6].
In the present work, we obtain the system of dynamical equations in the approximation in which besides the oneparticle anomalous averages, the two-particle correlations are also taken into account but the correlations of more number of particles are neglected. Small oscillations on the background of a spatially homogeneous equilibrium state are studied. It is shown that when considering pair correlations in Bose-Einstein condensate, there exist two branches of elementary excitations. One of them is acoustic, another one has an energy gap in the long wavelength limit.

II. EQUATIONS FOR AVERAGES OF THE FIELD OPERATORS
An arbitrary operator in the Heisenberg representation where in the second-quantized representation the Hamiltonian can be written as a sum of operators of the kinetic energy and the energy of pair interaction H = H 1 + H 2 , (2) Here m is the bose-particle mass, U 0 (r) -the energy of particle in the external field, U (|r 1 −r 2 |) -the particle interaction potential. The field operators satisfy the usual commutation relations. Let Ψ is the average of field operator. Then we can write down our field operator separating the c-number and operator parts: The operator part is defined so that it meets the obvious relations: Here the averaging is implied in the sense of the quasiaverages for systems with broken phase symmetry [7,8].
We will assume that the normal averages, which are invariant under the phase transformation of the field operators Ψ → Ψ ′ = e iα Ψ, and the anomalous averages, where this invariance is broken, are both nonzero. It is to be noted that the property of superfluidity is connected with just the existence of these anomalous averages. We introduce the the following notation for the anomalous average of the field operator η(r, t) = Ψ(r, t) , η * = Ψ † (r, t) .
The averages of the products of several field operators can be written in terms of the averages of the products of operators ξ, ξ † , which will be called as the overcondensate operators . For example, for the case of two field operators, using Eq.(5), we have Similarly, we can write down the averages of a greater number of field operators. They will also contain the averages of a greater number of the overcondensate operators of the form and etc. Setting successively the operator A in the Heisenberg equation (1) equal to Ψ † Ψ, ΨΨ, Ψ † Ψ † , . . . and carrying out averaging, we will get the infinite chain of coupled equations for the averages Ψ , Ψ † Ψ , ΨΨ , Ψ † Ψ † . . . , which is similar to the Bogoliubov-Born-Green-Kirkwood-Yvon chain [6] in the kinetic theory of classical gases. So, the equation for the average of the field operator (6) has the form: and the equations for the normal and anomalous pairwise correlations are written as follows In the following we will describe the condensate by using the one-particle averages (6) and restrict ourselves to considering only pair correlations of the overcondensate operators introduced by relations (4), having defined the following correlation functions: The averages of a greater number of the overcondensate operators will be neglected that seems to be acceptable for sufficiently diluted systems. Functions (11) have the obvious symmetry properties When only the pair correlations are taken into account, from (8) -(10) one can get the closed system of equations One should note that this system of equations is invariant under time-reversal transformation, because along with solutions η(r, t), g(r, r ′ , t), τ (r, r ′ , t) it also has the solutions η * (r, −t), g * (r, r ′ , −t), τ * (r, r ′ , −t). Neglecting the pair correlations τ (r, r ′ , t) and g(r, r ′ , t), equation (13) takes the form of the Gross-Pitaevskii equation [1,2]. In what follows, where it will not cause confusion, as in equations (13) -(15), for brevity the explicit timedependance of the averages will be omitted.
The average of the operator of the total number of particles N is given by the formula and the particle number density is, obviously, n(r, t) = η * (r, t)η(r, t) + g(r, r, t).

III. LOCAL FORM OF EQUATIONS
Equations (13) -(15) are integro-differential. While studying the states that slowly vary on the scales comparable to characteristic radius r 0 of action of the interparticle interaction potential U (|r − r ′ |), we can pass to differential equations. The pair correlation functions (11) depend on two coordinates r, r ′ . It is convenient to These functions slowly change depending on the pair's center-mass coordinate R on distances of the order of the interparticle potential radius r 0 . The correlation functions can be presented in the form: In the following, only the term with k = 0 will be taken into account in these sums. It means that instead of the exact functions τ (R, ρ) g(R, ρ), we will use the functions which are averaged over a macroscopic volume V 0 ∼ L 3 , where L ≫ r 0 : This approximation is acceptable if one considers perturbations on spatial scales that significantly exceed the radius of action of the interparticle potential.
It is worth noting that in obtained equations the behavior of the interparticle interaction potential at short distances plays an important role. The form of the potential is poorly known here. Moreover, for many model potentials such as, for example the Lennard-Jones potential, it is assumed that at short distances it goes to infinity. Note also that the use of model potentials which go to infinity at short distances in some cases leads to considerable difficulties, because such potentials do not have the Fourier representation. Meanwhile, the requirement of "impermeability" of atoms at arbitrary high pressures is very strict, since there should exist a pressure at which an atom will be "crushed" and stop existing as a separate structural unit. Therefore, in our opinion, it is physically reasonable and natural to use the potentials, which take a finite value at short distances. It should also be noted that the quantum-chemical calculations give potentials with a finite, albeit large, value at zero [9,10].
In this local approximation, the system of equations (13) -(15) takes the form: Here we used the designations: τ (R) ≡ τ (r), g 0 (R) ≡ g(r), U 0 ≡ U (r)dr. Since the magnitude of the potential energy of interaction of atoms at short distances U (0) is poorly known, it will be regarded as a phenomenological adjustable parameter. For specific calculations we will use a simple model potential of "semi-transparent sphere" form: The parameter u is assumed to be positive. In this case U (0) = u, U 0 = uv, where v ≡ 4π 3 r 3 0 is the "atomic volume". The potential (22) has earlier been used in the studies of the Bose systems (see e.g. [11]). Further, we will analyze the system of equations (19) -(21) and use the potential (22) for some estimations.

IV. EQUILIBRIUM SPATIALLY HOMOGENEOUS STATE
In this section we will consider the equilibrium state of a spatially homogeneous system. In this case, the functions η(r) ≡ η, τ (r) ≡ τ , g(r) ≡ g do not depend on coordinates and time. For sufficiently weak interaction, the most of particles will be in the single-particle condensate at zero temperature [12,13]. Therefore, we will assume the normal correlation function g = 0 at the equilibrium state, so that the equilibrium particle number density n = |η| 2 . Then, in the absence of the external field, from (19) - (21) there follow the equations which determine the equilibrium state of the system: Let's write down the complex quantities extracting their moduli and phases: η = η 0 e iα , τ = τ 0 e iβ . From (25) it follows that sin (2α − β) = 0. Thus, there are two possibilities 2α− β = 0 or 2α− β = π. We have to choose the second possibility, since only in this case equations (23),(24) have the physically correct solutions: With this choice of phase τ = −τ 0 e 2iα . After eliminating the chemical potential µ from these equations, we get Note that there does not exist the solution for the system of equations (26) and (27) such that only the singleparticle condensate exists η 0 = 0 and the pair condensate is absent τ 0 = 0. This feature of the Bose systems has been earlier pointed out [14]. From (28) it follows the relation between density and the pair correlation: Here we have used the relation U 0 = uv (v is the "atomic volume"), which was obtained above for the potential (22). Hence it follows the restriction vτ 0 < 1/2. The chemical potential proves to be positive. The dependence of the pair correlation on the particle density is given by the relation Here nv = v/Ω, where Ω is the volume per one particle. In a dilute system Ω ≫ v and nv ≪ 1. In this case n ≈ τ 0 .

V. SPECTRUM OF ELEMENTARY EXCITATIONS
In this section we will consider the propagation of small perturbations in a spatially homogeneous system. Setting η(r, t) = η 0 + δη(r, t), τ (r, t) = τ 0 + δτ, g(r, t) = δg(r, t), (32) we get from (19) -(21) the system of linearized equations: It is convenient to pass from the complex variables δη(r, t), δτ (r, t), δg(r, t) to the real ones δΨ(r, t) = δη(r, t) + δη * (r, t), δΦ(r, t) = i δη(r, t) − δη * (r, t) , δΘ(r, t) = δτ (r, t) + δτ * (r, t), δΛ(r, t) = i δτ (r, t) − δτ * (r, t) . (36) For real quantities the system of linearized equations takes the form: Assuming that the dependence of the fluctuations on the coordinates and time is of the form exp [ikr − ωt], one gets the system of homogeneous linear algebraic equations. From the condition of equality to zero of its determinant we obtain the biquadratic equation that determines the dispersion laws of possible excitations: Here where ε k = 2 k 2 2m is the free particle energy. The coefficients in Eq.(43) have the form (44) The system admits a solution in the form of spatially homogeneous oscillations with a frequency ω 0 , which is determined by the following formula: The dependence of this frequency on the density, determined with account of the relation (31), is shown in Fig.1 (curve 1). The biquadratic equation (42) has two solutions, which determine two excitation branches: The solution ω − at small wave numbers gives the sound branch ω 2 − = c 2 k 2 , where the square of speed of sound is determined by the formula: The dependence of the speed of sound on density is shown in Fig. 1 (curve 2). The solution ω + corresponds to the excitation branch with an energy gap. For small k the dependence of the frequency on the wave number has the form where α = ω −2 0 U 3 (0) + 10 U 2 (0)U 0 (η 2 0 + τ 0 ) + + 4 U(0)U 2 0 (8τ 2 0 + 3η 4 0 + 10η 2 0 τ 0 ) + 4U 3 0 τ 2 0 (η 2 0 + 8τ 0 ) . (49) Although, strictly speaking, these equations are applicable for the long-wave excitations, but the solutions (46) give the reasonable values also for large k: One branch turns into to the dispersion law of a free particle and the other branch gives the dispersion law of a pair of coupled particles. The branches of elementary excitations in the Bose system with account of pair correlations are shown in Fig.2. This figure also shows the Bogolyubov dispersion law ω B = ε k (ε k + 2uvn).

VI. CONCLUSION
We have obtained the system of differential equations (19) -(21) which describes the dynamics of the Bose-Einstein condensate with account of pair correlations. The spectrum of small oscillations in a spatially homogeneous system was studied. It was shown that there are two branches of elementary excitations: one branch with the sound dispersion law in the long wavelength limit and the second branch which has in this limit an energy gap.
It should be noted that the question of the possible existence in the Bose systems of excitations with an energy gap has a long history and has been discussed in many papers (see for e.g. [15]- [20]). The possible existence, in addition to the phonon branch of the spectrum, of the other branch with an energy gap was discussed at a qualitative level in [13, p. 322].
In the experimental paper [21] the absorption of microwave radiation was discovered in the superfluid helium at a frequency of about 180 GHz, which authors attributed to the creation of single rotons. However, due to the fact that the momentum of a roton is by orders of magnitude greater than the momentum of a photon at the given frequency, it was suggested in [22] that the observed absorption owes its existence to the presence of excitations with energy gap in the superfluid helium.
The branches of elementary excitations (one of which is sound and another has an energy gap), which are obtained theoretically in this work, can be considered as a confirmation of the qualitative arguments in favor of the possible existence of the gap excitations in the superfluid helium and the need for modification of the energy spectrum of He-II, which were formulated in [22].