I am following chapter 3 of X. G. Wen's book "Quantum Field Theory of Many-Body Systems". The following action for a weakly interacting Bose gas is derived: $$S[\varphi,\varphi^*] = \int dt \ d^D \textbf{r} \bigg[\frac{1}{2}i(\varphi^*\partial_t \varphi -\varphi \partial_t \varphi^*)-\frac{1}{2m}\nabla \varphi^* \nabla \varphi+\mu |\varphi|^2 -\frac{V_0}{2}|\varphi|^4\bigg]$$
The equations of motion can easily be found to be \begin{align} &\bigg(-i\partial_t - \frac{1}{2m}\nabla^2-\mu + V_0|\varphi|^2\bigg)\varphi=0\\ &\bigg(i\partial_t - \frac{1}{2m}\nabla^2-\mu + V_0|\varphi|^2\bigg)\varphi^*=0 \end{align}
Now if $\mu>0$ we are in the superfluid phase so the ground state is $|\varphi_0| = \sqrt{\frac{\mu}{V_0}} \equiv \sqrt{\rho_0}$. To get the low energy excitations I consider small fluctuations from this ground state by setting $\varphi = \varphi_0 + \delta \varphi$. After linearising the equations of motion I get that \begin{align} &\bigg(-i\partial_t - \frac{1}{2m}\nabla^2+\mu\bigg)\delta \varphi=-\mu \delta \varphi^*\\ &\bigg(i\partial_t - \frac{1}{2m}\nabla^2+\mu \bigg)\delta \varphi^*=-\mu \delta \varphi \end{align} which can be combined to give \begin{equation} \bigg|-i\partial_t - \frac{1}{2m}\nabla^2+\mu\bigg|^2 \delta \varphi = \mu^2 \delta \varphi \end{equation} The resulting dispersion relation is \begin{equation} \omega = \frac{1}{2m}k^2 \ \text{ and } \ \omega = \bigg|\frac{1}{2m}k^2-2\mu\bigg| \end{equation}
However the second of these looks very odd to me.
Furthermore, I was expecting to get a gapless, linearly propagating Goldstone mode $\omega \propto k$ (due to the spontaneously broken U(1) symmetry). Indeed section 3.3.3 in Wen's book explains how to get these modes by introducing the fields $\theta$ and $\delta \rho$ (the heavy field) as \begin{equation} \varphi = \sqrt{\rho_0+\delta \rho}e^{i\theta} \end{equation} and integrating $\delta \rho$ out to get a low energy effective field theory for $\theta$.
So what excitations has my method produced instead? Have I found some sort of higher-energy excitations, perhaps due to density fluctuations? To me it seems like setting $\rho = \rho_0+\delta \rho$ corresponds to taking $\theta=0$ and $\delta \phi \sim \delta \rho$, so that my method found the dispersion for the heavy $\delta \rho$ fluctuations.