Self-Consistent Derivation of the Modified Gross–Pitaevskii Equation with Lee–Huang–Yang Correction

Abstract

We consider a dilute and ultracold bosonic gas of weakly-interacting atoms. Within the framework of quantum field theory, we derive a zero-temperature modified Gross–Pitaevskii equation with beyond-mean-field corrections due to quantum depletion and anomalous density. This result is obtained from the stationary equation of the Bose–Einstein order parameter coupled to the Bogoliubov–de Gennes equations of the out-of-condensate field operator. We show that, in the presence of a generic external trapping potential, the key steps to get the modified Gross–Pitaevskii equation are the semiclassical approximation for the Bogoliubov–de Gennes equations, a slowly-varying order parameter and a small quantum depletion. In the uniform case, from the modified Gross–Pitaevskii equation, we get the familiar equation of state with Lee–Huang–Yang correction.

1. Introduction

In 1924, Bose and Einstein introduced the concept of Bose–Einstein statistics and also Bose–Einstein condensation, i.e., the macroscopic occupation of the lowest single-particle state of a system of bosons [1,2]. In 1938, London suggested that the normal-superfluid phase transition of ${}^4$He is related to the Bose–Einstein condensation and to the existence of a macroscopic wavefunction for the Bose condensate [3,4]. In 1947, Bogoliubov calculated, for a uniform weakly-interacting Bose gas, the quantum depletion, i.e., the fraction of bosons that are not in the Bose condensate at zero temperature due to a repulsive interaction strength [5]. In 1957, Lee, Huang and Yang evaluated the first correcting term to the mean-field equation of state of a uniform and weakly-interacting Bose gas [6]. In 1961, Gross and Pitaevskii derived the mean-field equation for the space-dependent macroscopic wavefunction of a weakly-interacting Bose gas in the presence of an external trapping potential [7,8]. The Gross–Pitaevskii equation is the main tool used to describe the properties of the Bose–Einstein condensates, which are now routinely produced with ultracold and dilute alkali-metal atoms [9].

Some years ago, experiments with atomic gases reported evidence of beyond-mean-field effects on the equation of state of repulsive bosons [10,11]. These experimental results are quite well reproduced [12] by a modified Gross–Pitaevskii equation, which includes a beyond-mean-field correction that is the local version of the Lee–Huang–Yang term. A few years ago, Petrov suggested theoretically the existence of self-bound quantum droplets in an attractive Bose-Bose mixture, where the collapse is suppressed by a beyond-mean-field term [13]. Very recent experiments [14,15] with two internal states of ${}^{39}$K atoms in a three-dimensional configuration substantially confirm these theoretical predictions based on a modified Gross–Pitaevskii equation.

Beyond-mean-field correcting terms into the Gross–Pitaevskii equation [16,17,18], or into similar nonlinear Schrödinger equations for superfluids [19,20], are usually introduced heuristically in the spirit of the density functional theory. Here, we derive the modified Gross–Pitaevskii equation in a self-consistent way, starting from the Heisenberg equation of motion of the bosonic field operator $\widehat{\psi }(\mathbf{r},t)$ and the familiar Bogoliubov prescription of writing the quantum field operator as the sum of a classical complex field ${\psi }_0\left(\mathbf{r}\right)$, which is the order parameter or macroscopic wavefunction of the Bose–Einstein condensate, and a quantum field $\widehat{\eta }(\mathbf{r},t)$, which takes into account quantum and thermal fluctuations [21,22]. The presence of a generic external trapping potential $U\left(\mathbf{r}\right)$ is circumvented by adopting a semiclassical approximation for the Bogoliubov–de Gennes equations of the fluctuating quantum field [23,24]. In this way, at zero temperature, we obtain the local density $\tilde{n}\left(\mathbf{r}\right)$ of the out-of-condensate bosons as a function of the classical field ${\psi }_0\left(\mathbf{r}\right)$ and the corresponding equation for ${\psi }_0\left(\mathbf{r}\right)$, that is the stationary modified Gross–Pitaevskii equation with beyond-mean-field terms. From these terms, we recover the Lee–Huang–Yang correction [6] in the case of a uniform and real Bose–Einstein order parameter.

2. Quantum Field Theory of Bosons

Let us consider the bosonic quantum field operator $\widehat{\psi }(\mathbf{r},t)$ describing a non-relativistic system of confined and interacting identical atoms in the same hyperfine state. Its Heisenberg equation of motion is given by [21]:

$$i\hslash \frac{\partial }{\partial t}\widehat{\psi }(\mathbf{r},t)=\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+U\left(\mathbf{r}\right)-\mu \right]\widehat{\psi }(\mathbf{r},t)+g{\widehat{\psi }}^{+}(\mathbf{r},t)\widehat{\psi }(\mathbf{r},t)\widehat{\psi }(\mathbf{r},t).$$

(1)

where m is the mass of the atom, $U\left(\mathbf{r}\right)$ is the confining external potential, g is the strength of the interatomic potential and $\mu$ is the chemical potential, which is fixed by the conservation of the particle number N, that is an eigenvalue of the number operator:

$$\widehat{N}=\int d^3\mathbf{r}{\widehat{\psi }}^{+}(\mathbf{r},t)\widehat{\psi }(\mathbf{r},t).$$

(2)

The bosonic field operator $\widehat{\psi }(\mathbf{r},t)$ satisfies the familiar equal-time commutation rules. In Equation (1), we have assumed that the system is very dilute and such that the scattering length and the range of the interatomic interaction are much smaller than the average interatomic distance. Thus, the true interatomic potential is approximated by a local pseudo-potential:

$$V(\mathbf{r},{\mathbf{r}}^{\prime })=g{\delta }^3(\mathbf{r}-{\mathbf{r}}^{\prime }),$$

(3)

where:

$$g=\frac{4\pi {\hslash }^2{a}_s}{m}$$

(4)

is the scattering amplitude of the spin triplet channel with $a_s$ the s-wave scattering length [21].

3. Bogoliubov Prescription and Quantum Fluctuations

In a bosonic system, one can separate Bose-condensed particles from non-condensed ones by using of Bogoliubov prescription [21,22]:

$$\widehat{\psi }(\mathbf{r},t)={\psi }_0\left(\mathbf{r}\right)+\widehat{\eta }(\mathbf{r},t),$$

(5)

where:

$${\psi }_0\left(\mathbf{r}\right)=\langle \widehat{\psi }(\mathbf{r},t)\rangle$$

(6)

is the time-independent, but space-dependent complex order parameter (macroscopic wavefunction) of the Bose–Einstein condensate with $\langle \dots \rangle$ the thermal average over an equilibrium ensemble. Notice that we work at thermal equilibrium, and consequently, the thermal averages are time independent. The field $\widehat{\eta }(\mathbf{r},t)$ is the operator of quantum and thermal fluctuations, which describes out-of-condensate bosons.

The Bogoliubov prescription for the field operator $\widehat{\psi }(\mathbf{r},t)$ enables us to write the three-body thermal average in the following way:

$$\langle {\widehat{\psi }}^{+}(\mathbf{r},t)\widehat{\psi }(\mathbf{r},t)\widehat{\psi }(\mathbf{r},t)\rangle ={\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2{\psi }_0\left(\mathbf{r}\right)+2\tilde{n}\left(\mathbf{r}\right){\psi }_0\left(\mathbf{r}\right)+\tilde{m}\left(\mathbf{r}\right){\psi }_0^{*}\left(\mathbf{r}\right)+\tilde{s}\left(\mathbf{r}\right),$$

(7)

where $\tilde{n}\left(\mathbf{r}\right)=\langle {\widehat{\eta }}^{+}(\mathbf{r},t)\widehat{\eta }(\mathbf{r},t)\rangle$ is the density of non-condensed particles, while $\tilde{m}\left(\mathbf{r}\right)=\langle \widehat{\eta }(\mathbf{r},t)\widehat{\eta }(\mathbf{r},t)\rangle$ is the anomalous density and $\tilde{s}\left(\mathbf{r}\right)=\langle {\widehat{\eta }}^{+}(\mathbf{r},t)\widehat{\eta }(\mathbf{r},t)\widehat{\eta }(\mathbf{r},t)\rangle$ is the anomalous correlation [22].

Now, we obtain an equation for ${\psi }_0\left(\mathbf{r}\right)$ by taking the thermal average on Equation (1). In this way, we find:

$$\mu {\psi }_0\left(\mathbf{r}\right)=\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+U\left(\mathbf{r}\right)+g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2+2g\tilde{n}\left(\mathbf{r}\right)\right]{\psi }_0\left(\mathbf{r}\right)+g\tilde{m}\left(\mathbf{r}\right){\psi }_0^{*}\left(\mathbf{r}\right)+g\tilde{s}\left(\mathbf{r}\right),$$

(8)

which is the exact equation of motion of the Bose–Einstein order parameter ${\psi }_0\left(\mathbf{r}\right)$ [22,23]. This is not a closed equation due to the presence of the non-condensed density $\tilde{n}\left(\mathbf{r}\right)$ and of the anomalous densities $\tilde{m}\left(\mathbf{r}\right)$ and $\tilde{s}\left(\mathbf{r}\right)$. Neglecting the non-condensed density and the anomalous densities, the previous equation becomes:

$$\mu {\psi }_0\left(\mathbf{r}\right)=\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+U\left(\mathbf{r}\right)+g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2\right]{\psi }_0\left(\mathbf{r}\right),$$

(9)

which is the familiar Gross–Pitaevskii equation [7,8]. A less drastic approximation, which is called the Bogoliubov–Popov–Beliaev approximation [22,23,24], neglects only the term $\tilde{s}\left(\mathbf{r}\right)$. Then, the equation of motion of the Bose–Einstein order parameter ${\psi }_0\left(\mathbf{r}\right)$ becomes:

$$\mu {\psi }_0\left(\mathbf{r}\right)=\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+U\left(\mathbf{r}\right)+g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2+2g\tilde{n}\left(\mathbf{r}\right)\right]{\psi }_0\left(\mathbf{r}\right)+g\tilde{m}\left(\mathbf{r}\right){\psi }_0^{*}\left(\mathbf{r}\right).$$

(10)

Furthermore, this equation is not closed. We must add an equation for the non-condensed density $\tilde{n}\left(\mathbf{r}\right)$ and the anomalous density $\tilde{n}\left(\mathbf{r}\right)$ by studying the fluctuation operator $\widehat{\eta }(\mathbf{r},t)$ [22,23].

The equation of motion of the fluctuation operator $\widehat{\eta }(\mathbf{r},t)$ is obtained by subtracting Equation (10) from Equation (1). The standard Bogoliubov–Popov approximation [22,23] neglects both the non-condensate density and the anomalous terms, and it takes only linear terms of $\widehat{\eta }(\mathbf{r},t)$ and ${\widehat{\eta }}^{+}(\mathbf{r},t)$. In this way, the linearized equation of motion of the fluctuation operator reads:

$$i\hslash \frac{\partial }{\partial t}\widehat{\eta }(\mathbf{r},t)=\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+U\left(\mathbf{r}\right)-\mu +2g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2\right]\widehat{\eta }(\mathbf{r},t)+g{\psi }_0{\left(\mathbf{r}\right)}^2{\widehat{\eta }}^{+}(\mathbf{r},t).$$

(11)

4. Bogoliubov–de Gennes Equations and Their Semiclassical Approximation

The fluctuation operator can be written as:

$$\widehat{\eta }(\mathbf{r},t)=\sum _j\left[u_j\left(\mathbf{r}\right)e^{-i{E}_{j}t/\hslash }{\widehat{a}}_j+v_j\left(\mathbf{r}\right)e^{i{E}_{j}t/\hslash }{\widehat{a}}_j^{+}\right],$$

(12)

where ${\widehat{a}}_j$ and ${\widehat{a}}_j^{+}$ are bosonic operators and the real functions $u_j\left(\mathbf{r}\right)$ and $v_j\left(\mathbf{r}\right)$ are the wavefunctions of the quasi-particle and quasi-hole excitations of energy $E_j$ [23]. As a consequence, one finds:

$$\tilde{n}\left(\mathbf{r}\right)=\sum _j\left[\left(u_j{\left(\mathbf{r}\right)}^2+v_j{\left(\mathbf{r}\right)}^2\right)\langle {\widehat{a}}_j^{+}{\widehat{a}}_j\rangle +v_j{\left(\mathbf{r}\right)}^2\right],$$

(13)

and:

$$\langle {\widehat{a}}_j^{+}{\widehat{a}}_j\rangle =\frac{1}{e^{E_j/k_{B}T}-1}$$

(14)

is the Bose factor at temperature T with $k_B$ the Boltzmann constant. We stress that at zero temperature, one gets:

$$\tilde{n}\left(\mathbf{r}\right)=\sum _j{v}_j{\left(\mathbf{r}\right)}^2,$$

(15)

and also:

$$\tilde{m}\left(\mathbf{r}\right)=\sum _j{u}_j\left(\mathbf{r}\right)v_j\left(\mathbf{r}\right).$$

(16)

By inserting Equation (12) into Equation (11), we obtain the Bogoliubov–de Gennes equations:

$$\begin{array}{ccc}\hfill \widehat{L}{u}_j\left(\mathbf{r}\right)+g{\psi }_0{\left(\mathbf{r}\right)}^2{v}_j\left(\mathbf{r}\right)& =& E_j{u}_j\left(\mathbf{r}\right),\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill \widehat{L}{v}_j\left(\mathbf{r}\right)+g{\psi }_0^{*}{\left(\mathbf{r}\right)}^2{u}_j\left(\mathbf{r}\right)& =& -E_j{v}_j\left(\mathbf{r}\right).\hfill \end{array}$$

(18)

where:

$$\widehat{L}=-\frac{{\hslash }^2}{2m}{\nabla }^2+U\left(\mathbf{r}\right)-\mu +2g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2.$$

(19)

The solution of these equation can be done numerically by choosing the external potential $U\left(\mathbf{r}\right)$. However, an analytical solution can be obtained within the semiclassical approximation, where $-i\nabla \to \mathbf{k}$ and ${\sum }_j\to \int d^3\mathbf{k}/{\left(2\pi \right)}^3$ [23]. It follows that the Bogoliubov differential equations become algebraic equations:

$$\begin{array}{ccc}\hfill L_{\mathbf{k}}\left(\mathbf{r}\right)u_{\mathbf{k}}\left(\mathbf{r}\right)+g{\psi }_0{\left(\mathbf{r}\right)}^2{v}_{\mathbf{k}}\left(\mathbf{r}\right)& =& E_{\mathbf{k}}\left(\mathbf{r}\right)u_{\mathbf{k}}\left(\mathbf{r}\right),\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill L_{\mathbf{k}}\left(\mathbf{r}\right)v_{\mathbf{k}}\left(\mathbf{r}\right)+g{\psi }_0^{*}{\left(\mathbf{r}\right)}^2{u}_{\mathbf{k}}\left(\mathbf{r}\right)& =& -E_{\mathbf{k}}\left(\mathbf{r}\right)v_{\mathbf{k}}\left(\mathbf{r}\right),\hfill \end{array}$$

(21)

where:

$$L_{\mathbf{k}}\left(\mathbf{r}\right)=\frac{{\hslash }^2{k}^2}{2m}+U\left(\mathbf{r}\right)-\mu +2g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2,$$

(22)

and that the zero-temperature non-condensed density reads:

$$\tilde{n}\left(\mathbf{r}\right)=\int \frac{d^3\mathbf{k}}{{\left(2\pi \right)}^3}{v}_{\mathbf{k}}{\left(\mathbf{r}\right)}^2.$$

(23)

This quantity is also called local quantum depletion of the Bose–Einstein condensate. In addition, the local anomalous density is given by:

$$\tilde{m}\left(\mathbf{r}\right)=\int \frac{d^3\mathbf{k}}{{\left(2\pi \right)}^3}{u}_{\mathbf{k}}\left(\mathbf{r}\right)v_{\mathbf{k}}\left(\mathbf{r}\right).$$

(24)

5. Local Quantum Depletion and Generalized Gross–Pitaevskii Equation

Assuming a slowly-varying order parameter, such that the gradient term can be neglected, but also a small quantum depletion, from Equation (10), the chemical potential $\mu$ can be approximated as:

$$\mu \simeq g|{\psi }_0{\left(\mathbf{r}\right)|}^2+U\left(\mathbf{r}\right).$$

(25)

It is then straightforward to derive the elementary excitations:

$$E_{\mathbf{k}}\left(\mathbf{r}\right)=\sqrt{\frac{{\hslash }^2{k}^2}{2m}\left(\frac{{\hslash }^2{k}^2}{2m}+2g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2\right)},$$

(26)

and the real quasi-particle amplitudes:

$$\begin{array}{ccc}\hfill u_{\mathbf{k}}\left(\mathbf{r}\right)& =& \frac{1}{\sqrt{2}}{\left(\frac{\frac{{\hslash }^2{k}^2}{2m}+g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2}{E_{\mathbf{k}}\left(\mathbf{r}\right)}+1\right)}^{1/2}\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill v_{\mathbf{k}}\left(\mathbf{r}\right)& =& -\frac{1}{\sqrt{2}}{\left(\frac{\frac{{\hslash }^2{k}^2}{2m}+g{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^2}{E_{\mathbf{k}}\left(\mathbf{r}\right)}-1\right)}^{1/2}.\hfill \end{array}$$

(28)

We can now insert Equation (28) into (23), and after integration over the linear momenta, we obtain:

$$\tilde{n}\left(\mathbf{r}\right)=\frac{\sqrt{2}}{12{\pi }^2}{\left(\frac{2mg}{{\hslash }^2}\right)}^{3/2}{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^3.$$

(29)

This is the local version of the familiar Bogoliubov term for the quantum depletion, originally obtained for a uniform bosonic system, i.e., with $U\left(\mathbf{r}\right)=0$. For the local anomalous average density, after dimensional regularization [12], we find instead:

$$\tilde{m}\left(\mathbf{r}\right)=3\tilde{n}\left(\mathbf{r}\right).$$

(30)

Finally, inserting this expression into Equation (10), we get:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mu {\psi }_0\left(\mathbf{r}\right)& =\left[-\frac{{\hslash }^2}{2m}{\nabla }^2+U\left(\mathbf{r}\right)+g|{\psi }_0{\left(\mathbf{r}\right)|}^2+\frac{\sqrt{2}}{6{\pi }^2}{\left(\frac{2m}{{\hslash }^2}\right)}^{3/2}{g}^{5/2}{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^3\right]{\psi }_0\left(\mathbf{r}\right)\hfill \\ & +\frac{\sqrt{2}}{4{\pi }^2}{\left(\frac{2mg}{{\hslash }^2}\right)}^{3/2}{g}^{5/2}{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^3{\psi }_0^{*}\left(\mathbf{r}\right).\hfill \end{array}\end{array}$$

(31)

This is a modified Gross–Pitaevskii equation containing beyond-mean-field corrections due the the presence of local quantum depletion and anomalous average density. The chemical potential $\mu$ of Equation (31) is fixed by the normalization condition:

$$N=\int d^3\mathbf{r}\left[|{\psi }_0{\left(\mathbf{r}\right)|}^2+\tilde{n}\left(\mathbf{r}\right)\right]=\int d^3\mathbf{r}\left[|{\psi }_0{\left(\mathbf{r}\right)|}^2+\frac{\sqrt{2}}{12{\pi }^2}{\left(\frac{2mg}{{\hslash }^2}\right)}^{3/2}{\left|{\psi }_0\left(\mathbf{r}\right)\right|}^3\right]$$

(32)

with N the total number of bosons.

It is important to stress that, in the case of a uniform and real Bose–Einstein condensate, Equation (31) gives:

$$\mu =g{n}_0+2g\tilde{n}+3g\tilde{n}=gn+4g\tilde{n}=gn+g{n}_0\frac{32}{3\sqrt{\pi }}\sqrt{n_0{a}_s^3},$$

(33)

which is the chemical potential with the familiar beyond-mean-field Lee–Huang–Yang correction [6] under the assumption of small quantum depletion. Clearly, one does not obtain this zero-temperature result neglecting the anomalous average density $\tilde{m}$.

6. Conclusions

We have shown that a modified Gross–Pitaevskii equation with local beyond-mean-field terms can be obtained in a straightforward way from a quantum-field-theory formulation without invoking the density functional theory. However, the derivation is not exact because one performs some approximations on the spectrum of elementary excitations and on the spatial dependence of the macroscopic wavefunction of the Bose–Einstein condensate. In [25], a different derivation of a zero-temperature stationary modified Gross–Pitaevskii equation without external confinement is shown, but which also takes into account anomalous averages. A more formal and mathematical approach to the Hartree–Fock–Bogoliubov methods to obtain time-dependent modified Gross–Pitaevskii equations can be found in [26]. In the near future, we want to derive and use coupled modified Gross–Pitaevskii equations for studying Bose-Bose mixtures under double-well confinement and spin-orbit coupling, extending our previous results [27,28].