## Abstract

We study collective modes of superfluid Bose gases in optical lattices at commensurate fillings. We focus on the vicinity of the quantum phase transition to the Mott insulator, where there exists the Higgs amplitude mode in addition to the Nambu–Goldstone phase mode associated with the spontaneous U(1) symmetry breaking. We analyze finite-temperature effects on the damping of the collective modes by using an effective spin-1 model and field-theoretical methods based on the finite-temperature Green’s function. We calculate the damping rates up to 1-loop order and evaluate them analytically and numerically. We show that the damping rate of the Higgs mode increases with increasing temperature but it remains underdamped up to the typical temperature achieved in experiments. Moreover, we find that the Nambu–Goldstone mode attenuates via a Landau damping process resulting from interactions with the Higgs mode and it can be overdamped at the typical temperature in a certain parameter region.

## Introduction

In a system with spontaneous breaking of continuous symmetry, there exists a gapless
collective mode that corresponds to a motion restoring the symmetry, namely the
Nambu–Goldstone (NG) mode [1,2]. If the system also has particle–hole
symmetry, a gapful collective mode emerges. It corresponds to fluctuations of the
order-parameter amplitude and is often referred to as the Higgs mode because of its
analogy with the Higgs scalar boson in particle physics [3]. In recent years, Higgs modes in condensed matter and
ultracold gases have attracted particular attention [4,5] thanks to
experimental developments for detecting these modes in various systems, including
superconductors $${\rm
NbSe_{2}}$$ [6–10] and Nb$$_{1-x}$$Ti$$_x$$N [11–13]^{1}, quantum antiferromagnets $$\rm TlCuCl_3$$ [16,17] and $$\rm KCuCl_3$$ [18], charge-density-wave materials $${\rm
K_{0.3}MoO_3}$$ [19,20] and $${\rm
TbTe_3}$$ [21,22], superfluid $$\rm
^3He$$ B-phase [23,24], and
superfluid Bose gases in optical lattices [25,26].

As for bosons in optical lattices, it has been predicted [27,28] that the Higgs mode emerges in the superfluid state near the quantum phase transition to the Mott insulator [29] at commensurate fillings, where the system has an approximate particle–hole symmetry [30]. The most careful experimental analyses regarding the Higgs modes thus far have been made in Ref. [26]; the response of 2D gases to temporal modulation of the lattice amplitude has been measured as a function of the modulation frequency by means of optical-lattice microscope techniques. Although the measured onset frequency of the response agrees with the theoretically computed energy gap of the Higgs mode, the response exhibits a broad continuum above the onset rather than a sharp peak. This means that the existence of the Higgs mode as a well-defined resonance has not yet been experimentally verified in this system.

The experiments of Ref. [26] have triggered extensive theoretical investigations on the Higgs modes of 2D bosons in optical lattices [31–39]. In particular, previous studies using the quantum Monte Carlo simulations of the 2D Bose–Hubbard model [31,37] have shown that the broadening of the spectral response is due to the combined effects of quantum fluctuations, finite temperatures, and spatial inhomogeneity caused by the confinement potential. Because the effects of quantum and thermal fluctuations are in general weaker in higher dimensions, it is expected that 3D systems are advantageous over 2D for observing a resonance peak of the Higgs mode in experiments. Altman and Auerbach [27] have indeed shown that the resonance peak of the Higgs mode in a homogeneous Bose–Hubbard system at zero temperature is significantly sharper in 3D than in 2D, i.e., the damping of the Higgs mode is weaker in 3D. As a next step, it is important to address the effects of finite temperatures in 3D.

It is worth noting that the optical-lattice microscope techniques, which have been used for measuring the Higgs gap in 2D [26], can be applied for analyzing some properties of 3D systems. In typical experiments using the optical-lattice microscope, a single 2D gas is extracted from a multilayer gas in 3D by selectively imposing an RF pulse on the single 2D layer and blasting out the other layers before actual manipulations of the system [40]. In order to analyze Higgs modes in 3D systems, one needs to perform the extraction protocol after preparing a 3D superfluid Bose gas in an optical lattice and performing a temporal modulation of the lattice amplitude. Indeed, coherent dynamics of bilayer systems (thin 3D systems) has been observed in such a manner [41].

Another interesting aspect of the Higgs mode is that its interaction with the NG mode may dramatically change properties of the latter mode. For instance, Nakayama et al. [42] recently predicted the Fano resonance of the NG mode mediated by a bound Higgs mode localized around potential barriers as a result of coupling between the two modes induced by the barriers. At finite temperatures, dynamically excited NG modes can interact with thermally excited Higgs modes and one of the natural consequences of such interactions should be the damping of the NG mode.

In this paper, we study damping of the Higgs and NG modes of 3D superfluid Bose gases in homogeneous optical lattices in the vicinity of the Mott transition at commensurate fillings, with a particular focus on the effects of finite temperatures. To calculate the damping rates, we use a low-energy effective model [27], which has the same form as the $$S=1$$ XY model with uniaxial single-ion anisotropy, and the field-theoretic approach based on the finite-temperature Green’s function. We analytically obtain approximate expressions of the damping rates for the two modes and also present numerical evaluations of the damping rates to clarify the validity region of the analytical formulae. We show that the damping rate of the Higgs mode at zero momentum increases with increasing temperature but it remains smaller than the oscillation frequency of the mode at the finite temperatures that can be realized in typical experiments. Moreover, we find that the interactions between the two modes allow for the Landau damping process of the NG mode, in which the damping rate can be even larger than the mode frequency at the typical temperatures, i.e., the NG mode can be overdamped.

This paper is organized as follows. In Sect. 2, we briefly review the mapping of the Bose–Hubbard model onto the effective spin-1 model. In Sect. 3, we explain our methods to calculate the damping rates at finite temperatures. In Sect. 4, using the formula obtained in Sect. 3, we calculate the damping rate of the Higgs mode to discuss its dependence on the temperature and the interaction strength. In Sect. 5, we discuss the properties of the damping of the NG mode. In Sect. 6, the results are summarized. Throughout the discussion, we set $$\hbar=k_{\rm B}=a=1$$, where $$\hbar$$, $$k_{\rm B}$$, and $$a$$ denote the reduced Plank constant, the Boltzmann constant, and the lattice constant.

## Model

We consider ultracold bosonic atoms in a cubic optical lattice. Assuming that the lattice is sufficiently deep, the system can be described by the Bose–Hubbard model within the tight-binding approximation [43,44]:

In the vicinity of the $${\bar n}$$th Mott phase, only three states {$$|{\bar n}-1\rangle_i$$, $$|{\bar n}\rangle_i$$, $$|{\bar n}+1\rangle_i$$} per site dominate the low-energy behavior of the superfluid near the Mott phase because the local fluctuations from the mean-field ground state in the Mott phase $$\prod_i|{\bar n}\rangle_i$$ are sufficiently suppressed so that we can ignore the high-energy excited states $$\{|{\bar n}\pm2\rangle_i, |{\bar n}\pm3\rangle_i, \ldots\}$$ in the complete Hilbert space. In the reduced Hilbert subspace, the operators in Eq. (1) take an approximated form. In order to obtain the specific expressions, it is convenient to introduce Schwinger bosons [45] in the following manner:

Substituting these expressions into the Bose–Hubbard model (1), we obtain an effective model that describes the low-energy behavior in the reduced space and has the same form as the XY model with uniaxial single-ion anisotropy and magnetic coupling [27]:

## Methods

*Canonical transformation*

Let us explain how to describe the collective excitations of the strongly correlated superfluid by using the effective model (4), within a mean-field approximation developed by Altman and Auerbach [27]. In the following discussion, we deal with the particle–hole symmetric case, i.e., $$h=\delta \mu =0$$. First, we define a canonical transformation as

### Holstein–Primakoff expansion

We assume that all of the fluctuations from the mean-field ground state are small. Then, the effective model (4) represented by the new basis $$b_{mi}$$$$(b_{mi}^{\dagger})$$$$(m=0,1,2)$$ can be simplified by employing a Holstein–Primakoff expansion [46]. With the local constraint $$\sum_{m=0}^{m=2}b^{\dagger}_{mi}b_{mi}=1$$, this expansion eliminates $$b_{0i}$$$$(b^{\dagger}_{0i})$$ from the model (4) as

Substituting the expansion (7) into the model (4), we obtain the simplified model, which has sequent terms with each order, and would show them up to the third order:

*Bogoliubov transformation*

To diagonalize $$H_{\rm eff}^{(2)}$$, let us perform a Bogoliubov transformation [51] in each branch defined as

The operator $$\beta^{\dagger}_{m{\bf k}}$$ creates a quasiparticle with the excitation energy $${\cal E}_{m\bf k}$$, which is given by

It is easily seen that the dispersion of the Higgs mode $${\cal E}_{1\bf k}$$ has an energy gap $$\Delta \equiv 2J{\bar n}z\sqrt{1-u^2}$$ at zero momentum, while that of the NG mode $${\cal E}_{2\bf k}$$ is gapless. Obviously, the gap of the Higgs mode closes at the critical point $$u=u_c(=1)$$.

Equation (16) has third-order terms characterizing the interactions among the three excitations. Specifically, $$H_{\rm eff}^{(3)}$$ has five types of interaction term, which are shown in Fig. 1, and their Hermite conjugates. The interaction terms shown in Figs. 1(a) and (b) generate scattering processes closed only in the Higgs branch, while the terms shown in Figs. 1(c), (d), and (e) couple the two different branches with one Higgs mode and two NG modes. As we will see in the following subsection, these terms cause damping of the elementary excitations. It is worth noting that $$H_{\rm eff}^{(3)}$$ does not include vertices consisting only of propagators of the NG mode. This is a generic property of a superfluid with particle–hole symmetry [52] and is in stark contrast to the case of a weakly interacting Bose gas, which has such vertices [53–55].

When $$|{\bf k}| \ll 1$$, the above expressions are approximated with simpler forms. The simplification is done by the Taylor expansion of $$\gamma_{\bf k}$$ with respect to $$\bf k$$. For example, the excitation energies (17) and (18) become

These expressions are isotropic in the momentum space, thus allowing for analytical evaluations of the damping rates.

*Finite-temperature Green’s function*

To calculate the damping rates of the Higgs and NG modes, we use the field-theoretical approaches based on the finite-temperature Green’s function (see Refs. [56–58]). In our calculation, we define the Green’s functions of the Higgs ($$m=1$$) and NG ($$m=2$$) modes as

The third-order action $${\cal S}^{(3)}_{\rm eff}$$ can be obtained from $$H_{\rm eff}^{(3)}$$ with mere replacement of the operators with the field variables.

In field theory, the damping rates can be calculated as imaginary parts of the self-energies of the Higgs mode $$\Sigma_1(i\omega_n;{\bf k})$$ and of the NG mode $$\Sigma_2(i\omega_n;{\bf k})$$. We define them as $$\Gamma_{m{\bf k}}\equiv {\rm Im}\Sigma_{m}({\cal E}_{m\bf k}+i\epsilon;{\bf k})$$ where $$\epsilon$$ is a positive-valued infinitesimal quantity. Calculating these damping rates by perturbative expansion with respect to $${\cal S}_{{\rm eff}}^{(3)}$$ up to the second order, we obtain

The contributions to each damping rate consist of only one type of perturbative correction. Figure 2 shows the Feynman diagrams that provide nonzero contributions to the damping rates (25) and (26). In general, there are other diagrams with different structures from the ones shown in Fig. 2. However, we find that the actual contributions to the damping come from only the diagrams depicted in Fig. 2 within the 1-loop order, because of the energy–momentum conservations laws, which are represented as a $$\delta$$ function in Eqs. (25) and (26).

The diagram shown in Fig. 2(a) means that the Higgs mode attenuates by decaying into the two NG modes due to quantum and thermal fluctuations. This type of damping of collective modes is called Beliaev damping [53–59]. In general, Beliaev damping can occur even at zero temperature. The Beliaev damping of the Higgs mode at zero temperature was first predicted by Altman and Auerbach [27]. On the other hand, the diagram depicted in Fig. 2(b) means that the NG mode attenuates through processes in which the initial NG mode absorbs another thermally excited NG mode, and subsequently turns into the Higgs mode as the final state. This type of damping is called Landau damping [60–61]. At zero temperature, this damping cannot occur because there is no thermal excitation. Therefore, in the vicinity of the Mott phase, the NG mode does not attenuate at zero temperature within the approximations discussed above.

## Damping rate of the Higgs mode

In this section, we discuss the damping of the Higgs mode at finite temperatures by evaluating expression (25). In particular, we consider the case of the Higgs mode with zero momentum because a typical perturbation used for exciting the Higgs mode in experiments is the lattice-amplitude modulation with zero momentum [26].

First, let us evaluate the integrations of the formula (25) within a long-wavelength approximation, where $${\cal E}_{1{\bf k}}$$, $${\cal E}_{2{\bf k}}$$, $$u_{m{\bf k}}$$, and $$v_{m{\bf k}}$$ are approximated with Eqs. (19), (20), (21), and (22). This approximation is better justified in a closer vicinity of the critical point, $$u=u_c$$, where the energy of the NG mode dominant to the damping of the Higgs mode, $$\Delta/2$$, is smaller. Integrating with respect to $${\bf k}_1$$ and $${\bf k}_2$$ on the r.h.s. of Eq. (25) within the approximation, we obtain the simple formula

It is obvious from Eq. (28) that the dependence on temperature $$T$$ is determined by the factor $${\rm coth}(\beta\Delta/4)$$. At $$T=0$$ and $$|u-u_c|\ll 1$$ (recall that $$u_c=1$$), this factor becomes $$1$$ and the damping rate behaves as $$\Gamma_{1{\bf k}={\bf 0}}\sim J{\bar n}z\sqrt{1-u^2}$$. This consequence agrees with the previous one obtained by Altman and Auerbach [27]. In this sense, the analytic expression (28) gives the finite-temperature correction to the previous result at zero temperature.

In Fig. 3, we show the dependence on $$u$$ of the damping rate (28) at fixed temperatures and a comparison with the numerical calculation without the long-wavelength approximation. As expected, the analytic results coincide well with the numerical data in the vicinity of $$u=u_c$$, while they deviate as $$u$$ decreases from $$u=u_c$$. It is seen that, for a given $$u$$, the damping rate monotonically increases as the temperature increases. This tendency is stronger in closer vicinity of the critical point $$u=u_c$$.

The dimensionless quantity $$\Gamma_{1{\bf k}={\bf 0}}/\Delta$$ characterizes the behavior of the temporal oscillation of the Higgs mode and determines the width of the resonance peak in the spectral function. If $$\Gamma_{1{\bf k}={\bf 0}}/\Delta > 1$$, then the oscillation abruptly attenuates during one period, i.e., it is overdamped. From Fig. 4, where $$\Gamma_{1{\bf k}={\bf 0}}/\Delta$$ is plotted against $$u$$, we see that $$\Gamma_{1{\bf k}={\bf 0}}/\Delta$$ increases as the critical point is approached but the Higgs mode remains underdamped ($$\Gamma_{1{\bf k}={\bf 0}}/\Delta<1$$) even at $$T=0.3 J\bar{n}z$$. Provided that $$\bar{n}\gg1$$ is assumed and that the temperature can be as low as $$T/J=O(1)$$ in typical experiments with bosons in optical lattices [26,62], our result implies that the Higgs mode has a sharp resonance peak in the spectral function at the typical temperatures, at least in the absence of a trapping potential. Notice that in Figs. 3 and 4 we have not shown the data points in the parameter region, where $$T>\Delta \sim T_c$$, because our method is valid only in the superfluid phase. Here, $$T_c$$ denotes the transition temperature from a superfluid to a normal fluid.

Our calculations do not take into account the logarithmic correction to the damping rate, which stems from the renormalization of an effective coupling constant because $$d=3$$ is the upper critical dimension [63]. Due to the correction, $$\Gamma_{1{\bf k}={\bf 0}}/\Delta$$ at $$T=0$$ approaches zero as $$\sim 1/\ln|u-u_c|$$ in the limit that $$u\rightarrow u_c$$. Ignorance of the logarithmic correction is not problematic in the practical sense that the parameter region, where the correction is effective, is so narrow that it is very difficult to observe in either experiments or numerical simulations, especially for dynamical quantities like damping rates.

## Damping rate of the NG mode

In this section, we evaluate the damping rate of the NG mode expressed in Eq. (26) in specific parameter regions of interest. In cold-atom experiments, the NG mode can be dynamically excited by means of two-photon Bragg scattering techniques [28,64–66], which allow for wide control of the momentum $${\bf k}$$. Hence, we analyze the $${\bf k}$$-dependence of the damping rate in addition to the $$u$$-dependence and the $$T$$-dependence. First, let us obtain analytical expressions of the damping rate within the long-wavelength approximation. Substituting Eqs. (19), (20), (21), and (22) into Eq. (26) and integrating it with respect to $${\bf k}_1$$ and $${\bf k}_2$$ except for the variable $$|{\bf k}_2|$$ lead to

These parameters determine the maximum and minimum momenta of the thermally excited NG modes that are absorbed by the initial NG mode in the Landau damping process. Carrying out the remaining integration with respect to $$|{\bf k}_2|$$, we obtain

We emphasize that this approximate expression (31) is well justified when $$k_l \ll 1$$. The condition $$k_l \ll 1$$ can be converted to

When the conditions that $$k_l \ll k_u$$ and $$c_{\rm ng}k \ll T \ll c_{\rm ng}k_u$$ are satisfied, a further simplification of Eq. (31) can be made as follows:

Equation (33) clearly shows that the damping rate is proportional to the Bose distribution function, $$f_{\rm B}( c_{\rm ng}k_l )$$, of the thermally excited NG modes with the lowest momentum $$k_l$$ that is allowed by the conservation law.

In Fig. 5, we show the damping rate of the NG mode $$\Gamma_{2{\bf k}}/(J\bar{n}z)$$ as a function of $$u$$ for several values of the temperature $$T$$ and the initial momentum $${\bf k}$$. In addition to the results of the analytical expression (31), we plot the data points obtained from numerical integration of Eq. (26) with respect to $${\bf k}_1$$ and $${\bf k}_2$$ for comparison. In order to clarify the effects of the anisotropy of $$\gamma_{\bf k}$$ in the momentum space, which are ignored in Eq. (26), we show the cases of three different directions, namely $${\bf k}/k=(1,0,0)$$, $$(\sqrt{3}/2,1/2,0)$$, and $$(1/\sqrt{3},1/\sqrt{3},1/\sqrt{3})$$.

From Eq. (31) and Fig. 5, we see the following four generic tendencies. First, for given $$u$$ and $${\bf k}$$, the damping rate increases with increasing temperature. This is very natural in the sense that the number of thermally excited NG modes, which are a main source of Landau damping, is larger for a higher temperature. Second, the effects of the anisotropy are more noticeable for larger $$k$$, which is also natural. Third, the analytical results of Eq. (31) better agree with the numerical data at a closer vicinity of the critical point, as expected from the validity condition (32). Fourth, it is most remarkable that the damping rate is significantly large near the critical point and monotonically decays into zero from its peak as $$|u-u_c|$$ increases.

Let us explain the fourth tendency from the viewpoint of the energy–momentum conservation law. If the Higgs gap satisfies the condition that $$c_{\rm ng}k\simeq \Delta/2$$ near $$u=u_c$$, as illustrated in Fig. 6(a), the lowest energy of the thermal NG mode that is absorbed by the initial NG mode is given by $$c_{\rm ng}k_l \simeq \Delta/2$$ as a result of the energy–momentum conservation law. This means that the Bose distribution function in Eq. (33) can be $$f_{\rm B}( c_{\rm ng}k_l )=O(1)$$ even at $$T<\Delta \sim T_c$$. In other words, in this case there are sufficiently many thermal NG modes that the initial NG mode can absorb to attenuate. On the other hand, if $$|u-u_c|$$ increases for fixed $$T$$ such that $$c_{\rm ng}k\ll \Delta$$ (see Fig. 6(b)), the condition that $$c_{\rm ng}k_l > \Delta$$ is imposed by the energy–momentum conservation law. This implies that $$c_{\rm ng}k_l \gg T$$, thus leading to the exponential suppression of the damping rate.

In Fig. 7, we show the ratio $$\Gamma_{2{\bf k}}/(c_{\rm ng}k)$$ to characterize the behavior of the temporal oscillation of the NG mode. As seen in Fig. 7(a), when $$k=0.1$$ and $$T=0.3 J\bar{n}z$$, $$\Gamma_{2{\bf k}}/(c_{\rm ng}k)$$ exceeds unity near the critical point, i.e., the NG mode is overdamped. Figure 7(b) shows that the damping rate can be even larger by optimizing the wavelength of the initial NG mode $$k$$. One also sees from Fig. 7(b) that the analytical expression (31) completely fails in the limit of $$k\rightarrow 0$$, which is consistent with the validity condition (32).

We finally note the limitation of the 1-loop approximation regarding the predictability of overdamping of a collective mode. Because this approximation is a perturbative approach, $$\Gamma_{m{\bf k}}> \mathcal{E}_{m{\bf k}}$$, which is the condition of the overdamping, means that the perturbative correction is larger than the nonperturbative value. Therefore, one cannot judge whether or not the overdamping is an artifact of the 1-loop approximation until the higher-order corrections are evaluated.

## Summary

In this paper, we studied the damping of the Higgs and NG modes of Bose gases in a cubic optical lattice at finite temperatures. We calculated the damping rates by using the effective spin-1 model and field-theoretical methods. We derived analytic expressions of the damping rates within the long-wavelength approximation and confirmed their validity in the vicinity of the critical point through comparison with numerical calculations. We showed that, while the Higgs mode attenuates more significantly at higher temperatures, it is not overdamped at temperatures that can be achieved in typical experiments. This result indicates the feasibility of detecting the Higgs mode at 3D and finite temperatures as a resonance peak in a spectral function, at least when there is no trapping potential. As for the NG mode, we found parameter regions where the Landau damping process leads to overdamping of the NG mode, and discussed the origin of the strong damping, especially near the critical point.

In future studies addressing the detectability of the Higgs mode in 3D more quantitatively, it will be important to include the effects of the trapping potential and the breaking of the particle–hole symmetry. While a prescription to treat these effects at tree level has been presented in Ref. [28], one needs to extend it to the 1-loop level to take into account the effects of quantum and thermal fluctuations as well.

## Acknowledgements

The authors thank S. Tsuchiya for useful discussions. The authors also thank the Yukawa Institute for Theoretical Physics (YITP) at Kyoto University, where this work was initiated during the YITP workshop (YITP-W-14-02) on “Higgs Modes in Condensed Matter and Quantum Gases”. I.D. acknowledges Grants-in-Aid for Scientific Research from JSPS: Grant Nos. 25800228 and 25220711.

## Gutzwiller’s mean-field ansatz in the reduced Hilbert subspace

In this appendix, we briefly review the bosonic Gutzwiller’s mean-field ansatz in the reduced Hilbert subspace [27,28], which is used to construct the ground-state wave function near the critical point within the mean-field approximation, in order to supplement the canonical transformation (6) in Sect. 3.1. The main concern below is to understand the origin of the particular choice of the angle $$\theta_0={\rm tan}^{-1}(\sqrt{1-u^2}/u)$$.

In the bosonic Gutzwiller’s mean-field ansatz (in either complete or reduced Hilbert space), we begin with assuming a direct product state with complex-valued variational parameters $$(f_0,f_1,f_2,\ldots)$$ such as $$\prod_i \{ f_0 |0\rangle_i + f_1|1\rangle_i +f_2|2\rangle_i+\cdots \}$$ as the wave function of the system [5,27,28]. The state vector $$|n\rangle_i$$ is given by $$|n\rangle_i=(a^{\dagger}_i)^n/\sqrt{n!}|0\rangle$$ for $$n=0,1,2,\ldots$$ ($$|0\rangle$$ is the vacuum of bosons). Assuming that, in the vicinity of the critical point concerned with the $$\bar n$$th Mott phase, the dominant states are reduced into only three states ($$|{\bar n}-1\rangle$$, $$|{\bar n}\rangle$$, $$|{\bar n}+1\rangle$$), we can define the variational wave function as

Let us minimize the energy (A2) with respect to the variational parameters $$(\theta,\eta,\varphi,\chi)$$ to obtain the ground-state wave function and construct the canonical transformation (6). In the following discussion, we restrict ourselves to the commensurate large-filling case, i.e., $${\bar n}$$ is a large integer. Assuming that $$\delta \mu =0$$ and effecting the variation of the energy (A2) with respect to $$\eta$$, $$\varphi$$, $$\chi$$, and $$\theta$$, we find that $$\eta=0$$, $$\varphi=\varphi_0$$, which is any constant, $$\chi=\pi/2$$, and $${\rm sin}^2(\theta/2)=(1-u)/2 \equiv {\rm sin}^2(\theta_0/2)$$ and that the phase transition between the superfluid and Mott insulator phases occurs at $$u=1$$. Substituting these results into Eq. (A1), we obtain the ground-state wave function

Notice that, in the current discussion, $$\varphi_0$$ explicitly remains as a constant for convenience while we set it to be zero in Sect. 3. Similarly, we also obtain the order parameter characterizing the superfluid order on the ground state $$\Psi^* \approx \sqrt{{\bar n}/2}{\rm sin}\theta_0 e^{-i\varphi_0}$$. This expression reveals that the variational parameters $${\theta_0}$$ and $$\varphi_0$$ correspond to the amplitude and phase degrees of freedom of the order parameter.

Next, let us consider the amplitude and phase fluctuations on the ground state. Creating these modes corresponds to replacements of the parameters ($$\theta_0$$, $$\varphi_0$$) such as $$\theta_0 \rightarrow \theta_0+\delta \theta_i$$ and $$\varphi \rightarrow \varphi_0+\delta \varphi_i$$, where the local deviations $$\delta \theta_i$$ and $$\delta \varphi_i$$ are small quantities. Then, we find that the ground-state wave function (A4) becomes a direct product state of a linear combination spanned by