⇦Chinese Physics Letters, 2023, Vol. 40, No. 6, Article code 066701 Hydrodynamics of a Multi-Component Bosonic Superfluid Fan Zhang (张帆) and Lan Yin (尹澜)* Affiliations School of Physics, Peking University, Beijing 100871, China Received 18 March 2023; accepted manuscript online 4 May 2023; published online 16 May 2023 *Corresponding author. Email: yinlan@pku.edu.cn Citation Text: Zhang F and Yin L 2023 Chin. Phys. Lett. 40 066701    Abstract We obtain the superfluid hydrodynamic equations of a multi-component Bose gas with short-ranged interactions at zero temperature under the local equilibrium assumption and show that the quantum pressure is generally present in the nonuniform case. Our approach can be extended to systems with long-range interactions such as dipole-dipole interactions by treating the Hartree energy properly. For a highly symmetric superfluid, we obtain the excitation spectrum and show that except for the density phonon, all other excitations are all degenerate. The implication of our results is discussed. DOI:10.1088/0256-307X/40/6/066701 © 2023 Chinese Physics Society Article Text About ninety years ago, shortly after the discovery of superfluidity,[1] Tisza[2] and Landau[3] proposed a two-fluid model to describe the hydrodynamics of a bosonic superfluid. At zero temperature, the superfluid flow is described by a set of superfluid hydrodynamic equations.[4] In this phenomenological approach based on macroscopic conservation laws,[5-9] as well as in the effective-theory approach,[10] the superfluid acceleration is determined by the local chemical potential. For a dilute Bose gas, the superfluid dynamics are often studied by the Gross–Pitaevskii equation (GPE)[11,12] in which every boson is assumed to be in the condensate. The superfluid hydrodynamic equations can be derived from the GPE,[13-16] but in this approach, the superfluid acceleration is determined not only by the local chemical potential, but also by the quantum pressure which depends on the density gradient. This microscopic approach[17-29] is valid in the dilute region when the system is weakly interacting. The hydrodynamic equations of the bosonic superfluid can also be derived from the Schrödinger equation for short-ranged interactions by the interaction-radius expansion,[30,31] where the quantum pressure is found in the mean-field level, but higher order terms are inconsistent with the Lee–Huang–Yang (LHY) energy. Whether or not the quantum pressure is present generally in the superfluid remains a critical issue. In this Letter, we obtain the general form of the superfluid hydrodynamic equations of a multi-component bosonic system with short-ranged interactions in the path-integral formalism under the local equilibrium assumption (LEA) and shows that the quantum pressure generally exists in a nonuniform superfluid. The quantum-pressure term in the superfluid hydrodynamic equations is necessary to recover the free-particle energy at large momentum. Our approach can be extended to boson systems with long-ranged interactions by treating the Hartree energy as an external potential. The superfluid hydrodynamic equations fully describe the motion of atoms in and outside the condensate in contrast to the GPE-type equations which are in terms of the condensate wavefunctions. The Andreev–Bashkin effect of a multi-component superfluid, the drag between two components with different superfluid velocities, can be described in our approach. We further study the case of a highly symmetric superfluid and obtain its excitation spectrum, showing that except for the density phonon, all other excitations are all degenerate. The implication to quantum droplets is discussed. Superfluid Hydrodynamic Equations under LEA. We first consider a multi-component bosonic system with short-ranged interactions. Its action is given by \begin{align} S=\,&\int dt d {\boldsymbol r} \sum_{\sigma} \varPsi^*_{\sigma} ({\boldsymbol r},t)\Big[i \hbar \partial_t+\frac{\hbar^2\nabla^{2}}{2m_{\sigma}}-V_\sigma({\boldsymbol r})\Big]\varPsi_{\sigma}({\boldsymbol r},t)\notag\\ &+S_{\rm int} , \tag {1} \\ S_{\rm int}=\,&-\frac{1}{2}\int dt d {\boldsymbol r} d {\boldsymbol r}' \sum_{\sigma {\sigma}'} \varPsi^*_{\sigma} ({\boldsymbol r},t)\varPsi^*_{\sigma'}({\boldsymbol r}',t) U_{{\sigma} {\sigma}'}({\boldsymbol r}-{\boldsymbol r}') \notag\\ &\times\varPsi_{\sigma'}({\boldsymbol r}',t)\varPsi_{\sigma}({\boldsymbol r},t), \tag {2} \end{align} where $\varPsi_{\sigma}$ and $m_{\sigma}$ are the boson field and mass of the $\sigma$ component, $U_{{\sigma} {\sigma}'}({\boldsymbol r}-{\boldsymbol r}')$ is the short-ranged interaction, and $V_\sigma({\boldsymbol r})$ is the external potential. The boson field $\varPsi_{\sigma}$ can be written in terms of the density $\rho_{\sigma}$ and phase $\varPhi_{\sigma}$ as $\varPsi_{\sigma}=\sqrt{\rho_{\sigma}}e^{i\varPhi_{\sigma}}$. In the superfluid phase, the density and phase fields can be separated into their mean values and fluctuations, \begin{align} \rho_{\sigma}=n_{\sigma}+\delta n_{\sigma},~\varPhi_{\sigma}=\phi_{\sigma}+\delta \phi_{\sigma} , \tag {3} \end{align} where $n_{\sigma}$ is the average density of the $\sigma$ component, $\phi_{\sigma}$ is the mean phase of the $\sigma$ component, and $\delta n_{\sigma}$ and $\delta \phi_{\sigma}$ are fluctuations in the density and phase.[7,16,32,33] For a homogeneous system with translational invariance, the density $n_{\sigma}$ is the superfluid density. When the translational invariance is broken, the superfluid density can be smaller than $n_{\sigma}$ due to Leggett's bound.[34-36] In the extreme case such as the Mott-insulator phase of a Bose gas in the optical lattice, the superfluid fraction can vanish. In the following, we focus on the case that the external potential $V$ is not strong and slowly varies in space and time, so that the system is in the superfluid phase. We consider the dynamical process with temporal and spatial scales of variances much larger than the intrinsic scales of the system so that the LEA can be applied. In this case, the system can be divided into many small blocks in space and time with scales much smaller than those of the spatial and temporal variance and much larger than the intrinsic scales, e.g., healing lengths, interparticle distance, and the range of the interaction potential in space, local relaxation time and period of the boson de Broglie wave. The subsystem within each block can be approximated as a macroscopic and homogeneous system in its local equilibrium, with constant densities and linear phases. In the path integral, the contribution to the time-evolution matrix element from each block is given by \begin{align} \mathcal{Z}_j=\int_j \mathcal{D}\varPsi \mathcal{D}\varPsi^* e^{\frac{i}{\hbar}S_j}, \tag {4} \end{align} where $j=(j_x,j_y,j_z,j_t)$ is the overall index of the 4-dimensional block, the subscript $j$ of the integral stands for the integration over the fields within the block, and $S_j$ is the action of the block. Within the block, the average density $n_{\sigma}({\boldsymbol r}, t)$ varies linearly in space and time, whereas this variance is very small due to the size of the block, and the average phase behaves similarly. Thus the action of the block $S_j$ is approximately given by the action of a uniform system $S_j^U$ plus a term $S_{0j}$ containing derivatives of the average density and phase, $S_j \approx S_j^U+S_{0j}$, where \begin{align} S_{0j} \approx - \varDelta_j \sum_{\sigma}n_{j\sigma}\Big[\hbar\partial_t\phi_{j\sigma}+\frac{\hbar^2}{2m_{\sigma}}\frac{|\nabla n_{j\sigma}|^2}{4n_{j\sigma}}\Big] , \tag {5} \end{align} $n_{j\sigma}=n_{\sigma}({\boldsymbol r}_j,t_j)$ is the mean density, $\phi_{j\sigma}=\phi_{\sigma}({\boldsymbol r}_j,t_j)$ is mean phase, $({\boldsymbol r}_j,t_j)$ and $\varDelta_j$ are the center and volume of the block, and $S_j^U$ is the action of the corresponding a uniform system within the same block with average density $n_{j\sigma}$ and phase $\phi_{j\sigma}$, and phase gradient $\nabla \phi_{j\sigma}$. Thus in LEA, integrating out the boson fields is equivalent to putting the time-evolution operator on the lowest-energy state of a uniform system, \begin{align} \mathcal{Z}_j \sim e^{\frac{i}{\hbar}[S_{0j}-\varDelta_j\mathcal{E}_j ]}, \tag {6} \end{align} where $\mathcal{E}_j$ is the energy density of the lowest-energy state of the uniform system with constant densities $n_{j\sigma}$ and phase gradient $\nabla \phi_{j\sigma}$. After performing the integration in all the blocks, we obtain an effective action for the whole system \begin{align} S_{\rm eff}=\,&\int d t \int d {\boldsymbol r} \Big\{ -\sum_{\sigma} n_{\sigma}\Big[\hbar \partial_t\phi_{\sigma }+V_{\sigma}\notag\\ & +\frac{\hbar^2}{2m_\sigma}\Big(\frac{|\nabla n_{\sigma}|^2}{4n_{\sigma}}+|\nabla \phi_{\sigma}|^2\Big)\Big]-\mathcal{E}_{\rm I}\Big\}, \tag {7} \end{align} where $\mathcal{E}_{\rm I}$ is energy density due to interaction which shall be obtained from studies on uniform systems. The equation of motion can be obtained by the variational conditions \begin{align} \frac{\delta S_{\rm eff}}{\delta \phi_{\sigma}}=0,~~\frac{\delta S_{\rm eff}}{\delta n_{\sigma}}=0, \tag {8} \end{align} leading to the superfluid hydrodynamic equations \begin{align} &\partial_t n_{\sigma}=-\nabla\cdot(n_{\sigma}{\boldsymbol v}_{\sigma}+{\boldsymbol j}_{\sigma}'), \notag\\ &m_{\sigma}\partial_t {\boldsymbol v}_{\sigma}=-\nabla\Big(\frac{1}{2}m_{\sigma} v_{\sigma}^2+\frac{\partial\mathcal{E}_{\rm I}}{\partial n_{\sigma}}+V_{\sigma}-\frac{\hbar^2\nabla^2\sqrt{n_{\sigma}}}{2m_{\sigma}\sqrt{n_{\sigma}}}\Big), \tag {9} \end{align} where ${\boldsymbol v}_{\sigma}=\hbar \nabla\phi_{\sigma}/m_{\sigma}$ is the $\sigma$-component superfluid velocity, ${\boldsymbol j}_{\sigma}'=(j'_{\sigma x},j'_{\sigma y},j'_{\sigma z})$, $j'_{\sigma\alpha}=\partial \mathcal{E}_{\rm I}/(\hbar \partial \phi'_{\sigma\alpha}) $, and $\phi'_{\sigma\alpha}=\partial_{\alpha} \phi_{\sigma}$, $\alpha=x,y,z$. The anomalous current ${\boldsymbol j}_{\sigma}'$ is related to the Andreev–Bashkin effect of a multi-component superfluid, which will be discussed in the latter part. In our approach, the superfluid hydrodynamic equations are obtained without the dilute condition and are applicable to strongly interacting bosonic systems as well as those with nonnegligible beyond-mean-field effects. The quantum-pressure term shows up as the last term on the right-hand side of Eq. (9), consistent with the quantum potential in systems with nonuniform probability densities.[37] The Single-Component Case. To understand the significance of the quantum pressure, and compare the hydrodynamics equations from the GPE, we first look into the single-component case. Due to Galilean invariance, the energy density due to interaction, $\mathcal{E}_{\rm I}(n,\nabla \phi)$, is exactly the energy density of uniform ground-state $\mathcal{E}_0(n)$. The superfluid hydrodynamic equations can be simplified to the familiar form, \begin{align} &\partial_t n=-\nabla\cdot(n{\boldsymbol v}), \notag\\ &m\partial_t {\boldsymbol v}=-\nabla\Big(\frac{1}{2}m v^2+\mu-\frac{\hbar^2\nabla^2\sqrt{n}}{2m\sqrt{n}}\Big), \tag {10} \end{align} where $\mu=V+\partial\mathcal{E}_0/\partial n$ is the local chemical potential. If we define the order parameter $\psi=\sqrt{n}e^{i\phi}$, it is easy to obtain its equation of motion from Eq. (10), \begin{align} i\hbar\partial_t \psi({\boldsymbol r},t)= \Big[-\frac{\hbar^2\nabla^{2}}{2\,m}+\mu \Big]\psi({\boldsymbol r},t) , \tag {11} \end{align} which looks similar to the GPE, but it is obtained without the assumption that every atom is in the condensate. Following a standard textbook treatment, the phonon energy of a uniform superfluid can be obtained from Eq. (10), \begin{align} \omega_{\boldsymbol k}=\sqrt{\epsilon_{\boldsymbol k}(\epsilon_{\boldsymbol k}+2\mu)}, \tag {12} \end{align} where $\epsilon_{\boldsymbol k}=\hbar^2 k^2/2\,m $ is the kinetic energy of a free boson. In contrast, if the quantum pressure term is missing, the phonon energy is just linearly dispersed, \begin{align} \omega'_{\boldsymbol k}=\sqrt{2\mu\epsilon_{\boldsymbol k}}. \tag {13} \end{align} Thus the quantum pressure is responsible for recovering the free-particle kinetic energy at large momentum. Although the order parameter and the condensate wavefunction have the same phase, it relates to the total density $n$, not just the condensate density $n_0$. In the dilute limit, to the leading order of the gas parameter, the superfluid density is equal to the condensate density, $n \approx n_{0}$, and Eq. (11) is the same as the GPE. It is a legitimate question whether or not by modifying GPE, the higher-order effects can be captured. To answer this question, we rewrite the continuity equation as \begin{align} \partial_t n_{0}+\partial_t \delta n=-\nabla\cdot({\boldsymbol j}_{0}+\delta {\boldsymbol j}), \tag {14} \end{align} where $\delta n= n-n_{0}$ is the quantum depletion and is a function of $n_0$ in LEA, ${\boldsymbol j}_{0}=n_{0} {\boldsymbol v}$ is the condensate current, and $\delta {\boldsymbol j}=\delta n{\boldsymbol v}$ is the current of the quantum depletion. To the leading order, the continuity equation holds for the condensate, $\partial_t n_{0}=-\nabla\cdot {\boldsymbol j}_{0}$, which can be put into the second term on the left-hand side of Eq. (14) for the next order, and we obtain \begin{align} \partial_t n_{0}=-\Big(1-\frac{d \delta n}{d n_{0}}\Big)\nabla\cdot{\boldsymbol j}_{0}-\nabla \frac{ \delta n}{n_{0}} \cdot {\boldsymbol j}_{0}. \tag {15} \end{align} Equation (15) together with Eq. (10) cannot be rewritten in the form of a GPE-type equation, \begin{align} i\hbar\partial_t \psi_{0}({\boldsymbol r},t)=\Big[-\frac{\hbar^2\nabla^{2}}{2\,m}+\mu'\Big]\psi_{0}({\boldsymbol r},t), \tag {16} \end{align} for any modified chemical potential $\mu'$. Thus we conclude that the GPE-type equations cannot completely describe the effect of the quantum depletion. The Andreev–Bashkin Effect of a Multi-Component Superfluid. Andreev and Bashkin predicted that in a 3He–4He mixture the supercurrent of one component will drag the other component into superflow.[38] Although it has not been observed so far, there are theoretical proposals[39-41] that such effect is also present in a two-component Bose superfluid. Here we show that the Andreev–Bashkin effect of a multi-component Bose superfluid can be generally described by the superfluid hydrodynamic equations given by Eq. (9). We consider the case with small supervelocities in the linear-response regime so that the interaction-energy density can be expanded to the quadratic order in phase gradients, \begin{align} &\,\mathcal{E}_{\rm I} (\{n_{\sigma}\},\{\nabla\phi_{\sigma}\})\notag\\ \approx&\, \mathcal{E}_0(\{n_{\sigma}\})+\frac{1}{2}\sum_{\sigma,\sigma'}\chi_{\sigma\sigma'}\nabla \phi_{\sigma}\cdot\nabla \phi_{\sigma'}, \tag {17} \end{align} where $\mathcal{E}_0(\{n_{\sigma}\})$ is the energy density of the uniform ground state, and the expansion coefficient $\chi_{\sigma\sigma'}$ is related to the static current-current correlation function,[42] $\chi_{\sigma\sigma'}=\chi_{\sigma'\sigma}$. Note that we assume the inversion symmetry which eliminates the transverse coupling between phase gradients in the quadratic order. For the case of a uniform system with the same superfluid velocity ${\boldsymbol v}$ for all the components, the Galilean invariance holds \begin{align} \mathcal{E}_{\rm I}\Big(\{n_{\sigma}\},\Big\{\nabla\phi_{\sigma}=\frac{m_{\sigma}}{\hbar}{\boldsymbol v}\Big\}\Big)= \mathcal{E}_0(\{n_{\sigma}\}), \tag {18} \end{align} leading to the identity \begin{align} \sum_{\sigma,\sigma'}\chi_{\sigma\sigma'}\nabla \phi_{\sigma}\cdot\nabla \phi_{\sigma'} =\sum_{\sigma,\sigma'} f_{\sigma\sigma'} |{\boldsymbol v}_{\sigma}-{\boldsymbol v}_{\sigma'}|^2, \tag {19} \end{align} where the coefficient $f_{\sigma\sigma'}$ depends on densities. In this case, the continuity equation is simply given by \begin{align} \partial_t n_{\sigma}=-\nabla\cdot n_{\sigma}{\boldsymbol v}, \tag {20} \end{align} which yields \begin{align} \sum_{\sigma'} m_{\sigma'}\chi_{\sigma\sigma'}=0. \tag {21} \end{align} The total current is given by ${\boldsymbol j}=\sum_{\sigma} n_{\sigma}{\boldsymbol v}$, in agreement with the two-component case.[42] The superfluid hydrodynamic equations can be further written as \begin{align} \partial_t n_{\sigma}=\,&-\nabla\cdot\Big(n_{\sigma}{\boldsymbol v}_{\sigma}+\sum_{\sigma'} \frac{m_{\sigma'}}{\hbar^2} \chi_{\sigma \sigma'}{\boldsymbol v}_{\sigma'}\Big) , ~~~~~~ \tag {22a}\\ m_{\sigma}\partial_t {\boldsymbol v}_{\sigma}=\,&-\nabla\Big(\frac{1}{2}m_{\sigma} v_{\sigma}^2+\mu_{\sigma}-\frac{\hbar^2\nabla^2\sqrt{n_{\sigma}}}{2m_{\sigma}\sqrt{n_{\sigma}}} \notag\\ &+\sum_{\sigma_1,\sigma_2}\frac{1}{2}({\boldsymbol v}_{\sigma_1}-{\boldsymbol v}_{\sigma_2})^2\nabla \frac{\partial f_{\sigma_1\sigma_2}}{\partial n_{\sigma}}\Big),\,~~~~~~ \tag {22b} \end{align} where $\mu_{\sigma}=\partial \mathcal{E}_0/\partial n_{\sigma}+V_{\sigma}$ is the local chemical potential. In the last term on the right-hand side of Eq. (22b), the derivative $\partial f_{\sigma_1\sigma_2}/\partial n_{\sigma}$ should be significant only for $\sigma=\sigma_1$ or $\sigma=\sigma_2$, because the density change in other components barely affects the coupling between velocities of these two components. Therefore in general the drag force is present when there is a finite velocity difference between two components. The continuity Eq. (22a) shows that generally the current of the $\sigma$-component is different from $n_{\sigma}{\boldsymbol v}_{\sigma}$ due to the anomalous current, ${\boldsymbol j}_{\sigma}'=\sum_{\sigma'} m_{\sigma'} \chi_{\sigma \sigma'}{\boldsymbol v}_{\sigma'}/\hbar^2$. Boson Systems with Long-Ranged Interactions. For a uniform many-body system with long-ranged interactions, the Hartree energy is always divergent in the thermodynamic limit, but this divergence is artificial for physical systems. In solids, the electrostatic energy from the ions cancels out the Hartree energy of the electrons, as demonstrated in the jellium model.[43] In ultracold quantum gases, although the dipole-dipole interaction (DDI) is long-ranged, the experimental systems are always finite in size. For these finite systems, the Hartree energy as well as the rest interaction energies can be expressed in density functionals, as in electron systems.[44] In the following, by properly treating the Hartree energy, we derive the superfluid hydrodynamic equations of bosons with long-ranged interactions within LEA. We consider a multi-component bosonic system with the action given by Eq. (1), except that the interactions are now long-ranged. We first single out the Hartree energy given by \begin{align} E_H=\frac{1}{2}\sum_{\sigma,\sigma'}\int d {\boldsymbol r} d {\boldsymbol r}' U_{{\sigma} {\sigma}'}({\boldsymbol r}-{\boldsymbol r}')n_{\sigma'}({\boldsymbol r}') n_\sigma({\boldsymbol r}) .~~ \tag {23} \end{align} The rest steps are essentially the same as in the above section. Under LEA, we divide the system into many small blocks. Neighboring blocks are now coupled by the Hartree energy. In each block, the Hartree energy is treated as a constant potential, and fluctuations can be integrated out locally. Thus we obtain the effective action given by \begin{align} S_{\rm eff}=\,&\int d t \int d {\boldsymbol r} \Big\{ \sum_{\sigma}n_{\sigma}\Big[-\hbar \partial_t \phi_{\sigma}-V_\sigma\notag\\ & -\frac{\hbar^2}{2m_{\sigma}}\Big(\frac{|\nabla n_{\sigma}|^2}{4n_{\sigma}}+|\nabla \phi_{\sigma}|^2\Big)\notag\\ &+\frac{1}{2} \sum_{\sigma'}\int d {\boldsymbol r}' U_{{\sigma} {\sigma}'}({\boldsymbol r}-{\boldsymbol r}')n_{\sigma'}({\boldsymbol r}')\Big]-\mathcal{E}_{\rm I}\Big\},~~ \tag {24} \end{align} where $\mathcal{E}_{\rm I}$ is now the interaction-energy density of the lowest-energy state for uniform densities and phase gradients with the Hartree energy excluded, and can be calculated in the thermodynamic limit by adding a proper background potential to cancel out the Hartree energy as in the jellium model. The superfluid hydrodynamic equations can be obtained from the equation of motion as follows: \begin{align} &\partial_t n_{\sigma}=-\nabla\cdot(n_{\sigma}{\boldsymbol v}_{\sigma}+{\boldsymbol j}_{\sigma}'), \notag\\ &m_{\sigma}\partial_t {\boldsymbol v}_{\sigma}=-\nabla\Big(\frac{1}{2}m_{\sigma} v_{\sigma}^2+\frac{\partial\mathcal{E}_{\rm I}}{\partial n_{\sigma}}+V'_{\sigma}-\frac{\hbar^2\nabla^2\sqrt{n_{\sigma}}}{2m_{\sigma}\sqrt{n_{\sigma}}}\Big),~~ \tag {25} \end{align} where the modified potential is given by \begin{align} V'_{\sigma}({\boldsymbol r},t)=V_{\sigma}({\boldsymbol r},t)+\sum_{\sigma'}\int d {\boldsymbol r}'U_{{\sigma} {\sigma}'}({\boldsymbol r}-{\boldsymbol r}')n_{\sigma'}({\boldsymbol r}',t).\notag \end{align} Thus the superfluid hydrodynamic equations have the same form as those in the case with the short-ranged interactions except for the treatment of the Hartree energy. Excitations of a Symmetric Superfluid. The dynamics of a multi-component boson superfluid are rather complex. In this section, we focus on a highly symmetric case where all $N$ components have the same density $n$ and mass $m$, all the intraspecies interactions are the same, and all the interspecies interactions are the same. In the absence of superfluid currents, the excitations in a uniform system can be described by linearizing the superfluid hydrodynamic Eq. (22), \begin{align} &\partial_t \delta n_{\sigma}=-\sum_{\sigma'}D_{\sigma\sigma'}\nabla \cdot \delta {\boldsymbol v}_{\sigma'}, \notag \\ &\partial_t \delta {\boldsymbol v}_{\sigma}=-\sum_{\sigma'}[F_{\sigma\sigma'} \nabla \delta n_{\sigma'}-G_{\sigma\sigma'} \nabla^3 \delta n_{\sigma'}],~~ \tag {26} \end{align} where $\delta n_{\sigma}$ and $\delta {\boldsymbol v}_{\sigma}$ are the fluctuations in superfluid density and velocity, $D_{\sigma\sigma'}=n\delta_{\sigma\sigma'}+m \chi_{\sigma\sigma'}/\hbar^2$, $F_{\sigma\sigma'}=\partial \mu_{\sigma}/m\partial n_{{\sigma}'}$, and $G_{\sigma\sigma'}=\delta_{\sigma\sigma'}\hbar^2/(4m^2n)$. We look for plane-wave excitations with $\delta n_{\sigma}=X_{\sigma}e^{i{\boldsymbol q} \cdot {\boldsymbol r}-i\omega_q t}$ and $\delta {\boldsymbol v} _{\sigma}= Y_{\sigma} {\boldsymbol q} e^{i{\boldsymbol q} \cdot {\boldsymbol r}-i\omega_q t}$. From Eq. (26), we obtain \begin{align} &\omega_q X_{\sigma}=\sum_{\sigma'}D_{\sigma\sigma'} q^2 Y_{\sigma'}, \notag \\ & \omega_q Y_{\sigma}=\sum_{\sigma'}\big[F_{\sigma\sigma'}+G_{\sigma\sigma'} q^2\big] X_{\sigma'},~~ \tag {27} \end{align} leading to the equation for ${\boldsymbol X}$, \begin{align} \omega_q^2 {\boldsymbol X}=q^2 {\boldsymbol D} ({\boldsymbol F}+q^2 {\boldsymbol G}) {\boldsymbol X}.~~ \tag {28} \end{align} The excitation frequency $\omega_q$ can be solved from the secular equation \begin{align} |q^2 {\boldsymbol D} ({\boldsymbol F}+q^2 {\boldsymbol G})- \omega_q^2 {\boldsymbol I}|=0.~~ \tag {29} \end{align} In Eq. (29), the matrix ${\boldsymbol G}$ comes from the quantum pressure and is diagonal. Matrixes ${\boldsymbol D}$, ${\boldsymbol G}$ and their product ${{\boldsymbol D}} ({\boldsymbol F}+q^2 {{\boldsymbol G}})$ are highly symmetric and share the same property that in each matrix all the diagonal terms are the same, and all the off-diagonal terms are the same. From Eq. (29), we obtain that in the long-wavelength limit the density-phonon frequency is given by \begin{align} \omega_{q} \approx c_{\rm p} q.~~ \tag {30} \end{align} All the other excitations are degenerate, with the excitation energy in the long-wavelength limit given by \begin{align} \omega'_{q} \approx c_m q ,~~ \tag {31} \end{align} where $c_{\rm p}^2=[(N-1)D_{\rm o}+D_{\rm d}][(N-1)F_{\rm o}+F_{\rm d}]$, $c_m^2=(D_{\rm d}-D_{\rm o})(F_{\rm d}-F_{\rm o})$; $D_{\rm d}$ and $D_{\rm o}$ are the diagonal and off-diagonal matrix-elements of ${\boldsymbol D}$, respectively; and $F_{\rm d}$ and $F_{\rm o}$ are the diagonal and off-diagonal matrix-elements of ${\boldsymbol F}$, respectively. The density-phonon speed $c_{\rm p}$ does not depend on the susceptibility $\chi_{\sigma \sigma'}$ based on Eq. (21), consistent with the two-component case.[42] For a highly-symmetric and weakly-interacting Bose gas, to the lowest order in the gas parameter, the susceptibility $\chi_{\sigma \sigma'}$ can be ignored, $D_{\rm o} \approx 0$, $D_{\rm d} \approx n$, $F_{\rm d} \approx g_{1}/m$, $F_{\rm o} \approx g_{2}/m$, and the phonon speeds are given by, \begin{align} &c_{\rm p} \approx \sqrt{[g_1+(N-1)g_2]n/m},\notag \\ &c_m \approx \sqrt{(g_1-g_2)n/m},~~ \tag {32} \end{align} where $g_1$ and $g_2$ are the intra- and inter-species coupling constants, respectively. We can infer from Eq. (32) that the mean-field stability condition of the superfluid state is given by $g_1>g_2$ and $g_1+(N-1)g_2>0$. Our results can be tested if the highly symmetric Bose gas can be realized in experiments in the future. Although for a general asymmetric Bose gas the degeneracy in the excitation spectrum is broken, its excitation spectrum can still be obtained from the superfluid hydrodynamic equations. Discussions and Conclusion. Quantum droplets of ultracold atoms are formed in the mean-field unstable region and stabilized by quantum fluctuations. They are perfect examples to demonstrate the beyond-mean-field effects. So far most experiments on quantum droplets are performed on nonmagnetic binary boson mixtures, such as the homonuclear mixture of $^{39}\rm{K}$,[45-47] and the heteronuclear $^{39}\rm{K}$–$^{87}\rm{Rb}$,[48] and the single-component dipolar Bose gas, such as $^{164}\rm{Dy}$[49-53] and $^{166}\rm{Er}$[54] atoms. In the binary boson mixture, the mean-field energy is small and attractive, and the system is stabilized by the LHY energy.[55] In the Bogoliubov theory, the LHY energy has not only a dominant real part but also a small imaginary part due to phonon instability which sparks the research on the ground state of the droplet.[56-58] In a Beliav approach,[58] it was found that higher-order fluctuations restore the phonon stability removing the imaginary part of the LHY energy, which was also confirmed in a path-integral approach.[59] For this binary droplet without phase separation,[60] the density is approximately given by the condensate, $n_\sigma \approx n_{0\sigma}$, and the interaction-energy density of the droplet is dominated by the mean-field and LHY energies,[55] \begin{align} \mathcal{E}_{\rm I} \approx \sum_{\sigma{\sigma}'}\frac{g_{\sigma{\sigma}'}}{2}n_{\sigma}n_{{\sigma}'} +\frac{8m^{3/2}}{15\pi^2\hbar^3}(g_{11}n_1+g_{22}n_2)^{5/2},~~ \tag {33} \end{align} where $m=m_1=m_2$, and $a_{\sigma\sigma'}=\frac{mg_{\sigma\sigma'}}{4\pi \hbar^2}$ is the scattering length between a $\sigma$-component atom and a $\sigma'$-component atom. From Eq. (11), we obtain the EGPE of the two components as in Ref. [48], \begin{align} i\hbar\partial_t \psi_{1}({\boldsymbol r},t)= \Big[-\frac{\hbar^2\nabla^{2}}{2\,m}+\frac{\partial{\mathcal E}_{\rm I}}{\partial n_1}\Big]\psi_{1}({\boldsymbol r},t),\notag \\ i\hbar\partial_t \psi_{2}({\boldsymbol r},t)= \Big[-\frac{\hbar^2\nabla^{2}}{2\,m}+\frac{\partial{\mathcal E}_{\rm I}}{\partial n_2}\Big]\psi_{2}({\boldsymbol r},t),~~ \tag {34} \end{align} where $\psi_{\sigma}=\sqrt{n_{\sigma}}e^{i\phi_\sigma}$. For the binary quantum droplet, $\delta g =g_{12}+\sqrt{g_{11}g_{22}} < 0$, the densities satisfy $n_1/n_2\approx\sqrt{g_{22}/g_{11}}$ and $\frac{\partial{\mathcal E}_{\rm I}}{\partial n_1}=\frac{\partial{\mathcal E}_{\rm I}}{\partial n_2}=\mu$,[59] and in the solution the wavefunctions of the two components are proportional to each other. Thus EGPE of the binary quantum droplet is reduced to a single-mode equation \begin{align} i\hbar\partial_t \psi=\,&\Big[-\frac{\hbar^2\nabla^{2}}{2\,m}+\frac{2\sqrt{g_{11}g_{22}}\delta g|\psi|^2}{(\sqrt{g_{11}}+\sqrt{g_{22}})^2}\notag\\ &+\frac{4m^{3/2}}{3\pi^2\hbar^3}(g_{11}g_{22})^{\frac{5}{4}}{|\psi|}^3\Big]\psi,~~ \tag {35} \end{align} where $\psi=\sqrt{n_{\rm tot}}e^{i\phi}$ and $n_{\rm tot}=n_1+n_2$. This equation is consistent with the result in Ref. [55]. A similar situation occurs in a single-component dipolar Bose gas, where the quantum droplet is formed in the region with the DDI strength larger than that of the repulsive s-wave interaction, and stabilized by the LHY energy as demonstrated from EGPE.[61-64] In the Bogoliubov theory, the LHY energy of the dipolar droplet also has an imaginary part due to phonon instability in certain propagating directions. This imaginary energy is also removed by higher-order fluctuations.[65] In this case, the EGPE can be also constructed from the superfluid hydrodynamic equations. In both cases, the EGPE is successful in the dilute region by taking into account the LHY energy which is the dominant beyond-mean-field effect. In general, however, it is more advantageous to study the superfluid hydrodynamic equations for their broader range of validity, especially when the quantum depletion is significant. In conclusion, we obtain the superfluid hydrodynamic equations of a multi-component bosonic system in the path-integral formalism under the local equilibrium assumption and show that the quantum pressure term is generally present. Our method is valid for nonuniform and strongly interacting systems with short-ranged interactions as well as with long-ranged interactions. The superfluid hydrodynamic equations are superior to the GPE-type equations with the condensate wavefunctions as variables, as the latter cannot fully describe the dynamics of quantum depletion. Our approach provides a general description of the Andreev and Bashkin effect of a multi-component Bose superfluid. In the case of a symmetric superfluid, we find that except for the density phonon all other excitations are degenerate. The implications on quantum droplets are discussed. Acknowledgments. We would like to thank Z. Q. Yu for his helpful discussion. References Viscosity of Liquid Helium below the λ-PointTransport Phenomena in Helium IITheory of the Superfluidity of Helium IINonlinear Schrödinger equation and superfluid hydrodynamicsTwo-fluid hydrodynamic modes in a trapped superfluid gasBeyond mean-field properties of binary dipolar Bose mixtures at low temperaturesLifshitz hydrodynamicsStructure of a quantized vortex in boson systemsCollective Excitations of a Trapped Bose-Condensed GasQuantized hydrodynamic model and the dynamic structure factor for a trapped Bose gasHydrodynamic quantization approach to Bose-Einstein condensationFluidity aspects of Bose-Einstein condensed systemsTwo-Fluid Dynamics for a Bose-Einstein Condensate out of Local Equilibrium with the NoncondensateLandau-Khalatnikov two-fluid hydrodynamics of a trapped Bose gasBose-Condensed Gases at Finite TemperaturesCollective modes and generation of a new vortex in a trapped Bose gas at finite temperatureSuperfluid Bose-Fermi mixture from weak coupling to unitarityFrom Bose condensation to quantum gravity and backGeneralized quantum hydrodynamics of a trapped dilute Bose gasCounter-flow instability of a quantum mixture of two superfluidsProblem with the single-particle description and the spectra of intrinsic modes of degenerate boson-fermion systemsHydrodynamics of the atomic Bose–Einstein condensate beyond the mean-field approximationQuantum hydrodynamics of the spinor Bose–Einstein condensate at non-zero temperaturesLow dimensional Bose gasesCan a Solid Be "Superfluid"?On the Superfluid Fraction of an Arbitrary Many-Body System at T=0Superfluid fraction in an interacting spatially modulated Bose-Einstein condensateA Suggested Interpretation of the Quantum Theory in Terms of "Hidden" Variables. INondissipative drag of superflow in a two-component Bose gasAndreev–Bashkin effect in superfluid cold gases mixturesLinear response study of collisionless spin dragSelf-Consistent Equations Including Exchange and Correlation EffectsQuantum 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 mixtureObserving 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 Mechanical Stabilization of a Collapsing Bose-Bose MixtureConsistent 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 mixtureEffective single-mode model of a binary boson mixture in the quantum droplet regionQuantum Criticality of Liquid-Gas Transition in a Binary Bose MixtureQuantum filaments in dipolar Bose-Einstein condensatesPath-Integral Monte Carlo Study on a Droplet of a Dipolar Bose–Einstein Condensate Stabilized by Quantum FluctuationGround-state phase diagram of a dipolar condensate with quantum fluctuationsDroplet Crystal Ground States of a Dipolar Bose GasPhonon Stability of Quantum Droplets in Dipolar Bose Gases
[1] Kapitza P 1938 Nature 141 74
[2] Tisza L 1938 Nature 141 913
[3] Landau L 1941 Phys. Rev. 60 356
[4]Landau L D and Lifshitz E M 1959 Fluid Mechanics 2nd edn (Oxford: Butterworth–Heinemann)
[5] Coste C 1998 Eur. Phys. J. B 1 245
[6] Taylor E and Griffin A 2005 Phys. Rev. A 72 053630
[7] Pastukhov V 2017 Phys. Rev. A 95 023614
[8] Hoyos C, Kim B S, and Oz Y 2013 J. High Energy Phys. 2013(11) 145
[9]Schmitt A 2014 Introduction to Superfluidity—Field-Theoretical Approach and Applications. Lecture Notes in Physics Vol 888
[10]Svistunov B V, Babaev E S, and Prokof'ev N V 2015 Superfluid States of Matter (New York: CRC Press)
[11] Gross E P 1961 Il Nuovo Cimento 20 454
[12]Pitaevskii L P 1961 Sov. Phys.-JETP 13 451
[13] Stringari S 1996 Phys. Rev. Lett. 77 2360
[14]Pethick C J and Smith H 2002 Bose-Einstein Condensation in Dilute Gases (Cambridge: Cambridge University Press)
[15]Leggett A J 2006 Quantum Liquids Bose condensation and Cooper Pairing in Condensed-Matter Systems (Oxford: Oxford University Press)
[16]Pitaevskii L and Stringari S 2016 Bose-Einstein Condensation and Superfluidity (Oxford: Oxford University Press)
[17] Wu W C and Griffin A 1996 Phys. Rev. A 54 4204
[18] Marques G C and Bagnato V S 2000 Phys. Rev. A 61 053607
[19] Marques G C, Bagnato V S, and Spehler D 2001 Phys. Rev. A 63 43607
[20] Nikuni T, Zaremba E, and Griffin A 1999 Phys. Rev. Lett. 83 10
[21] Nikuni T and Griffin A 2001 Phys. Rev. A 63 033608
[22] Griffin A, Nikuni T, and Zaremba E 2009 Bose-Condensed Gases at Finite Temperatures (Cambridge: Cambridge University Press)
[23] Boudjemâa A 2013 Phys. Rev. A 88 23619
[24] Adhikari S K and Salasnich L 2008 Phys. Rev. A 78 043616
[25]Ueda M 2010 Fundamentals and New Frontiers of Bose-Einstein Condensation (Singerpore: World Scientific)
[26] Ilinski K N and Stepanenko A S 1998 arXiv:cond-mat/9803233 [cond-mat.stat-mech]
[27] Minguzzi A, Chiofalo M L, and Tosi M P 1997 Phys. Lett. A 236 237
[28] Abad M, Recati A, Stringari S, and Chevy F 2015 Eur. Phys. J. D 69 126
[29] Andreev P A and Kuzmenkov L S 2008 Phys. Rev. A 78 053624
[30] Andreev P A 2021 Laser Phys. Lett. 18 055501
[31] Andreev P A, Mosaki I, and Trukhanova M I 2021 Phys. Fluids 33 067108
[32]Popov V N 1987 Functional Integrals and Collective Excitations (Cambridge: Cambridge University Press)
[33] Al K U, Andersen J O, Proukakis N P, and Stoof H T C 2002 Phys. Rev. A 66 013615
[34] Leggett A J 1970 Phys. Rev. Lett. 25 1543
[35] Leggett A J 1998 J. Stat. Phys. 93 927
[36] Chauveau G, Maury C, Rabec F, Heintze C, Brochier G, Nascimbene S, Dalibard J, Beugnon J, Roccuzzo S, and Stringari S 2023 arXiv:2302.01776 [cond-mat.quant-gas]
[37] Bohm D 1952 Phys. Rev. 85 166
[38]Andreev A and Bashkin E 1975 Sov. J. Exp. Theor. Phys. 42 164
[39]Khalatnikov I 1957 Sov. Phys.-JETP 5 542
[40] Fil D and Shevchenko S 2005 Phys. Rev. A 72 013616
[41] Nespolo J, Astrakharchik G E, and Recati A 2017 New J. Phys. 19 125005
[42] Romito D, Lobo C, and Recati A 2021 Phys. Rev. Res. 3 023196
[43]Fetter A L and Walecka J D 2012 Quantum Theory of Many-Particle Systems (Courier Corporation)
[44] Kohn W and Sham L J 1965 Phys. Rev. 140 A1133
[45] Cabrera C R, Tanzi L, Sanz J, Naylor B, Thomas P, Cheiney P, and Tarruell L 2018 Science 359 301
[46] Cheiney P, Cabrera C R, Sanz J, Naylor B, Tanzi L, and Tarruell L 2018 Phys. Rev. Lett. 120 135301
[47] 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
[48] D'Errico C, Burchianti A, Prevedelli M, Salasnich L, Ancilotto F, Modugno M, Minardi F, and Fort C 2019 Phys. Rev. Res. 1 033155
[49] Kadau H, Schmitt M, Wenzel M, Wink C, Maier T, Ferrier-Barbut I, and Pfau T 2016 Nature 530 194
[50] Ferrier-Barbut I, Kadau H, Schmitt M, Wenzel M, and Pfau T 2016 Phys. Rev. Lett. 116 215301
[51] Ferrier-Barbut I, Schmitt M, Wenzel M, Kadau H, and Pfau T 2016 J. Phys. B 49 214004
[52] Schmitt M, Wenzel M, Böttcher F, Ferrier-Barbut I, and Pfau T 2016 Nature 539 259
[53] Wenzel M, Böttcher F, Tim L, Ferrier-Barbut I, and Pfau T 2017 Phys. Rev. A 96 053630
[54] Chomaz L, Baier S, Petter D, Mark M J, Wächtler F, Santos L, and Ferlaino F 2016 Phys. Rev. X 6 041039
[55] Petrov D S 2015 Phys. Rev. Lett. 115 155302
[56] Hu H and Liu X J 2020 Phys. Rev. Lett. 125 195302
[57] Wang Y Q, Guo L F, Yi S, and Shi T 2020 Phys. Rev. Res. 2 043074
[58] Gu Q and Yin L 2020 Phys. Rev. B 102 220503
[59] Xiong Y C and Yin L 2022 Phys. Rev. A 105 053305
[60] He L, Li H, Yi W, and Yu Z Q 2022 arXiv:2209.13559 [cond-mat.quant-gas]
[61] Wächtler F and Santos L 2016 Phys. Rev. A 93 061603
[62] Saito H 2016 J. Phys. Soc. Jpn. 85 053001
[63] Bisset R N, Wilson R M, Baillie D, and Blakie P B 2016 Phys. Rev. A 94 033619
[64] Baillie D and Blakie P B 2018 Phys. Rev. Lett. 121 195301
[65] Zhang F and Yin L 2022 Chin. Phys. Lett. 39 060301