Damping of the Higgs and Nambu-Goldstone modes of superfluid Bose gases at finite temperatures

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 the 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 the temperature but it remains underdamped up to a 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.

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 PTEP 2016, 063I01 K. Nagao and I. Danshita 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][32][33][34][35][36][37][38][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 lowenergy 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 = k B = a = 1, where , k 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]: where a † i and a i denote the creation and annihilation operators for a boson at site i and satisfy the bosonic commutation relations i a i is the filling factor. The symbol i, j denotes the nearest-neighbor sites. The parameters J , U , and δμ are the hopping energy, the on-site interaction energy, and the chemical potential. Note that at a sufficiently large fillingn 1, δμ = 0 corresponds to the commensurate case, i.e., integer fillingn ∈ N. This definition simplifies the following derivation of the effective model to describe the system near the Mott phase with fillingn. In this paper, we focus on the high-filling case.
In the vicinity of thenth Mott phase, only three states {|n − 1 i , |n i , |n + 1 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 i |n i are sufficiently suppressed so that we can ignore the high-energy excited states {|n ± 2 i , |n ± 3 i , . . .} 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: where |vac is a vacuum of the Schwinger bosons. These operators satisfy the bosonic commutation relations [t mi , t † nj ] = δ m,n δ i,j (m, n = −1, 0, 1) and we impose a local constraint m=1 m=−1 t † mi t mi = 1 on them to eliminate nonphysical states, e.g., t † 1i t † 0i |vac . The operator t 0i represents the local mean field in the Mott phase and t 1i (t −1i ) is the single particle (hole) excitation at site i. The bosonic operator a † i can be expressed in terms of the Schwinger bosons (2): 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]: where h = δμ. The spin operators are written in terms of the Schwinger bosons (2): and obey the standard SU(2) algebra. FornJz/U 1, the ground state is the XY ferromagnetic ordered state where the U(1) symmetry is spontaneously broken. Here z = 2d is the coordination 3/16 Downloaded from https://academic.oup.com/ptep/article-abstract/2016/6/063I01/2605974 by guest on 28 July 2018 PTEP 2016, 063I01 K. Nagao and I. Danshita number and d = 3 is the spatial dimension of the system. FornJz/U 1, the ground state has U(1) symmetry and the mean-field wave function is given by i |S = 1, m z = 0 i ≡ i |n i . Notice that, in the h = δμ = 0 case, the effective model (4) has particle-hole symmetry, which corresponds to that emerging in the model (1) near the superfluid-Mott insulator transition at commensurate fillings. In fact, except the last term, the effective model (4) is invariant under the interchange of particle excitations t † 1i (t 1i ) and hole excitations t † −1i (t −1i ).

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 = δμ = 0. First, we define a canonical transformation as which generates a rotation of the old basis spanned by three Schwinger bosons (t 1i , t 0i , t −1i ) into a new basis. The new operators also satisfy the bosonic commutation relations is a dimensionless parameter (the reason for this choice will be explained in Appendix A). The transformation (6) can be derived by the bosonic Gutzwiller mean-field ansatz in the reduced space [27,28]. The operators b 0i and b 1i (b 2i ) stand for the mean-field ground state and the excitations in the Higgs (NG) branch. Note that, in the mean-field theory and the high-filling limitn 1, the superfluid to Mott insulator transition occurs at u = 1 ≡ u c . In other words, a quantity |u − u c | = |u − 1| measures the distance from the critical point u = 1.

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 ) (m = 0, 1, 2) can be simplified by employing a Holstein-Primakoff expansion [46]. With the local constraint m=2 where m = {1, 2}. The higher-order terms are ignored if we take into account the leading vertex terms with third order with respect to the excitations. This approximation is similar to the spin-wave expansion in localized spin systems with long-range orders (more details can be found in, e.g., Ref. [47]). Notice that, in the case of finite temperatures, the above truncation is valid only at d ≥ 3. At lower dimensions, the thermal fluctuations are so strong that they destroy the long-range order of the mean-field ground state [48][49][50] and one cannot justify ignorance of the higher terms in the expansion (7). 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: where the index l in H (l) eff means lth order with respect to the operators b † mi and b mi (m = 1, 2). Here, we perform the Fourier transformation of the operators defined as where N is the system size and the vector x i denotes site i. The summation with respect to momentum k is taken over the first Brillouin zone. Then, each term H (l) eff can be written as the following equations, respectively: where c ≡ cos(θ 0 /2) and s ≡ sin(θ 0 /2). The explicit expressions of some coefficients in Eqs. (11) and (12) are given by where γ k = e e −ik·e /z = d s=1 cos(k s )/d.

Bogoliubov transformation
To diagonalize H (2) eff , let us perform a Bogoliubov transformation [51] in each branch defined as where m = {1, 2}. The new operators obey the same commutation relations [β mp , β † nq ] = δ n,m δ p,q as the previous operator such that the coefficients u mk and v mk satisfy a relation |u mk | 2 − |v mk | 2 = 1. Assuming that u mk and v mk are real numbers and imposing the condition that H (2) eff in the new representation has no anomalous term, we determine the coefficients where sgn(x) is the sign function. After the Bogoliubov transformation, Eq. (4) reads The operator β † mk creates a quasiparticle with the excitation energy E mk , which is given by It is easily seen that the dispersion of the Higgs mode E 1k has an energy gap ≡ 2Jnz √ 1 − u 2 at zero momentum, while that of the NG mode E 2k 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 eff 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 eff 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][54][55].
When |k| 1, the above expressions are approximated with simpler forms. The simplification is done by the Taylor expansion of γ k with respect to k. For example, the excitation energies (17) and (18) become  where c h = 2Jn √ zu and c ng = Jn √ z (1 + u). Similarly, if we assume that u = u c , the coefficients (14) and (15) 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][57][58]). In our calculation, we define the Green's functions of the Higgs (m = 1) and NG (m = 2) modes as where ω n = 2πn/β (n = 0, ±1, ±2, . . .) is the Matsubara frequency, β = T −1 is the inverse temperature, β mk (iω n ) and the conjugateβ mk (iω n ) = (β mk (iω n )) * are complex-valued field variables at (ω n , k), and D(β,β) is a measure of the integrations. The effective action S eff = S (2) eff + S eff is derived from the effective model (16) and the quadratic action S (2) eff is given by The third-order action S eff can be obtained from H eff 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 1 (iω n ; k) and of the NG mode 2 (iω n ; k). We define them as mk ≡ Im m (E mk + i ; k) where is a positive-valued infinitesimal quantity. Calculating these damping rates by perturbative expansion with respect to S (3) eff up to the second order, we obtain  Fig. 2. Contribution to the self-energy of (a) the Higgs mode 1 (iω n ; k) and (b) the NG mode 2 (iω n ; k).
Other diagrams do not yield a nonzero contribution to the damping rates.
where δ(x) is the δ function, f B (x) = 1/(e βx − 1) is the Bose distribution function, and the matrix elements M k 1 ,k 2 ,k 3 are given by 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 δ 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 E 1k , E 2k , u mk , and v mk 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, /2, is smaller. Integrating with respect to k 1 8/16 PTEP 2016, 063I01 K. Nagao and I. Danshita It is obvious from Eq. (28) that the dependence on temperature T is determined by the factor coth(β /4). At T = 0 and |u − u c | 1 (recall that u c = 1), this factor becomes 1 and the damping rate behaves as 1k=0 ∼ Jnz √ 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 1k=0 / characterizes the behavior of the temporal oscillation of the Higgs mode and determines the width of the resonance peak in the spectral function. If 1k=0 / > 1, then the oscillation abruptly attenuates during one period, i.e., it is overdamped. From Fig. 4, where 1k=0 / is plotted against u, we see that 1k=0 / increases as the critical point is approached but the Higgs mode remains underdamped ( 1k=0 / < 1) even at T = 0.3Jnz. Provided thatn 1 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 > ∼ 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, 1k=0 / at T = 0 approaches zero as ∼ 1/ ln |u − u c | in the limit that u → 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 9  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][65][66], which allow for wide control of the momentum k. Hence, we analyze the k-dependence of the damping rate in addition to the udependence 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 k 1 and k 2 except for the variable |k 2 | lead to where k ≡ |k|. The upper and lower bounds of the integration k u and k l are given by 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 |k 2 |, we obtain 2k = We emphasize that this approximate expression (31) is well justified when k l 1. The condition k l 1 can be converted to  When the conditions that k l k u and c ng k T c 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 B (c 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 2k /(Jnz) as a function of u for several values of the temperature T and the initial momentum 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 k 1 and k 2 for comparison. In order to clarify the effects of the anisotropy of γ k in the momentum space, which are ignored in Eq. (26), we show the cases of three different directions, namely k/k = (1, 0, 0), ( √ 3/2, 1/2, 0), and (1/ √ 3, 1/ √ 3, 1/ √ 3). From Eq. (31) and Fig. 5, we see the following four generic tendencies. First, for given u and 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 ng k /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 ng k l /2 as a result of the energy-momentum conservation law. This means that the Bose distribution function in Eq. (33) can be f B (c ng k l ) = O(1) even at T < ∼ 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 ng k (see Fig. 6(b)), the condition that c ng k l > is imposed by the energy-momentum conservation law. This implies that c ng k l T , thus leading to the exponential suppression of the damping rate. In Fig. 7, we show the ratio 2k /(c 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.3Jnz, 2k /(c 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 → 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, mk > E mk , which is the condition of the overdamping, means that the perturbative correction is larger than the PTEP 2016, 063I01 K. Nagao and I. Danshita 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 longwavelength 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.