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]:

(1)
$H=−J∑⟨i,j⟩(ai†aj+h.c.)+U2∑i(ai†ai−n¯)2−δμ∑i(ai†ai−n¯),$
where $$a_i^{\dagger}$$ and $$a_i$$ denote the creation and annihilation operators for a boson at site $$i$$ and satisfy the bosonic commutation relations $$[a_i,a_j^{\dagger}]=\delta_{i,j}$$, $$[a_i,a_j]=[a_i^{\dagger},a_j^{\dagger}]=0$$. Here, $${\bar n} \equiv \langle a_i^{\dagger}a_i \rangle$$ is the filling factor. The symbol $$\langle i,j\rangle$$ denotes the nearest-neighbor sites. The parameters $$J$$, $$U$$, and $$\delta \mu$$ are the hopping energy, the on-site interaction energy, and the chemical potential. Note that at a sufficiently large filling $${\bar n} \gg 1$$, $$\delta \mu =0$$ corresponds to the commensurate case, i.e., integer filling $${\bar n}\in{\bf N}$$. This definition simplifies the following derivation of the effective model to describe the system near the Mott phase with filling $$\bar n$$. In this paper, we focus on the high-filling case.

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:

(2)
$|n¯+1⟩i=t1i†|vac⟩,|n¯⟩i=t0i†|vac⟩,|n¯−1⟩i=t−1i†|vac⟩,$
where $$|{\rm vac}\rangle$$ is a vacuum of the Schwinger bosons. These operators satisfy the bosonic commutation relations $$[t_{mi},t_{nj}^{\dagger}]=\delta_{m,n}\delta_{i,j}$$$$(m,n=-1,0,1)$$ and we impose a local constraint $$\sum_{m=-1}^{m=1} t_{mi}^{\dagger}t_{mi}=1$$ on them to eliminate nonphysical states, e.g., $$t_{1i}^{\dagger}t_{0i}^{\dagger}|{\rm vac}\rangle$$. 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^{\dagger}$$ can be expressed in terms of the Schwinger bosons (2):
(3)
$ai†=n¯+1t1i†t0i+n¯t0i†t−1i.$

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]:

(4)
$Heff=−Jn¯2∑⟨i,j⟩(Si+Sj−+h.c.)+U2∑i(Siz)2−h∑iSiz,$
where $$h=\delta \mu$$. The spin operators are written in terms of the Schwinger bosons (2):
(5)
$Si+=2(t1i†t0i+t0i†t−1i),Si−=(Si+)†,Siz=t1i†t1i−t−1i†t−1i,$
and obey the standard SU(2) algebra. For $${\bar n}Jz/U \gg 1$$, the ground state is the XY ferromagnetic ordered state where the U(1) symmetry is spontaneously broken. Here $$z=2d$$ is the coordination number and $$d=3$$ is the spatial dimension of the system. For $${\bar n}Jz/U \ll 1$$, the ground state has U(1) symmetry and the mean-field wave function is given by $$\prod_i | S=1,m_z=0\rangle_i \equiv \prod_i|{\bar n}\rangle_i$$. Notice that, in the $$h=\delta \mu =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}^{\dagger}$$ ($$t_{1i}$$) and hole excitations $$t_{-1i}^{\dagger}$$ ($$t_{-1i}$$).

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

(6)
$b0i=cos(θ0/2)t0i+sin(θ0/2)(t1i+t−1i)/2,b1i=sin(θ0/2)t0i−cos(θ0/2)(t1i+t−1i)/2,b2i=(t1i−t−1i)/2,$
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 $$[b_{mi},b^{\dagger}_{nj}]=\delta_{m,n}\delta_{i,j}$$, $$[b_{mi},b_{nj}]=[b^{\dagger}_{mi},b^{\dagger}_{nj}]=0$$ and the local constraint $$\sum_{m=0}^{m=2}b^{\dagger}_{mi}b_{mi}=1$$. The angle of transformation is given by $$\theta_0={\rm tan}^{-1}(\sqrt{1-u^2}/u)$$, where $$u\equiv U/(4J\bar{n}z)$$ 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 limit $${\bar n} \gg 1$$, the superfluid to Mott insulator transition occurs at $$u=1\equiv 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}^{\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

(7)
$bm†b0=bm†1−b1†b1−b2†b2,≈bm†−12bm†b1†b1−12bm†b2†b2+⋯,$
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\ge3$$. At lower dimensions, the thermal fluctuations are so strong that they destroy the long-range order of the mean-field ground state [48–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:

(8)
$Heff≡Heff(0)+Heff(1)+Heff(2)+Heff(3)+⋯,$
where the index $$l$$ in $$H_{\rm eff}^{(l)}$$ means $$l$$th order with respect to the operators $$b_{mi}^{\dagger}$$ and $$b_{mi}$$ ($$m=1,2$$). Here, we perform the Fourier transformation of the operators defined as
(9)
$bmi†=1N∑kbmk†e−ixi⋅k,bmi=1N∑kbmkeixi⋅k,$
where $$N$$ is the system size and the vector $${\bf x}_i$$ denotes site $$i$$. The summation with respect to momentum $$\bf k$$ is taken over the first Brillouin zone. Then, each term $$H_{\rm eff}^{(l)}$$ can be written as the following equations, respectively:
(10)
$Heff(0)=N(U2s2−2Jn¯zs2c2),Heff(1)=2Jn¯zscN(c2−s2−u)(b10†+b10)=0,$

(11)
$Heff(2)=∑kB1(k)b1k†b1k+∑kB2(k)(b1k†b1−k†+b1kb1−k+b1kb1k†+b1k†b1k)+∑kC1(k)b2k†b2k+∑kC2(k)(b2k†b2−k†+b2kb2−k−b2kb2k†−b2k†b2k),$

(12)
$Heff(3)=∑k1∑k2∑k3D1(k1,k2,k3)(b2k1†b2k2b1k3+h.c.)+∑k1∑k2∑k3D2(k1,k2,k3)(b1k1†b1k2b1k3+h.c.)+∑k1∑k2∑k3D3(k1,k2,k3)(b1k1†b2k2b2k3−b2k1b1k2†b2k3†+h.c.),$
where $$c\equiv{\rm cos}(\theta_0/2)$$ and $$s\equiv{\rm sin}(\theta_0/2)$$. The explicit expressions of some coefficients in Eqs. (11) and (12) are given by
$B1(k)=2Jn¯z,B2(k)=−Jn¯z2u2γkC1(k)=Jn¯z(1+u),C2(k)=Jn¯z4(1+u)γkD1(k1,k2,k3)=−1NJn¯zu1−u2γk3δk1−k2+k3D2(k1,k2,k3)=−1N2Jn¯zu1−u2γk3δk1−k2+k3D3(k1,k2,k3)=1NJn¯z21−u2γk3δk1−k2−k3$
where $$\gamma_{\bf k}=\sum_{\bf e}e^{-i{\bf k}\cdot{\bf e}}/z=\sum_{s=1}^d{\rm cos}(k_s)/d$$.

Bogoliubov transformation

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

(13)
$bmk=umkβmk+vm−k∗βm−k†,bm−k†=um−k∗βm−k†+vmkβmk,$
where $$m=\{1,2\}$$. The new operators obey the same commutation relations $$[\beta_{m{\bf p}},\beta^{\dagger}_{n{\bf q}}]=\delta_{n,m}\delta_{{\bf p},{\bf q}}$$ as the previous operator such that the coefficients $$u_{m{\bf k}}$$ and $$v_{m {\bf k}}$$ satisfy a relation $$|u_{m{\bf k}}|^2-|v_{m{\bf k}}|^2=1$$. Assuming that $$u_{m{\bf k}}$$ and $$v_{m {\bf k}}$$ are real numbers and imposing the condition that $$H_{\rm eff}^{(2)}$$ in the new representation has no anomalous term, we determine the coefficients
(14)
$u1k=2−u2γk41−u2γk+12,v1k=2−u2γk41−u2γk−12,$

(15)
$u2k=2−γk41−γk+12,v2k=−sgn(γk)2−γk41−γk−12,$
where $${\rm sgn}(x)$$ is the sign function. After the Bogoliubov transformation, Eq. (4) reads
(16)
$Heff=∑m=12∑kEmkβmk†βmk+Heff(3).$

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

(17)
$E1k=2Jn¯z1−u2γk,$

(18)
$E2k=Jn¯z(1+u)1−γk.$

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].

Fig 1.

Independent interaction terms contained in $$H_{\rm eff}^{(3)}$$. The solid and dashed lines represent the propagator of the Higgs mode and that of the NG mode. The incoming lines into and outgoing lines from the vertex correspond to the annihilation and creation of quasiparticles.

Fig 1.

Independent interaction terms contained in $$H_{\rm eff}^{(3)}$$. The solid and dashed lines represent the propagator of the Higgs mode and that of the NG mode. The incoming lines into and outgoing lines from the vertex correspond to the annihilation and creation of quasiparticles.

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

(19)
$E1k≈Δ2+ch2|k|2,$

(20)
$E2k≈cng|k|,$
where $$c_{\rm h}=2J{\bar n}\sqrt{z}u$$ and $$c_{\rm ng}=J{\bar n}\sqrt{z}(1+u)$$. Similarly, if we assume that $$u \neq u_c$$, the coefficients (14) and (15) become
(21)
$u1k≈2−u2+Δ2Δ+O(k2),v1k≈2−u2−Δ2Δ+O(k2),$

(22)
$u2k≈u+14cngk(1+cngku+1+O(k2)),v2k≈−u+14cngk(1−cngku+1+O(k2)).$

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

(23)
$⟨βmk(iωn)β¯mk(iωn)⟩≡∫D(β,β¯)βmk(iωn)β¯mk(iωn)exp(−Seff)∫D(β,β¯)exp(−Seff),$
where $$\omega_n=2\pi n/\beta$$$$(n=0,\pm1,\pm2,\ldots)$$ is the Matsubara frequency, $$\beta = T^{-1}$$ is the inverse temperature, $$\beta_{m{\bf k}}(i \omega_n)$$ and the conjugate $${\bar \beta_{m{\bf k}}}(i \omega_n)=(\beta_{m{\bf k}}(i \omega_n))^*$$ are complex-valued field variables at $$(\omega_n,{\bf k})$$, and $${\cal D}(\beta,{\bar \beta})$$ is a measure of the integrations. The effective action $${\cal S}_{\rm eff}={\cal S}^{(2)}_{\rm eff}+{\cal S}^{(3)}_{\rm eff}$$ is derived from the effective model (16) and the quadratic action $${\cal S}^{(2)}_{\rm eff}$$ is given by
(24)
$Seff(2)=∑m=12∑n∑k(−iωn+Emk)β¯mk(iωn)βmk(iωn).$

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

(25)
$Γ1k=π2∑k1∑k2|Mk,k1,k2|2(1+fB(E2k1)+fB(E2k2))δ(E1k−E2k1−E2k2),$

(26)
$Γ2k=π∑k1∑k2|Mk1,k,k2|2(fB(E2k2)−fB(E1k1))δ(E2k−E1k1+E2k2),$
where $$\delta(x)$$ is the $$\delta$$ function, $$f_{\rm B}(x)=1/(e^{\beta x}-1)$$ is the Bose distribution function, and the matrix elements $${\cal M}_{{\bf k}_1,{\bf k}_2,{\bf k}_3}$$ are given by
(27)
$Mk1,k2,k3=1Nδk1,k2+k3[−Jn¯zu1−u2γk1(u1k1+v1k1)(u2k2v2k3+v2k2u2k3)+Jn¯z21−u2γk2(u1k1u2k3−v1k1v2k3)(u2k2−v2k2)+Jn¯z21−u2γk3(u1k1u2k2−v1k1v2k2)(u2k3−v2k3)].$

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).

Fig 2.

Contribution to the self-energy of (a) the Higgs mode $$\Sigma_{1}(i\omega_n;{\bf k})$$ and (b) the NG mode $$\Sigma_{2}(i\omega_n;{\bf k})$$. Other diagrams do not yield a nonzero contribution to the damping rates.

Fig 2.

Contribution to the self-energy of (a) the Higgs mode $$\Sigma_{1}(i\omega_n;{\bf k})$$ and (b) the NG mode $$\Sigma_{2}(i\omega_n;{\bf k})$$. Other diagrams do not yield a nonzero contribution to the damping rates.

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

(28)
$Γ1k=0=33/2Jn¯z232π(1+u)1−u2cothβΔ4.$

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$$.

Fig 3.

Dependence on $$u$$ of the damping rate $$\Gamma_{1{\bf k}={\bf 0}}/(J\bar{n}z)$$ at fixed temperatures. The solid lines represent the analytic expression (28) and the dotted lines the numerical data, where $$N=500^3$$.

Fig 3.

Dependence on $$u$$ of the damping rate $$\Gamma_{1{\bf k}={\bf 0}}/(J\bar{n}z)$$ at fixed temperatures. The solid lines represent the analytic expression (28) and the dotted lines the numerical data, where $$N=500^3$$.

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.

Fig 4.

Dependence on $$u$$ of the ratio $$\Gamma_{1{\bf k=0}}/\Delta$$ at fixed temperatures. The solid and dotted lines stand for the analytic and numerical data. The choice of parameters is the same as that in Fig. 3.

Fig 4.

Dependence on $$u$$ of the ratio $$\Gamma_{1{\bf k=0}}/\Delta$$ at fixed temperatures. The solid and dotted lines stand for the analytic and numerical data. The choice of parameters is the same as that in Fig. 3.

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

(29)
$Γ2k=−(1+u)4(1−u2)8πcng2ch2k2∫klkud|k2|(fB(E2k+E2k2)−fB(E2k2)),$
where $$k\equiv|{\bf k}|$$. The upper and lower bounds of the integration $$k_u$$ and $$k_l$$ are given by
(30)
$ku=−k+Δ2cng2−ch2,kl=1cng2−ch2{−(cng2+ch2)k+4ch2cng2k2+Δ2(cng2−ch2)}.$

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

(31)
$Γ2k=35/2(1+u)2(1−u)4β2πu2k2log1−e−βcng(kl+k)1−e−βcng(ku+k)1−e−βcngku1−e−βcngkl.$

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

(32)
$z|u−uc|≪k,$
meaning that Eq. (31) is valid only near the critical point. It is obvious that in the zero-temperature limit ($$\beta \rightarrow \infty$$) the damping rate (31) vanishes because there is no thermal excitation.

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:

(33)
$Γ2k≈35/2(1+u)2(1−u)42πu2kcngfB(cngkl).$

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.

Fig 5.

Dependence on $$u$$ of the damping rate $$\Gamma_{2{\bf k}}/(J\bar{n}z)$$ at fixed temperatures and initial momenta. The solid lines represent the analytic results of Eq. (31) and the dots represent the numerical evaluation of Eq. (26) without the long-wavelength approximation. The system size taken in the numerical calculations is $$N=500^3$$. We set $$(T/(J{\bar n}z),k)$$ = $$(0.1,0.1)$$ (a), $$(0.3,0.3)$$ (b), $$(0.1,0.3)$$ (c), and $$(0.3,0.3)$$ (d). Notice that we do not show the region where $$T > \Delta \sim T_c$$.

Fig 5.

Dependence on $$u$$ of the damping rate $$\Gamma_{2{\bf k}}/(J\bar{n}z)$$ at fixed temperatures and initial momenta. The solid lines represent the analytic results of Eq. (31) and the dots represent the numerical evaluation of Eq. (26) without the long-wavelength approximation. The system size taken in the numerical calculations is $$N=500^3$$. We set $$(T/(J{\bar n}z),k)$$ = $$(0.1,0.1)$$ (a), $$(0.3,0.3)$$ (b), $$(0.1,0.3)$$ (c), and $$(0.3,0.3)$$ (d). Notice that we do not show the region where $$T > \Delta \sim T_c$$.

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.

Fig 6.

Schematic illustration of the Landau damping process of the NG mode. The red and blue lines stand for the Higgs and NG branches. The red and blue circles represent the excited Higgs and NG modes. (a) corresponds to a vicinity of the critical point such that $$c_{\rm ng}k \simeq \Delta/2$$. (b) corresponds to a region far from the critical point such that $$c_{\rm ng}k \ll \Delta$$.

Fig 6.

Schematic illustration of the Landau damping process of the NG mode. The red and blue lines stand for the Higgs and NG branches. The red and blue circles represent the excited Higgs and NG modes. (a) corresponds to a vicinity of the critical point such that $$c_{\rm ng}k \simeq \Delta/2$$. (b) corresponds to a region far from the critical point such that $$c_{\rm ng}k \ll \Delta$$.

Fig 7.

(a) Dependence on $$u$$ of the ratio $$\Gamma_{2{\bf k}}/(c_{\rm ng}k)$$ at $$T=0.3J{\bar n}z$$ and two values of $${\bf k}$$. (b) Dependence on $$k$$ of the ratio $$\Gamma_{2{\bf k}}/(c_{\rm ng}k)$$ at $$T=0.3J{\bar n}z$$, $${\bf k}/k=(1,0,0)$$, and three values of $$u$$. For both cases, we set $$N=500^3$$ for the numerical calculations.

Fig 7.

(a) Dependence on $$u$$ of the ratio $$\Gamma_{2{\bf k}}/(c_{\rm ng}k)$$ at $$T=0.3J{\bar n}z$$ and two values of $${\bf k}$$. (b) Dependence on $$k$$ of the ratio $$\Gamma_{2{\bf k}}/(c_{\rm ng}k)$$ at $$T=0.3J{\bar n}z$$, $${\bf k}/k=(1,0,0)$$, and three values of $$u$$. For both cases, we set $$N=500^3$$ for the numerical calculations.

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

(A1)
$|Ω(θ,η,φ,χ)⟩=∏i{cos(θ2)t0i†+eiηsin(θ2)(eiφsin(χ2)t1i†+e−iφcos(χ2)t−1i†)}|vac⟩,$
where $$\theta$$, $$\eta$$, $$\varphi$$, and $$\chi$$ are variational parameters [27,28]. Notice that the variational wave function $$|\Omega(\theta,\eta,\varphi,\chi)\rangle$$ is normalized such as $$\langle\Omega|\Omega\rangle =1$$. The operators $$t_{0i}^{\dagger}$$, $$t_{1i}^{\dagger}$$, and $$t_{-1i}^{\dagger}$$ are the Schwinger bosons defined in Eq. (2) and represent the state vectors $$|{\bar n}\rangle_i=t_{0i}^{\dagger}|{\rm vac}\rangle$$, $$|{\bar n}+1\rangle_i=t_{1i}^{\dagger}|{\rm vac}\rangle$$, and $$|{\bar n}-1\rangle_i=t_{-1i}^{\dagger}|{\rm vac}\rangle$$, respectively. The wave function (A1) allows us to calculate averages of the physical quantities within the mean-field approximation. Indeed, the average of the Hamiltonian (1) $$E \equiv \langle\Omega|H|\Omega\rangle$$ and order parameter $$\Psi^* \equiv \langle\Omega|a^{\dagger}|\Omega\rangle$$ can be easily calculated as
(A2)
$E(θ,η,φ,χ)N=(U2+δμcosχ)sin2(θ2)−Jn¯z4sin2θ(1+n¯−1sin2(χ2)+1+n¯−1sinχcos2η),$

(A3)
$Ψ∗=12sinθe−iφ(n¯+1e−iηsin(χ2)+n¯eiηcos(χ2)),$
where $$N$$ is the system size. Equations (A1), (A2), and (A3) play an essential role in identifying the canonical transformation (6).

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

(A4)
$|Ω(θ0,0,φ0,π/2)⟩=∏i{cos(θ02)t0i†+12sin(θ02)(eiφ0t1i†+e−iφ0t−1i†)}|vac⟩.$

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

$b0i†|vac〉={cos(θ02)t0i†+12sin(θ02)(eiφ0t1i†+e−iφ0t−1i†)}|vac〉,b1i†|vac〉={sin(θ02)t0i†−12cos(θ02)(eiφ0t1i†+e−iφ0t−1i†)}|vac〉,b2i†|vac〉=12(eiφ0t1i†−e−iφ0t−1i†)|vac〉,$
in the leading order of $$\delta \theta_i$$, $$\delta \varphi_i$$. The operators $$b_{1i}^{\dagger}$$ and $$b_{2i}^{\dagger}$$ represent excitations on the Higgs and NG branches. As one may easily confirm, the vectors $$b^{\dagger}_{0i}|{\rm vac}\rangle$$, $$b^{\dagger}_{1i}|{\rm vac}\rangle$$, and $$b^{\dagger}_{2i}|{\rm vac}\rangle$$ are orthogonal to each other and make the canonical transformation (6), where the angle of rotation is given by the parameter $$\theta_0={\rm tan}^{-1}(\sqrt{1-u^2}/u)$$, determined by the variation of the energy (A2).

References

[1]
Nambu
Y.
,
Phys. Rev.

117
,
648
(
1960
).
()
[2]
Goldstone
J.
,
Nuovo Cimento

19
,
154
(
1961
).
()
[3]
Higgs
P.W.
,
Phys. Rev. Lett.

13
,
508
(
1964
).
[4]
Volovik
G.E.
and
Zubkov
M.A.
,
J. Low Temp. Phys.

175
,
486
(
2014
).
[5]
Pekker
D.
and
Varma
C.M.
,
Annu. Rev. Cond. Matter Phys.

6
,
269
(
2015
).
[6]
Sooryakumar
R.
and
Klein
M.V.
,
Phys. Rev. Lett.

45
,
660
(
1980
).
[7]
Sooryakumar
R.
and
Klein
M.V.
,
Phys. Rev. B

23
,
3213
(
1981
).
[8]
Littlewood
P.B.
and
Varma
C.M.
,
Phys. Rev. Lett.

47
,
811
(
1981
).
[9]
Littlewood
P.B.
and
Varma
C.M.
,
Phys. Rev. B

26
,
4883
(
1982
).
[10]
Méasson
M.-A.
,
Gallais
Y.
,
Cazayous
M.
,
Clair
B.
,
Rodière
P.
,
Cario
L.
, and
Sacuto
A.
,
Phys. Rev. B

89
,
060503
(
2014
).
[11]
Matsunaga
R.
,
Y.I.
,
Makise
K.
,
Uzawa
Y.
,
Terai
H.
,
Wang
Z.
, and
Shimano
R.
,
Phys. Rev. Lett.

111
,
057002
(
2013
).
[12]
Matsunaga
R.
,
Tsuji
N.
,
Fujita
H.
,
Sugioka
A.
,
Makise
K.
,
Uzawa
Y.
,
Terai
H.
,
Wang
Z.
,
Aoki
H.
, and
Shimano
R.
,
Science

345
,
6201
(
2014
).
[13]
Sherman
D.
,
Pracht
U.S.
,
Gorshunov
B.
,
Poran
S.
,
Jesudasan
J.
,
Chand
M.
,
Raychaudhuri
P.
,
Swanson
M.
,
Trivedi
N.
,
Auerbach
A.
,
Schefﬂer
M.
,
Frydman
A.
, and
Dressel
M.
,
Nat. Phys.

11
,
188
(
2015
).
[14]
Cea
T.
,
Castellani
C.
, and
Benfatto
L.
, arXiv:1512.02544 [cond-mat.supr-con] [Search INSPIRE].
[15]
Cea
T.
,
Castellani
C.
,
Seibold
G.
, and
Benfatto
L.
,
Phys. Rev. Lett.

115
,
157002
(
2015
).
[16]
Rüegg
Ch.
,
Normand
B.
,
Matsumoto
M.
,
Furrer
A.
,
McMorrow
D.F.
,
Kramer
K.W.
,
Gudel
H.U.
,
Gvasaliya
S.N.
,
Mutka
H.
, and
Boehm
M.
,
Phys. Rev. Lett.

100
,
205701
(
2008
).
[17]
Merchant
P.
,
Normand
B.
,
Krämer
K.W.
,
Boehm
M.
,
McMorrow
D.F.
, and
Rüegg,
Ch.
Nat. Phys.

10
,
373
(
2014
).
[18]
Kuroe
H.
,
Takami
N.
,
Niwa
N.
,
Sekine
T.
,
Matsumoto
M.
,
F.
,
Tanaka
H.
, and
Takemura
K.
,
J. Phys.: Conf. Ser.

400
,
032042
(
2012
).
[19]
Demsar
J.
,
Biljaković
K.
, and
Mihailovic
D.
,
Phys. Rev. Lett.

83
,
800
(
1999
).
[20]
Schaefer
H.
,
Kabanov
V.V.
, and
Demsar
J.
,
Phys. Rev. B

89
,
045106
(
2014
).
[21]
Yusupov
R.
,
Mertelj
T.
,
Kabanov
V.V.
,
Brazovskii
S.
,
Kusar
P.
,
Chu
J.-H.
,
Fisher
I.R.
, and
Mihailovic
D.
,
Nat. Phys.

6
,
681
(
2010
).
[22]
Mertelj
T.
,
Kusar
P.
,
Kabanov
V.V.
,
Giraldo-Gallo
P.
,
Fisher
I.R.
, and
Mihailovic
D.
,
Phys. Rev. Lett.

110
,
156401
(
2013
).
[23]
Avenel
O.
,
Varoquaux
E.
, and
Ebisawa
H.
,
Phys. Rev. Lett.

45
,
1952
(
1980
).
[24]
Collett
C.A.
,
Pollanen
J.
, J.
Li
I.A.
,
Gannon
W.J.
, and
Halperin
W.P.
,
J. Low Temp. Phys.

171
,
214
(
2013
).
[25]
Bissbort
U.
,
Götze
S.
,
Li
Y.
,
Heinze
J.
,
Krause
J.S.
,
Weinberg
M.
,
Becker
C.
,
Sengstock
K.
, and
Hofstetter
W.
,
Phys. Rev. Lett.

106
,
205303
(
2011
).
[26]
Endres
M.
,
Fukuhara
T.
,
Pekker
D.
,
Cheneau
M.
,
Schauß
P.
,
Gross
C.
,
Demler
E.
,
Kuhrm
S.
, and
Bloch
I.
,
Nature

487
,
454
(
2012
).
[27]
Altman
E.
and
Auerbach
A.
,
Phys. Rev. Lett.

89
,
250404
(
2002
).
[28]
Huber
S.
,
Altman
E.
,
Büchler
H.P.
, and
Blatter
G.
,
Phys. Rev. B

75
,
085106
(
2007
).
[29]
Greiner
M.
,
Mandel
O.
,
Esslinger
T.
,
Hänsch
T.W.
, and
Bloch
I.
,
Nature

(London)
415
,
39
(
2002
).
[30]
Sachdev
S.
,
Quantum Phase Transition

(Cambridge University Press
,
Cambridge, UK
,
2011
),
2nd ed
.
[31]
Pollet
L.
and
Prokof’ev
N.
,
Phys. Rev. Lett.

109
,
010401
(
2012
).
[32]
Podolsky
D.
and
Sachdev
S.
,
Phys. Rev. B

86
,
054508
(
2012
).
[33]
Gazit
S.
,
Podolsky
D.
, and
Auerbach
A.
,
Phys. Rev. Lett.

110
,
140401
(
2013
).
[34]
Gazit
S.
,
Podolsky
D.
, and
Auerbach
A.
,
Phys. Rev. B

88
,
235108
(
2013
).
[35]
Chen
K.
,
Liu
L.
,
Deng
Y.
,
Pollet
L.
, and
Prokof’ev
N.
,
Phys. Rev. Lett.

110
,
170403
(
2013
).
[36]
Rançon
A.
and
Dupuis
N.
,
Phys. Rev. B

89
,
180501
(
2014
).
[37]
Liu
L.
,
Chen
K.
,
Deng
Y.
,
Endres
M.
,
Pollet
L.
, and
Prokof’ev
N.
,
Phys. Rev. B

92
,
174521
(
2015
).
[38]
Katan
Y.T.
and
Podolsky
D.
,
Phys. Rev. B

91
,
075132
(
2015
).
[39]
Rose
F.
,
Léonard
F.
, and
Dupuis
N.
,
Phys. Rev. B

91
,
224501
(
2015
).
[40]
Sherson
J.F.
,
Weitenberg
C.
,
Endres
M.
,
Cheneau
M.
,
Bloch
I.
, and
Kuhr
S.
,
Nature

(London)
467
,
68
(
2010
).
[41]
Preiss
P.M.
,
Ma
R.
,
Tai
M.E.
,
Simon
J.
, and
Greiner
M.
,
Phys. Rev. A

91
,
041602
(
2015
).
[42]
Nakayama
T.
,
Danshita
I.
,
Nikuni
T.
, and
Tsuchiya
S.
,
Phys. Rev. A

92
,
043610
(
2015
).
[43]
Fisher
P.A.M.
,
Weichman
P.B.
,
Grinstein
G.
, and
Fisher
D.S.
,
Phys. Rev. B

40
,
546
(
1989
).
[44]
Jaksch
D.
,
Bruder
C.
,
Cirac
J.I.
,
Gardiner
C.W.
, and
Zoller
P.
,
Phys. Rev. Lett.

81
,
3108
(
1998
).
[45]
Schwinger
J.
,
On Angular Momentum
(
Dover
,
New York
,
2015
).
[46]
Holstein
T.
and
Primakoff
H.
,
Phys. Rev.

58
,
1908
(
1940
).
[47]
Auerbach
A.
,
Interacting Electrons and Quantum Magnetism
(
Springer
,
New York
,
1994
).
[48]
Hohenberg
P.C.
,
Phys. Rev.

158
,
383
(
1967
).
[49]
Mermin
N.D.
and
Wanger
H.
,
Phys. Rev. Lett.

17
,
1133
(
1966
).
[50]
Coleman
S.
,
Commun. Math. Phys.

31
,
259
(
1973
).
[51]
Pitaevskii
L.
and
Stringari
S.
,
Bose–Einstein Condensation
(
Oxford University Press
,
New York
,
2003
).
[52]
Podolsky
D.
,
Auerbach
A.
, and
Arovas
D.P.
,
Phys. Rev. B

84
,
174522
(
2011
).
[53]
Beliaev
S.T.
,
Sov. Phys. JETP

34
,
299
(
1958
).
[54]
Tsuchiya
S.
and
Grifﬁn
A.
,
Phys. Rev. A

70
,
023611
(
2004
).
[55]
Tsuchiya
S.
and
Grifﬁn
A.
,
Phys. Rev. A

72
,
053621
(
2005
).
[56]
Abrikosov
A.A.
,
Gorkov
L.P.
, and
Dzyaloshinski
I.E.
,
Methods of Quantum Field Theory in Statistical Physics
(
Dover
,
New York
,
1975
).
[57]
Popov
V.N.
,
Functional Integrals and Collective Excitations
(
Cambridge University Press
,
Cambridge, UK
,
1987
).
[58]
Negele
J.W.
and
Orland
H.
,
Quantum Many-Particle Systems
(
Westview
,
1998
).
[59]
Lifshitz
E.M.
and
Pitaevskii
L.
,
Statistical Physics, Part 2
(
Pergamon, Oxford
,
UK
,
1980
).
[60]
Pethick
C.J.
and
Smith
H.
,
Bose–Einstein Condensation in Dilute Gases
(
Cambridge University Press
,
Cambridge, UK
,
2008
).
[61]
Lifshitz
E.M.
and
Pitaevskii
L.
,
Physical Kinetics
(
Pergamon, Oxford
,
UK
,
1981
).
[62]
Trotzky
S.
,
Pollet
L.
,
Gerbier
F.
,
Schnorrberger
U.
,
Bloch
I.
,
Prokof’ev
N.V.
,
Svistunov
B.
, and
Troyer
M.
,
Nat. Phys.

6
,
998
(
2010
).
[63]
Afﬂeck
I.
and
Wellman
G.F.
,
Phys. Rev. B

46
,
8934
(
1992
).
[64]
Kozuma
M.
,
Deng
L.
,
Hagley
E.W.
,
Wen
J.
,
Lutwak
R.
,
Helmerson
K.
,
Rolston
S.L.
, and
Phillips
W.D.
,
Phys. Rev. Lett.

82
,
871
(
1999
).
[65]
Stenger
J.
,
Inouye
S.
,
Chikkatur
A.P.
,
Stamper-Kurn
D.M.
,
Pritchard
D.E.
, and
Ketterle
W.
,
Phys. Rev. Lett.

82
,
4569
(
1999
).
[66]
Ernst
P.T.
,
Götze
S.
,
Krauser
J.S.
,
Pyka
K.
,
Kuhmannn
D.-S.
,
Pfannkuche
D.
, and
Sengstock
K.
,
Nat. Phys.

6
,
56
(
2009
).
Whether or not the experimental signals observed in the superconductor Nb$$_{1-x}$$Ti$$_x$$N can be interpreted as the Higgs mode remains controversial. For more details, see Refs. [14,15].
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.