Is density fluctuation gapless in superfluid? Deep in the superfluid phase, the superfluid order parameter $\phi$ can be decomposed into the amplitude (density) mode $\rho$ and the phase mode $\theta$ as
$$\phi=\sqrt{\rho} e^{\mathrm{i}\theta}.$$
It is believed that the density fluctuations are gapped as they correspond to "climbing up" the potential $V(\phi)$. What remains gapless at low energy are the phase fluctuations, or the Goldstone modes, described by
$$\mathcal{L}[\theta]=\frac{1}{2g}\big((\partial_{t}\theta)^2-(\partial_{\boldsymbol{x}}\theta)^2\big).$$
However, there is an emergent density operator $\rho=-\partial_t\theta$ and current operator $\boldsymbol{j}=\partial_{\boldsymbol{x}}\theta$ in terms of the phase field $\theta$, such that the continuity equation $\partial_{t}\rho+\partial_{\boldsymbol{x}}\boldsymbol{j}=0$ is satisfied on-shell. This means that the density fluctuation is actually gapless, as seen from the density-density correlation function in the momentum-frequency space:
$$\langle\rho_{k}\rho_{-k}\rangle=-\omega^2\langle\theta_{k}\theta_{-k}\rangle=\frac{\omega^2}{\omega^2-\boldsymbol{k}^2}.$$
If the density fluctuation in the superfluid is indeed gapless, how can we ignore them and claim $\mathcal{L}[\theta]$ alone as the effective description of the superfluid dynamics at low energy? But on the other hand, the classical picture of "claiming up" the Maxican-hat potential $V(\phi)$ does imply that the density fluctuation should be gapped. How to reconcile this contradiction?
 A: The density mode in a superfluid is indeed gapless (in fact, the density mode in a normal phase is gapless, too. This mode is called sound). 
The confusion arises because in order to arrive at the effective lagrangian for $\theta$, you have to integrate out the amplitude mode. As a result, the $\theta$ parameter in the effective lagrangian couples to the density. 
Postscript: The amplitude is not directly related to the superfluid density $\rho_s$. The superfluid density is defined by 
$$
\vec\pi = \rho_s v_s + \rho_n v_n
$$
where $\vec\pi$ is the momentum density and $v_s= i\hbar\nabla\theta/m$ is the superfluid velocity. This means that $\rho_s$ governs the response in momentum to gradients of the phase. Experimentally, $\rho_s$ is extracted by measuring the velocities of first and second sound.
A: Superfluids are Galilean invariant therefore it is a good idea to start from a Galilean invariant model when trying to understand their  dynamics.  For example the  Gross-Pitaevski (GP) model that comes from the  action integral 
$$
S[\phi, \phi^\dagger]= \int d^3x dt\left\{ \phi^\dagger (i \partial_t + \frac 1{2m} \nabla^2) \phi + \mu\phi^\dagger \phi -\frac 12 \lambda (\phi^\dagger\phi)^2\right\}.
$$
is Galilean invariant.
This action contains a Mexican hat potential
$$
V(\phi)= \frac 12 \lambda  (\phi^\dagger\phi)^2-\mu \phi^\dagger\phi
$$
which is minimized at $\phi^\dagger\phi = \mu/\lambda$.
The possible stationary solutions are therefore 
$$
\langle\phi \rangle= \phi_c= e^{i\theta} \sqrt{\frac \mu \lambda}
$$ 
In the GP model the $\rho$ in $\phi= \sqrt{\rho} e^{i\theta}$ really is the particle density because this model applies at very low temperature when essentially every particle is in the condensate. So the equilibrium particle density is $\rho= \mu/\lambda$.
If we look for small oscillations $\phi=\phi_c+\eta$ then
$$
V(\phi+\eta) \sim const+ \mu \eta^\dagger \eta +\frac 12\mu (\eta^2+(\eta^\dagger)^2)+O(\eta^3)
$$
and the linearized  equations of motion become
$$
i\partial_t \eta = -\frac1{2m} \nabla^2 \eta +\mu\eta +\mu \eta^\dagger,\\
-i\partial_t \eta = -\frac1{2m} \nabla^2 \eta^\dagger +\mu\eta +\mu \eta^\dagger.
$$
If we seek a solution 
$$
\eta= a e^{ikx-i\omega t} +b^\dagger e^{-ikx+i\omega t}
$$
we find that $(a,b)^T$ must obey
$$
\left[\begin{matrix} k^2/2m =\omega +\mu & \mu\cr \mu & k^2/2m +\omega+\mu\end{matrix}\right] \left[\begin{matrix} a \cr b\end{matrix}\right]=0,
$$
so the allowed frequencies are given by
$$
\omega^2 =(k^2/2m +\mu)^2 -\mu^2.
$$
At small $k$ this becomes becomes $\omega^2=c^2k^2$ with $ c^2 =\lambda \rho_0 /m$. These modes are  the gapless sound waves. During the motion the tip of the $\phi$ vector describes an ellipse about the equilibrium $\phi_c$. These sound modes are therefore a combination of Goldstone-like motion along the bottom of the Mexican hat potential well and an out of phase  radial ``Higgs- like'' radial  oscillation. There are  no separate circumferential "Goldstone" and radial "Higgs" modes in the non-relativistic bose-condensed superfluid. The coupled modes are sound waves with a   density fluctuation $\rho_0\to \rho_0+\delta \rho$ and a simultaneous (in-phase) back-and-forth velocity given by $v=\nabla\theta$.    
We had two equations for   $\eta$ that were first order in time. We can, if we like,  eliminate $\rho$ to get a second order wave  equation involving only $\theta$, or eliminate $\theta$ to get a wave equation  involving only $\rho$ --- but we won't get  a second-order-in-time equation involving both variables.  Beware, however, the linearized wave equations are not Galilean invariant. Futher, focusing only on the equations of motion   risks discarding   the $i\rho_0 \partial_t\theta$ in the action integral on the grounds that  it is total derivative.  This topological winding-number term  is essential for vortex dynamics where  it is responsible for the Magnus effect. In  its absence (as asked recently on this site) a propeller would not work in a superfluid. 
