⇦Chinese Physics Letters, 2022, Vol. 39, No. 6, Article code 060301 Phonon Stability of Quantum Droplets in Dipolar Bose Gases Fan Zhang (张帆) and Lan Yin (尹澜)* Affiliations School of Physics, Peking University, Beijing 100871, China Received 7 April 2022; accepted 15 May 2022; published online 29 May 2022 *Corresponding author. Email: yinlan@pku.edu.cn Citation Text: Zhang F and Yin L 2022 Chin. Phys. Lett. 39 060301    Abstract Stabilized by quantum fluctuations, dipolar Bose–Einstein condensates can form self-bound liquid-like droplets. However in the Bogoliubov theory, there are imaginary phonon energies in the long-wavelength limit, implying dynamical instability of this system. A similar instability appears in the Bogoliubov theory of a binary quantum droplet, and is removed due to higher-order quantum fluctuations as shown recently [Gu Q and Yin L 2020 Phys. Rev. B 102 220503(R)]. We study the excitation energy of a dipolar quantum droplet in the Beliaev formalism, and find that quantum fluctuations significantly enhance the phonon stability. We adopt a self-consistent approach without the problem of complex excitation energy in the Bogoliubov theory, and obtain a stable anisotropic sound velocity which is consistent with the superfluid hydrodynamic theory, but slightly different from the result of the extended Gross–Pitaevskii equation due to quantum depletion. A modified Gross–Pitaevskii equation in agreement with the Beliaev theory is proposed, which takes the effect of quantum fluctuations into account more completely.
cpl-39-6-060301-fig1.png
 DOI:10.1088/0256-307X/39/6/060301 © 2022 Chinese Physics Society Article Text Quantum droplets have been realized in different systems, including single-component dipolar Bose gases of $^{164}\rm{Dy}$,[1–5] $^{166}\rm{Er}$ atoms,[6] nonmagnetic binary $^{39}\rm{K}$ mixture,[7–9] and heteronuclear $^{39}\rm{K}$–$^{87}\rm{Rb}$ mixtures.[10] There are also proposals to create quantum droplets in binary dipolar mixtures.[11,12] Usually an ultracold bosonic gas is a dilute weakly interacting system in which the contribution from quantum fluctuations is typically small compared with the mean-field contribution, but quantum droplets are just exceptions. Due to the repulsive Lee–Huang–Yang (LHY) energy[13] which overcomes the net attractive mean-field energy near the mean-field unstable point, a nonmagnetic binary bosonic mixture displays a self-bound liquid-like state,[14] as observed in experiments.[7–10] The LHY energy[15] also stabilizes the single-component dipolar quantum droplet in the dipole-dominated regime where the dipolar interaction dominates over the repulsive s-wave interaction. Phase-coherent droplet arrays have also been found both theoretically and experimentally,[16–19] where broken translational symmetry and superfluid order appear at the same time, consistent with a supersolid.[20,21] In a nonmagnetic binary boson mixture, the droplet density can be tuned by spin-orbit coupling.[22] In a boson-fermion mixture, fermion pairing is helpful to the formation of a quantum droplet.[23] In a nonmagnetic boson mixture, quantum droplets are formed in the mean-field unstable region where the total mean-field energy is attractive. As a result, in the Bogoliubov theory, there are imaginary excitation energies in the long wavelength limit, implying the dynamical instability. However, these excitations have little contribution to the LHY energy,[14] and it is postulated that these excitations are stable after integrating out high-energy excitations. The LHY energy is thus obtained in the Bogoliubov theory with the imaginary part neglected, and used in the construction of the extended Gross–Pitaevskii equation (EGPE). The EGPE approach has been quite successful in interpreting experimental results,[24–27] with support to some extent by quantum Monte Carlo calculations.[28–31] Despite the success of the EGPE, how the dynamical instability is avoided in quantum droplets, or even, whether the Bose–Einstein-condensation (BEC) state is a stable ground state has become a pressing issue. The pairing[32] and Gaussian states[33] have been proposed as ground states. Recently in a Beliaev approach, it is found that the phonon energy is always positive after considering the interaction between spin and density fluctuations,[34] showing that the BEC state is a stable ground state for a nonmagnetic binary quantum droplet. In the Bogoliubov theory of a uniform dipolar Bose gas,[15] in the quantum droplet region the excitation energy in the direction perpendicular to the polarization is imaginary in the long-wavelength limit, suggesting dynamic instability. The Gaussian state was also proposed to be the ground state of dipolar quantum droplets.[33] In this work, we study the excitation stability of dipolar quantum droplets by going beyond the Bogoliubov theory and calculating the excitation energy in the Beliaev formalism.[35] In the quantum droplet region, as there are imaginary excitation energies in the Bogoliubov theory, we propose a self-consistent method to take account of the fluctuation effect from the start and avoid the complex excitation energy. In the long-wavelength limit the excitation energy is stable with relative dipole-dipole interaction (DDI) strength $\epsilon_{\rm{dd}} < 1.35$ for the $^{166}\rm{Er}$ system. We obtain the anisotropic phonon velocity consistent with the superfluid hydrodynamic theory, but slightly different from the EGPE result due to the quantum depletion. We propose a modified EGPE to produce the correct phonon velocity. Finally we discuss implication of our results in experiment. Beliaev Theory of a Dipolar Bose Gas. We consider a homogeneous gas of bosonic particles with mass $m$ and a finite magnetic dipole moment $d$ at zero temperature, with condensate density $n_0$. The magnetic dipoles align along the $z$ axis with the interaction potential given by $$\begin{align} V_{\rm{int}}({\boldsymbol{r}})=g\Big[\delta({\boldsymbol{r}}) +\frac{3\epsilon_{\rm{dd}}}{4\pi|{\boldsymbol{r}}|^3} \Big(1-3\frac{z^2}{|{\boldsymbol{r}}|^2}\Big)\Big],~~ \tag {1} \end{align} $$ where $g$ is the s-wave coupling constant, $g=4\pi\hbar^2a_{\rm{s}}/m$ with $a_{\rm{s}}$ being the s-wave scattering length, $\epsilon_{\rm{dd}}=a_{\rm{dd}}/a_{\rm{s}}$ is the relative interaction strength with $a_{\rm{dd}}=\mu_0d^2m/12\pi\hbar^2$ being the characteristic length of DDI and $\mu_0$ the magnetic permeability in vacuum. The Hamiltonian of a dipolar Bose gas is given by $$\begin{alignat}{1} \hat{H}={}&\int d^3r\hat{\psi}^†({\boldsymbol{r}})\Big(-\frac{\hbar^2\nabla^2}{2m}\Big)\hat{\psi}({\boldsymbol{r}}) \\ &+\frac{1}{2}\int d^{3}r d^{3}r'\hat{\psi}^†({\boldsymbol{r}})\hat{\psi}^†({\boldsymbol{r}'})\\ &\cdot V_{\rm{int}}({\boldsymbol{r}-\boldsymbol{r}'})\hat{\psi}({\boldsymbol{r}'})\hat{\psi}({\boldsymbol{r}}),~~ \tag {2} \end{alignat} $$ where $\hat{\psi}(\boldsymbol{r})$ is the boson field operator. The Fourier transform of the interaction potential is given by $$\begin{align} U({\boldsymbol{p}})=g[1+\epsilon_{\rm{dd}}(3\rm{{\cos}}^2\theta-1)],~~ \tag {3} \end{align} $$ where $\theta$ is the angle between $\boldsymbol{p}$ and the $z$ axis.[36] In the BEC phase, there are a macroscopic number of atoms in the zero-momentum state, and the boson field operators of the zero-momentum state $\hat{a}_0$ and $\hat{a}_0^†$ can be replaced by a $c$-number $N_0^{1/2}$, where $N_0=n_0V$ is the number of bosons in the condensate and $V$ is the volume.[37] The Boson field operator can be separated into two parts, $$ \begin{alignat}{1} &\hat{\psi}({\boldsymbol{r}})\equiv \hat{a}_0/V^{1/2}+\widetilde{\psi}({\boldsymbol{r}})\approx n_0^{1/2}+\widetilde{\psi}({\boldsymbol{r}}),~~ \tag {4}\\ & \hat{\psi}^†({\boldsymbol{r}})\equiv \hat{a}_0^†/V^{1/2}+\widetilde{\psi}^†({\boldsymbol{r}})\approx n_0^{1/2}+\widetilde{\psi}^†({\boldsymbol{r}}),~~ \tag {5} \end{alignat} $$ where $\widetilde{\psi}^†(\boldsymbol{r})$ is the field operator of bosons outside the condensate. In this work, we consider the dilute case in which the interaction term in the Hamiltonian can be expanded in terms of $n_0^{1/2}$. The Bogoliubov Hamiltonian describes the quadratic fluctuation, given by $$\begin{align} \hat{H}_{\rm B}={}&\frac{1}{2}n_0^2 VU(0)+\frac{1}{2}n_0\sum\limits_{\boldsymbol{k}}U({\boldsymbol{k}}) (\hat{a}_{\boldsymbol{k}}^†\hat{a}_{{-\boldsymbol{k}}}^†+\hat{a}_{\boldsymbol{k}}\hat{a}_{{-\boldsymbol{k}}})\\ & +\sum\limits_{\boldsymbol{k}}[\varepsilon_{\boldsymbol{k}}+n_0 U(0)+n_0U({\boldsymbol{k}})]\hat{a}_{\boldsymbol{k}}^†\hat{a}_{\boldsymbol{k}}.~~ \tag {6} \end{align} $$ Here, $\hat{a}_{\boldsymbol{k}}$ is the boson annihilation operator in the momentum space. The quasi-particle excitation energy is obtained by diagonalizing the Bogoliubov Hamiltonian, $$ \varepsilon_{\boldsymbol{p}}=\sqrt{\varepsilon^{(0)}_{\boldsymbol{p}}(2n_0U(\boldsymbol{p}) +\varepsilon^{(0)}_{\boldsymbol{p}})},~~ \tag {7} $$ where $\varepsilon^{(0)}_{\boldsymbol{p}}=\hbar^2 p^2/(2m)$ is the kinetic energy of the atom. In the BEC phase, the boson Green's function is a matrix, defined as[35] $$ {\rm G}({\boldsymbol{p}},t_1-t_2)=-i\langle T\{ \varPsi_{\boldsymbol{p}}(t_1)\varPsi_{\boldsymbol{p}}^†(t_2)\}\rangle,~~ \tag {8} $$ where $\varPsi^†_{\boldsymbol{p}}(t)=[\hat{a}^†_{\boldsymbol{p}}(t),\hat{a}_{-\boldsymbol{p}}(t)]$, and $T\{\dots\}$ is the time-ordering operator. In the momentum and energy space, Dyson's equation is given by $$ {\rm G}(p)={\rm G}^{(0)}(p)+{\rm G}^{(0)}(p){\rm \varSigma}(p){\rm G}(p),~~ \tag {9} $$ where $p=(\boldsymbol{p},p^0)$, and ${\rm \varSigma}(p)$ is the self-energy. In the matrices of Green's function and self-energy, $$\begin{align} &{\rm G}(p)=\begin{bmatrix} G_{11}(p) & G_{12}(p) \\ G_{21}(p) & G_{11}(-p) \end{bmatrix},~~ \tag {10} \end{align} $$ $$\begin{align} &{\rm G}^{(0)}(p)=\begin{bmatrix} G_0(p) & 0 \\ 0 & G_0(-p) \end{bmatrix},~~ \tag {11} \end{align} $$ $$\begin{align} &{\rm \varSigma}(p)=\begin{bmatrix} \varSigma_{11}(p) & \varSigma_{12}(p) \\ \varSigma_{21}(p) & \varSigma_{11}(-p) \end{bmatrix},~~ \tag {12} \end{align} $$ the subscript $11$ refers to an ingoing and an outgoing external line in the Feynman diagrams, the subscript $12$ refers to two outgoing lines, and the superscript $21$ of $\varSigma_{21}$ refers to two ingoing lines. The free-boson Green's function is given by $$ G_0(p)=1/(p^0-\varepsilon^{(0)}_{\boldsymbol{p}}-\mu+i0^+),~~ \tag {13} $$ where $\mu$ is the chemical potential. From the pole-equation of Green's function $\det|{\rm G}^{-1}(p)|=0$, the excitation spectrum can be solved. For a dilute Bose gas, the self-energy matrix ${\rm \varSigma}({\rm p})$ can be expanded in terms of the condensate density $n_0$, $$ {\rm \varSigma}({\rm p})={\rm \varSigma}^{(1)}({\rm p})+{\rm \varSigma}^{(2)}({\rm p})+\cdots,~~ \tag {14} $$ where the first-order self-energy matrix is given by $$ {\rm \varSigma}^{(1)}({\rm p})=\begin{bmatrix} n_0U(0)+n_0U(\boldsymbol{p}) & n_0U(\boldsymbol{p}) \\ n_0U(\boldsymbol{p}) & n_0U(0)+n_0U(\boldsymbol{p}) \end{bmatrix}.~~ \tag {15} $$ The chemical potential can be expanded as well, and the first-order term is given by $\mu^{(1)}=n_0 U(0)$. The first-order Green's function is given by $$\begin{align} G_{11}(p)={}&\frac{A_{\boldsymbol{p}}}{p^0-\varepsilon_{\boldsymbol{p}}+i0^+}-\frac{B_{\boldsymbol{p}}}{p^0+\varepsilon_{\boldsymbol{p}}-i0^+},~~ \tag {16} \end{align} $$ $$\begin{align} {G}_{12}(p)={}&-C_{\boldsymbol{p}}\Big[\frac{1}{p^0-\varepsilon_{\boldsymbol{p}}+i0^+}-\frac{1}{p^0+\varepsilon_{\boldsymbol{p}}-i0^+}\Big] \\ ={}&{G}_{21}(p),~~ \tag {17} \end{align} $$ where $$\begin{align} A_{\boldsymbol{p}}={}&[\varepsilon_{\boldsymbol{p}}+\varepsilon_{\boldsymbol{p}}^{(0)} +n_0U({\boldsymbol{p}})]/2\varepsilon_{\boldsymbol{p}},\\ B_{\boldsymbol{p}}={}&[-\varepsilon_{\boldsymbol{p}}+\varepsilon_{\boldsymbol{p}}^{(0)} +n_0U({\boldsymbol{p}})]/2\varepsilon_{\boldsymbol{p}},\\ C_{\boldsymbol{p}}={}&n_0U({\boldsymbol{p}})/2\varepsilon_{\boldsymbol{p}}.~~ \tag {18} \end{align} $$ From the pole of Green's function, the quasiparticle excitation energy can be obtained, as same as that in Eq. (7) from the Bogoliubov theory. In the droplet regime with $\epsilon_{\rm{dd}}>1$, the interaction $U(\boldsymbol{p})$ is attractive in the direction perpendicular to the polarization with $\theta=\pi/2$, and the excitation energy becomes imaginary in the long wavelength limit, indicating the dynamical instability.
cpl-39-6-060301-fig1.png
Fig. 1. Feynman diagrams of second-order self-energies. A filled rectangle denotes a sum of two rectangles, one for a direct interaction and the other for an exchange interaction. In the diagrams with four external lines, two lines are with zero momentum. The internal lines are first-order Green's functions. The line with one arrow is the normal part of Green's function, and that with two arrows is the anomalous part. A dashed line represents subtracting the contribution of $G^{(0)}$ to avoid double-counting.
In Ref. [35], Beliaev calculated the second-order self-energy from one-loop diagrams and obtained the correction to excitation spectrum of a Bose gas with an s-wave interaction. Here for a dipolar Bose gas, we adopt Belieav's approach to consider the fluctuation effect beyond the Bogoliubov theory, and obtain the second-order self-energy shown in Fig. 1, given by $$\begin{align} \varSigma_{11}^{(2)}&(p)=\frac{1}{2}n_0\int \frac{d \boldsymbol{q}}{(2\pi)^3} \frac{1}{p^0-\varepsilon_{\boldsymbol{k}}-\varepsilon_{\boldsymbol{q}}+i0^+}\\ &\cdot\{[U({\boldsymbol{p}})+U({\boldsymbol{k}})]^2(A_{{\boldsymbol{k}}};B_{\boldsymbol{q}}) +2[U({\boldsymbol{p}})+U({\boldsymbol{k}})]\\ &\cdot[U({\boldsymbol{p}})+U({\boldsymbol{q}})]C_{\boldsymbol{k}}C_{\boldsymbol{q}} +2[U^2({\boldsymbol{k}})+U({\boldsymbol{k}})U({\boldsymbol{q}})]\\ &\cdot A_{\boldsymbol{k}}A_{\boldsymbol{q}} -2[U({\boldsymbol{p}})+U({\boldsymbol{q}})][U({\boldsymbol{k}})+U(\boldsymbol{q})](A_{\boldsymbol{k}};C_{\boldsymbol{q}})\}\\ &-\frac{1}{2}n_0\int \frac{d \boldsymbol{q}}{(2\pi)^3} \frac{1}{p^0+\varepsilon_{\boldsymbol{k}}+\varepsilon_{\boldsymbol{q}}-i0^+}\\ &\cdot\{[U({\boldsymbol{p}})+U({\boldsymbol{k}})]^2(A_{\boldsymbol{k}};B_{\boldsymbol{q}})+2[U({\boldsymbol{p}})+U({\boldsymbol{k}})]\\ &\cdot[U({\boldsymbol{p}})+U({\boldsymbol{q}})]C_{\boldsymbol{k}}C_{\boldsymbol{q}}+2[U^2(\boldsymbol{k})+U({\boldsymbol{k}})U({\boldsymbol{q}})] B_{\boldsymbol{k}}B_{\boldsymbol{q}}\\ & -2[U({\boldsymbol{p}})+U({\boldsymbol{q}})][U(\boldsymbol{k})+U({\boldsymbol{q}})](B_{\boldsymbol{k}};C_{\boldsymbol{q}})\}\\ & +\frac{1}{2}n_0\int \frac{d \boldsymbol{q}}{(2\pi)^3} U^2(\boldsymbol{q})\Big(\frac{1}{\varepsilon_{\boldsymbol{k}}}+\frac{1}{\varepsilon_{\boldsymbol{q}}}\Big)\\ &+\int \frac{d \boldsymbol{q}}{(2\pi)^3}[U(0)+U({\boldsymbol{k}})]B_{\boldsymbol{q}}\\ &-\frac{1}{2}n_0\int \frac{d \boldsymbol{q}}{(2\pi)^3} U^2(\boldsymbol{q}) \Big[\frac{4}{\varepsilon^{(0)}_{\boldsymbol{p}}-\varepsilon^{(0)}_{\boldsymbol{q}}-\varepsilon^{(0)}_{\boldsymbol{k}}+i0^+}\\ &+\Big(\frac{1}{\varepsilon_{\boldsymbol{k}}}+\frac{1}{\varepsilon_{\boldsymbol{q}}}\Big)\Big], \end{align} $$ $$\begin{align} \varSigma_{12}^{(2)}&(p)=\frac{1}{2}n_0\int \frac{d \boldsymbol{q}}{(2\pi)^3}\Big(\frac{1}{p^0-\varepsilon_{\boldsymbol{k}}-\varepsilon_{\boldsymbol{q}}+i0^+}\\ &-\frac{1}{p^0+\varepsilon_{\boldsymbol{k}}+\varepsilon_{\boldsymbol{q}}-i0^+}\Big) \cdot\{[U({\boldsymbol{p}})+U(\boldsymbol{k})]\\ &\cdot[U({\boldsymbol{p}})+U({\boldsymbol{q}})](A_{\boldsymbol{k}};B_{\boldsymbol{q}}) +2([U({\boldsymbol{p}})+U({\boldsymbol{q}})]^2\\ &+U({\boldsymbol{q}})U({\boldsymbol{k}})+U({\boldsymbol{q}})^2)C_{\boldsymbol{k}}C_{\boldsymbol{q}}-[U({\boldsymbol{k}})+U({\boldsymbol{q}})]\\ &\cdot[U({\boldsymbol{k}})+U({\boldsymbol{p}})][(A_{\boldsymbol{k}}+B_{\boldsymbol{k}});C_{\boldsymbol{q}}]\}\\ & +\frac{1}{2}n_0\int \frac{d \boldsymbol{q}}{(2\pi)^3} U^2(\boldsymbol{q}) \Big(\frac{1}{\varepsilon^{(0)}_{\boldsymbol{q}}}-\frac{1}{\varepsilon_{\boldsymbol{q}}}\Big),~~ \tag {19} \end{align} $$ where ${\boldsymbol{k}}={\boldsymbol{p}} - {\boldsymbol{q}}$, the symbol (;) denotes a symmetrized product, $(A_{\boldsymbol{k}};B_{\boldsymbol{q}})=A_{\boldsymbol{k}}B_{\boldsymbol{q}}+B_{\boldsymbol{k}}A_{\boldsymbol{q}}$. The second-order chemical potential satisfies the Hugenholtz–Pines theorem,[38] $$\begin{align} \mu^{(2)}={}&\varSigma_{11}^{(2)}(0)-\varSigma_{12}^{(2)}(0)=\int \frac{d \boldsymbol{q}}{(2\pi)^3}[U(0)+U({\boldsymbol{q}})]B_{\boldsymbol{q}}\\ &+\frac{1}{2}n_0\int \frac{d \boldsymbol{q}}{(2\pi)^3} U^2({\boldsymbol{q}})\Big(\frac{1}{\varepsilon^{(0)}_{\boldsymbol{q}}}-\frac{1}{\varepsilon_{\boldsymbol{q}}}\Big).~~ \tag {20} \end{align} $$ The pole equation of Green's function up to the second order is given by $$\begin{align} &\Big[p^0-\frac{1}{2}(\varSigma_{11}^{+(2)}-\varSigma_{11}^{-(2)})\Big]^2\\ ={}&\Big[\varepsilon^{(0)}_{\boldsymbol{p}}+2n_0U({\boldsymbol{p}})+\frac{1}{2}(\varSigma_{11}^{+(2)}+\varSigma_{11}^{-(2)})-\mu^{(2)}+\varSigma_{12}^{(2)}\Big]\\ &\times \Big[\varepsilon^{(0)}_{\boldsymbol{p}}+\frac{1}{2}(\varSigma_{11}^{+(2)}+\varSigma_{11}^{-(2)})-\mu^{(2)} -\varSigma_{12}^{(2)}\Big],~~ \tag {21} \end{align} $$ where $\varSigma_{11}^\pm(p)=\varSigma_{11}(\pm p)$. In the long-wavelength limit, the second-order self-energy is given by $$\begin{align} &\varSigma_{11}^{+(2)}-\varSigma_{11}^{-(2)}=2\alpha p^0,\\ &\frac{1}{2}(\varSigma_{11}^{+(2)}+\varSigma_{11}^{-(2)})-\mu^{(2)}+\varSigma_{12}^{(2)}=\beta n_0U({\boldsymbol{p}}),\\ &\frac{1}{2}(\varSigma_{11}^{+(2)}+\varSigma_{11}^{-(2)})-\mu^{(2)}-\varSigma_{12}^{(2)}=\gamma \varepsilon^{(0)}_{\boldsymbol{p}}, \end{align} $$ where $\alpha$, $\beta$, and $\gamma$ are coefficients. Since the second-order self-energies have one more power of $\sqrt{n_0}$ than the first-order self-energies, these coefficients can be treated as small quantities. To determine the leading correction to the phonon spectrum, we linearize the pole equation (21), and obtain in the long-wavelength limit $$ \varepsilon^{B}_{\boldsymbol{p}}=\sqrt{n_0U({\boldsymbol{p}})\varepsilon^{(0)}_{\boldsymbol{p}} [2+(4\alpha+2\gamma+\beta)]}.~~ \tag {22} $$ It is worth noting that each matrix element of the second-order self-energy is divergent in the long-wavelength limit, but the divergencies cancel each other in the final expression for the phonon energy in the combination $4\alpha+2\gamma+\beta$ as in the original Beliaeve theory.[35] We separate the divergent parts, $\alpha_{\rm{d}}, \beta_{\rm{d}}$, and $\gamma_{\rm{d}}$, from coefficients $\alpha, \beta$, and $\gamma$, $$\begin{align} &\alpha_{\rm{d}}=\frac{n_0^2}{2}\int d{{\boldsymbol{q}}}[U({\boldsymbol{p}})U^2({\boldsymbol{q}})]\Big(\frac{1}{\varepsilon_{{\boldsymbol{q}}}^3}\Big),\\ &\beta_{\rm{d}}=-{n_0^2}\int d{\boldsymbol{q}}[U({\boldsymbol{p}})U^2({\boldsymbol{q}})]\Big(\frac{1}{\varepsilon_{\boldsymbol{q}}^3}\Big),\\ &\gamma_{\rm{d}}=-\frac{n_0^2}{2}\int d{\boldsymbol{q}} U({\boldsymbol{p}})U^2({\boldsymbol{q}})\Big(\frac{1}{\varepsilon_{\boldsymbol{q}}^3}\Big).~~ \tag {23} \end{align} $$ which satisfy $4\alpha_{\rm{d}}+2\gamma_{\rm{d}}+\beta_{\rm{d}}=0$. The convergent parts of the coefficients are denoted by $\alpha_c, \beta_c$, and $\gamma_c$, $$\begin{align} &\alpha_c=-4\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_3(\epsilon_{\rm{dd}}), \\ &\beta_c=16\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_3(\epsilon_{\rm{dd}})+\frac{32g}{U(\boldsymbol{p})}\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_5(\epsilon_{\rm{dd}}), \\ &\gamma_c=\frac{8}{3}\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_3(\epsilon_{\rm{dd}}).~~ \tag {24} \end{align} $$ The functions $Q_3(\epsilon_{\rm{dd}})$ and $Q_5(\epsilon_{\rm{dd}})$ describe the dipolar enhancement,[15,39] given by $$\begin{align} Q_3(\epsilon_{\rm{dd}})={}&\frac{(3\epsilon_{\rm{dd}})^{3/2}}{8} \Big[(2+5y)\sqrt{1+y}\\ &+3y^2\ln\frac{1+\sqrt{1+y}}{\sqrt{y}}\Big] \\ ={}&1+\frac{3}{10}\epsilon_{\rm{dd}}^2+O(\epsilon_{\rm{dd}}^3),~~ \tag {25} \end{align} $$ $$\begin{align} Q_5(\epsilon_{\rm{dd}})={}&\frac{(3\epsilon_{\rm{dd}})^{5/2}}{48}\Big[(8+26y+33y^2)\sqrt{1+y} \\ &+15y^3\ln\frac{1+\sqrt{1+y}}{\sqrt{y}}\Big] \\ ={}&1+\frac{3}{2}\epsilon_{\rm{dd}}^2+O(\epsilon_{\rm{dd}}^4),~~ \tag {26} \end{align} $$ where $y=(1-\epsilon_{\rm{dd}})/3\epsilon_{\rm{dd}}$. The second-order chemical potential can be written explicitly in terms of these functions, $$\begin{align} \mu^{(2)}={}&\frac{32n_0g}{3}\sqrt{\frac{n_0a_{\rm{s}}^3}{\pi}}Q_5(\epsilon_{\rm{dd}})\\ &+\frac{8n_0U(0)}{3}\sqrt{\frac{n_0a_{\rm{s}}^3}{\pi}}Q_3(\epsilon_{\rm{dd}}).~~ \tag {27} \end{align} $$ Thus in the Beliaev theory of a dipolar Bose gas, the phonon energy in the long-wavelength limit is given by $$\begin{alignat}{1} \varepsilon^{B}_{\boldsymbol{p}}=\sqrt{n_0U({\boldsymbol{p}})\varepsilon^{(0)}_{\boldsymbol{p}} [2+(4\alpha_c+2\gamma_c+\beta_c)]}=v|{\boldsymbol{p}}|,~~~~~~~ \tag {28} \end{alignat} $$ with $$\begin{align} v={}&\Big(\frac{n_0}{m}\Big)^{1/2}\Big\{U({\boldsymbol{p}})\Big[1+\frac{8}{3}\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_3(\epsilon_{\rm{dd}})\Big] \\ &+16g\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_5(\epsilon_{\rm{dd}})\Big\}^{1/2}.~~ \tag {29} \end{align} $$ Comparatively, in the Bogoliubov theory, the phonon velocity is given by $v_{_{\scriptstyle \rm B}}=\sqrt{n_0 U({\boldsymbol{p}})/m}$. The main difference comes from the last term in the square root in Eq. (29), which effectively increases the s-wave interaction and changes the anisotropy of the total interaction. The Modified EGPE. The beyond-mean-field effect of a dipolar Bose gas has also been considered in the superfluid hydrodynamic and EGPE by taking account of the LHY energy. In the superfluid hydrodynamic theory, with the LHY energy considered, the sound velocity is given by[15] $$ v_{\rm{s}}=\sqrt{\frac{n}{m}\{U({\boldsymbol{p}})+16g\sqrt{n a_{\rm{s}}^3/\pi}Q_5(\epsilon_{\rm{dd}})\}},~~ \tag {30} $$ where $n$ is the total boson density. Up to the first-two orders, the total density is given by $n=n_0+n_{\rm{d}}$, where the quantum-depletion density is given by[15] $$ n_{\rm{d}}=\frac{8}{3}n_0\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_3(\epsilon_{\rm{dd}}).~~ \tag {31} $$ The last term inside the square root on the right-hand side of Eq. (30) is from the LHY energy, where the density $n$ can be replaced by the condensation density $n_0$ with negligible difference of higher orders. Thus the sound velocity given by Eq. (29) in the Beliaev theory is in agreement with the result in the superfluid hydrodynamic theory. The sound velocity of a uniform dipolar Bose gas can also be obtained from the EGPE,[24,25,40] $$ i\hbar\dot{\psi}={L}_{\rm GP}\psi, $$ where $\psi$ is the wavefunction of the condensate, and $$\begin{alignat}{1} {L}_{\rm GP}={}&-\frac{\hbar^2 \nabla^2}{2m}+\int d{\boldsymbol{r}'}V_{\rm{int}}(\boldsymbol{r}-\boldsymbol{r}')|\psi({\boldsymbol{r}'})|^2 \\ &+\frac{32}{3}g\sqrt{a_{\rm{s}}^3/\pi}Q_5(\epsilon_{\rm{dd}})|{\psi}|^3.~~ \tag {32} \end{alignat} $$ By linearizing both sides of the EGPE around the uniform solution $\psi_0=\sqrt{n_0}$, one can obtain the sound velocity given by $$ v_{_{\scriptstyle \rm GP}}=\sqrt{\frac{n_0}{m}\{U({\boldsymbol{p}})+16g\sqrt{n_0 a_{\rm{s}}^3/\pi}Q_5(\epsilon_{\rm{dd}})\}}.~~ \tag {33} $$ Compared with Eq. (30), the first term inside the square root of $v_{_{\scriptstyle \rm GP}}$ is the condensation density $n_0$, not the total density $n$. This is due to the fact that in the EGPE the mean-field energy is computed only in the condensate, whereas in superfluid hydrodynamics the mean-field energy also includes the contribution from the quantum depletion. As this difference is of the same order of the LHY term, in principle the EGPE should be modified to take it into account. It is negligible in the dilute system with $n \approx n_0$, but non-negligible in the system with significant quantum depletion. Therefore to be consistent with the Belieav and superfluid hydrodynamic theories, the EGPE should be modified as $$ i\hbar\dot{\psi}={L}'_{\rm GP}\psi, $$ where the mean-field term is modified, $$\begin{align} {L}'_{\rm GP}={}&-\frac{\hbar^2 \nabla^2}{2m}+\int d{\boldsymbol{r}'}V_{\rm{int}} ({\boldsymbol{r}-\boldsymbol{r}'}){|\psi({\boldsymbol{r}'})|}^2\\ &\cdot\Big[1+\frac{16}{9}\sqrt{ a_{\rm{s}}^3/\pi}Q_3(\epsilon_{\rm{dd}}){|\psi({\boldsymbol{r}'})|}\Big]\\ &+\frac{32}{3}g\sqrt{a_{\rm{s}}^3/\pi}Q_5(\epsilon_{\rm{dd}})|{\psi}|^3.~~ \tag {34} \end{align} $$ In Ref. [41], the mean-field term in the EGPE is also modified duo to quantum depletion, but the numerical prefactor is incorrect and not consistent with the Beliaev theory. In a nonmagnetic binary quantum droplet, this effect is of higher order and not important.[42] Application to Quantum Droplets. In the droplet state, $\epsilon_{\rm{dd}}>1$, the total interaction $U({\boldsymbol{p}})$ is attractive in the direction perpendicular to the polarization and the Bogoliubov excitation energy $\varepsilon_{\boldsymbol{p}}=\sqrt{\varepsilon^{(0)}_{\boldsymbol{p}}(2n_0U(\boldsymbol{p})+\varepsilon^{(0)}_{\boldsymbol{p}})}$ becomes imaginary in the long wavelength limit. In comparison, in the excitation energy given by Eq. (28) from the Belieav theory, the s-wave component of the interaction is effectively enhanced. Thus quantum fluctuation helps to stabilize the quantum droplet in region with $\epsilon_{\rm{dd}}>1$. However in the Belieav theory, Green's function with the Bogoliubov excitation energy is used in the computation, resulting in complex self-energy and excitation energy. The similar situation also occurs in the EGPE approach where the LHY energy of the droplet computed from the Bogoliubov theory is complex and the imaginary part is simply ignored. There is a problem in the logic of this method, as the Bogoliubov theory implies the dynamic instability of the uniform system and only the real part of its LHY energy is used to show the self-bound property of a finite-size system. Here we apply a self-consistent method to take account of the effective enhancement of the s-wave interaction up to the leading order, and show quantum fluctuations increase the window of dynamic stability in the uniform system. In the self-consistent approach, we assume that the quantum fluctuation effectively increases the s-wave coupling constant by $\delta g$ and determine it self-consistently. The effective interaction is thus given by $$ U'({\boldsymbol{p}})=g'[1+\epsilon'_{\rm{dd}}(3\rm{{\cos}}^2\theta-1)], $$ where $g'=g+\delta g$ and $g\epsilon_{\rm{dd}}=g'\epsilon'_{\rm{dd}}$. The difference between the bare and effective interactions, $-\delta g$, serves as a counter term. The rest computation procedures are essentially same as the Beliaev theory described in the above section except replacing $U$ by $U'$ and keeping track of the contribution from the counter term. In this approach, the self-energy up to the second order is then given by $$ {\rm \varSigma}'={\rm \varSigma}'^{(1)}+{\rm \varSigma}'^{(2)}-\delta g n_0 \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix},~~ \tag {35} $$ where the last term comes from the counter term, and the expression of ${\rm \varSigma}'^{(1)}$ and ${\rm \varSigma}'^{(2)}$ are the same as ${\rm \varSigma}^{(1)}$ in Eq. (15) and ${\rm \varSigma}^{(2)}$ in Eq. (19) except that $U$ is replaced by $U'$. The chemical potential in this approach is given by $\mu'=\mu'^{(1)}+\mu'^{(2)}-\delta g n_0$ with $\mu'^{(1)}$ and $\mu'^{(2)}$ defined in the same way as $\mu^{(1)}$ and $\mu^{(2)}$ except $U$ replaced by $U'$. The excitation is then determined from the pole equation of Green's function, and in the long-wavelength limit we obtain $$ \varepsilon_{\boldsymbol{p}}^{B'}=\sqrt{n_0U'({\boldsymbol{p}})\varepsilon^{(0)}_{\boldsymbol{p}} \Big[2+(4\alpha'_c+2\gamma'_c+\beta'_c)-\frac{2\delta g}{U'(\boldsymbol{p})}\Big]},~~ \tag {36} $$ where $\alpha'_c$, $\beta'_c$, and $\gamma'_c$ are defined similarly as $\alpha_c$, $\beta_c$, and $\gamma_c$. The phonon velocity in the self-consistent approach is given by $$\begin{align} v^2={}&\frac{n_0}{m}\Big\{U'({\boldsymbol{p}})\Big[1+\frac{8}{3}\sqrt{n_{0} {a'^{3}_{\rm{s}}}/\pi}Q_3(\epsilon'_{\rm{dd}})\Big] \\ &+16g'\sqrt{n_{0} a'^{3}_{\rm{s}}/\pi}Q_5(\epsilon'_{\rm{dd}})-\delta g\Big\},~~ \tag {37} \end{align} $$ where $g'=4\pi\hbar^2a'_{\rm{s}}/m$. We require that there is no s-wave enhancement to the effective interaction $U'$ which leads to the self-consistent condition $$ \frac{\delta g}{g'}=16\sqrt{\frac{n_{0} a'^{3}_{\rm{s}}}{\pi}}Q_5(\epsilon'_{\rm{dd}}).~~ \tag {38} $$ In the Bogoliubov theory, the dynamic instability occurs at $\epsilon_{\rm{dd}}=1$. In the self-consistent theory, $\epsilon'_{\rm{dd}} < 1$ at $\epsilon_{\rm{dd}}=1$, the system is dynamically stable. Thus, in general, quantum fluctuations help to stabilize the system. In experiments on $^{166}\rm{Er}$ atoms,[6] the peak density is about $n_0=33\times10^{20}\, \mathrm{m}^{-3} $ at $\epsilon_{\rm{dd}}=1.23$. We find that a uniform system with the same density is stable for $\epsilon_{\rm{dd}} < 1.35$. In experiments on $^{164}\rm{Dy}$ atoms,[2] $n_0=6.5\times10^{20}\, \mathrm{m}^{-3}$, the corresponding uniform system is stable with $\epsilon_{\rm{dd}} < 1.4$. These results of the uniform case are applicable to a three-dimensional droplet of very large size, which has not been realized experimentally so far. Preparing the droplets in a square potential well will make it easier to observe the linear phonon dispersion. In experiments,[1–6] the dipolar droplets are always finite-size systems, elongated in the polarization direction, and can be formed with larger dipolar strengths, probably due to the finite-size effect which effectively increases the kinetic energy and overcomes the interaction energy in the direction perpendicular to the polarization. In the presence of a trap, the system size increases and the linear phonon dispersion is found in an EGPE study.[43] We note that recently the sound velocity of the two-component dipolar Bose gas was derived in a hydrodynamic framework.[44] Concluding Remarks. We have gone beyond Bogoliubov approximation to obtain the phonon energy of a dipolar Bose gas in the Beliaev formalism and found that higher-order quantum fluctuations significantly increase the stability region of the quantum droplet. The quantum fluctuations effectively enhance the s-wave interaction and change the anisotropy of the total interaction. A modified EGPE in agreement with the Beliaev theory and superfluid hydrodynamics is proposed, which takes into account quantum fluctuation effects more completely. Acknowledgments. We would like to thank Z.-Q. Yu and Q. Gu for helpful discussions. This work was supported by the National Research and Development Program of China (Grant No. 2016YFA0301501). References Observing the Rosensweig instability of a quantum ferrofluidObservation of Quantum Droplets in a Strongly Dipolar Bose GasLiquid quantum droplets of ultracold magnetic atomsSelf-bound droplets of a dilute magnetic quantum liquidStriped states in a many-body system of tilted dipolesQuantum-Fluctuation-Driven Crossover from a Dilute Bose-Einstein Condensate to a Macrodroplet in a Dipolar Quantum FluidQuantum liquid droplets in a mixture of Bose-Einstein condensatesBright Soliton to Quantum Droplet Transition in a Mixture of Bose-Einstein CondensatesSelf-Bound Quantum Droplets of Atomic Mixtures in Free SpaceObservation of quantum droplets in a heteronuclear bosonic mixtureQuantum Droplets of Dipolar MixturesQuantum Droplet States of a Binary Magnetic GasEigenvalues and Eigenfunctions of a Bose System of Hard Spheres and Its Low-Temperature PropertiesQuantum Mechanical Stabilization of a Collapsing Bose-Bose MixtureBeyond mean-field low-lying excitations of dipolar Bose gasesDipolar Bose Supersolid StripesObservation of a Dipolar Quantum Gas with Metastable Supersolid PropertiesSupersolid behavior of a dipolar Bose-Einstein condensate confined in a tubeTransient Supersolid Properties in an Array of Dipolar Quantum DropletsA stripe phase with supersolid properties in spin–orbit-coupled Bose–Einstein condensatesSupersolid formation in a quantum gas breaking a continuous translational symmetrySelf-Bound Quantum Droplet with Internal Stripe Structure in One-Dimensional Spin-Orbit-Coupled Bose GasQuantum Droplets in a Mixture of Bose–Fermi SuperfluidsGround-state properties and elementary excitations of quantum droplets in dipolar Bose-Einstein condensatesSelf-bound dipolar droplet: A localized matter wave in free spaceUltradilute Low-Dimensional LiquidsMetastability of Quantum Droplet ClustersPath-Integral Monte Carlo Study on a Droplet of a Dipolar Bose–Einstein Condensate Stabilized by Quantum FluctuationDroplets of Trapped Quantum Dipolar BosonsSuperfluid Filaments of Dipolar Bosons in Free SpaceDilute dipolar quantum droplets beyond the extended Gross-Pitaevskii equationConsistent Theory of Self-Bound Quantum Droplets with Bosonic PairingTheory for self-bound states of dipolar Bose-Einstein condensatesPhonon stability and sound velocity of quantum droplets in a boson mixtureThe physics of dipolar bosonic quantum gasesMEAN-FIELD EXPANSION IN BOSE–EINSTEIN CONDENSATES WITH FINITE-RANGE INTERACTIONSGround-State Energy and Excitation Spectrum of a System of Interacting BosonsGround-state phase diagram of a dipolar condensate with quantum fluctuationsCollective Excitations of Self-Bound Droplets of a Dipolar Quantum FluidTemperature-dependent density profiles of dipolar dropletsEffective single-mode model of a binary boson mixture in the quantum droplet regionExcitation Spectrum of a Trapped Dipolar Supersolid and Its Experimental EvidenceBeyond mean-field properties of binary dipolar Bose mixtures at low temperatures
[1] Kadau H, Schmitt M, Wenzel M, Wink C, Maier T, Ferrier-Barbut I, and Pfau T 2016 Nature 530 194
[2] Ferrier-Barbut I, Kadau H, Schmitt M, Wenzel M, and Pfau T 2016 Phys. Rev. Lett. 116 215301
[3] Ferrier-Barbut I, Schmitt M, Wenzel M, Kadau H, and Pfau T 2016 J. Phys. B 49 214004
[4] Schmitt M, Wenzel M, Böttcher F, Ferrier-Barbut I, and Pfau T 2016 Nature 539 259
[5] Wenzel M, Böttcher F, Langen T, Ferrier-Barbut I, and Pfau T 2017 Phys. Rev. A 96 053630
[6] Chomaz L, Baier S, Petter D, Mark M, Wächtler F, Santos L, and Ferlaino F 2016 Phys. Rev. X 6 041039
[7] Cabrera C, Tanzi L, Sanz J, Naylor B, Thomas P, Cheiney P, and Tarruell L 2018 Science 359 301
[8] Cheiney P, Cabrera C, Sanz J, Naylor B, Tanzi L, and Tarruell L 2018 Phys. Rev. Lett. 120 135301
[9] Semeghini G, Ferioli G, Masi L, Mazzinghi C, Wolswijk L, Minardi F, Modugno M, Modugno G, Inguscio M, and Fattori M 2018 Phys. Rev. Lett. 120 235301
[10] D'Errico C, Burchianti A, Prevedelli M, Salasnich L, Ancilotto F, Modugno M, Minardi F, and Fort C 2019 Phys. Rev. Res. 1 033155
[11] Bisset R, Ardila L P N, and Santos L 2021 Phys. Rev. Lett. 126 025301
[12] Smith J C, Baillie D, and Blakie P 2021 Phys. Rev. Lett. 126 025302
[13] Lee T D, Huang K, and Yang C N 1957 Phys. Rev. 106 1135
[14] Petrov D S 2015 Phys. Rev. Lett. 115 155302
[15] Lima A R and Pelster A 2012 Phys. Rev. A 86 063609
[16] Bombin R, Boronat J, and Mazzanti F 2017 Phys. Rev. Lett. 119 250402
[17] Tanzi L, Lucioni E, Famà F, Catani J, Fioretti A, Gabbanini C, Bisset R N, Santos L, and Modugno G 2019 Phys. Rev. Lett. 122 130405
[18] Roccuzzo S M and Ancilotto F 2019 Phys. Rev. A 99 041601
[19] Böttcher F, Schmidt J N, Wenzel M, Hertkorn J, Guo M, Langen T, and Pfau T 2019 Phys. Rev. X 9 011051
[20] Li J R, Lee J, Huang W, Burchesky S, Shteynas B, Top F, Jamison A O, and Ketterle W 2017 Nature 543 91
[21] Léonard J, Morales A, Zupancic P, Esslinger T, and Donner T 2017 Nature 543 87
[22] Xiong Y and Yin L 2021 Chin. Phys. Lett. 38 070301
[23] Wang J B, Pan J S, Cui X, and Yi W 2020 Chin. Phys. Lett. 37 076701
[24] Wächtler F and Santos L 2016 Phys. Rev. A 94 043618
[25] Baillie D, Wilson R, Bisset R, and Blakie P 2016 Phys. Rev. A 94 021602
[26] Petrov D and Astrakharchik G 2016 Phys. Rev. Lett. 117 100401
[27] Kartashov Y V, Malomed B A, and Torner L 2019 Phys. Rev. Lett. 122 193902
[28] Saito H 2016 J. Phys. Soc. Jpn. 85 053001
[29] Macia A, Sánchez-Baena J, Boronat J, and Mazzanti F 2016 Phys. Rev. Lett. 117 205301
[30] Cinti F, Cappellaro A, Salasnich L, and Macrı̀ T 2017 Phys. Rev. Lett. 119 215302
[31] Böttcher F, Wenzel M et al. 2019 Phys. Rev. Res. 1 033088
[32] Hu H and Liu X J 2020 Phys. Rev. Lett. 125 195302
[33] Wang Y, Guo L, Yi S, and Shi T 2020 Phys. Rev. Res. 2 043074
[34] Gu Q and Yin L 2020 Phys. Rev. B 102 220503
[35]Beliaev S 1958 Sov. Phys.-JETP 7 299
[36] Lahaye T, Menotti C, Santos L, Lewenstein M, and Pfau T 2009 Rep. Prog. Phys. 72 126401
[37] Schützhold R, Uhlmann M, Xu Y, and Fischer U R 2006 Int. J. Mod. Phys. B 20 3555
[38] Hugenholtz N and Pines D 1959 Phys. Rev. 116 489
[39] Bisset R, Wilson R, Baillie D, and Blakie P 2016 Phys. Rev. A 94 033619
[40] Baillie D, Wilson R, and Blakie P 2017 Phys. Rev. Lett. 119 255302
[41] Aybar E and Oktel M 2019 Phys. Rev. A 99 013620
[42] Xiong Y and Yin L 2022 Phys. Rev. A 105 053305
[43] Natale G, van Bijnen R, Patscheider A, Petter D, Mark M, Chomaz L, and Ferlaino F 2019 Phys. Rev. Lett. 123 050402
[44] Pastukhov V 2017 Phys. Rev. A 95 023614