Manipulating Goldstone modes via the superradiant light in a bosonic lattice gas inside a cavity

The mechanism of spontaneous symmetry breaking is crucial for understanding phase transitions and is also widely used to study the associated emergence of new particles and excitations Huang (1987). It is well known that when a continuous symmetry is spontaneously broken in nonrelativistic theories, there appear Nambu-Goldstone (NG) modes Nambu (1960); Goldstone (1961). The dispersion relations of which are either linear (type I) or quadratic (type II), where the numbers of distinct types of NG modes satisfy the Nielsen-Chadha inequality Nielsen and Chadha (1976). Thanks to recent experimental developments, the Goldstone modes have been studied in various condensed matter Bloch (1930); Chubukov et al. (1994); Podolsky and Sachdev (2012); Ye and Jiang (2007) and ultra-cold atomic systems  Ernst et al. (2010a); Papp et al. (2008); Steinhauer et al. (2002); Stenger et al. (1999); Stamper-Kurn et al. (1999); Kozuma et al. (1999); Yu et al. (2013); Pan et al. (2020).

A gas of bosonic atoms in an optical lattice has been reversibly tuned between superfluids (SF) and insulating ground states by varying the strength of the periodic potential Greiner et al. (2002); Stöferle et al. (2004). It provides an ideal platform to study the spontaneous symmetry breaking induced elementary excitations. Not only has the gapless Goldstone mode been found to exhaust all of the spectral weight in the weakly-interacting limit, but also the Higgs amplitude mode has been detected in the strongly-interacting regime near the transition between SF and Mott state at commensurate fillings Sachdev (2011); Endres et al. (2012); Pekker and Varma (2015).

Recently, a complimentary approach using ultracold atoms inside a cavity unveils the effect of cavity-mediated interactions, which results in the observation of a rich phase diagram with Mott insulator, superfluid, supersolid (SS) and charge-density-wave (CDW) phases Landig et al. (2016); Klinder et al. (2015). The low-energy excitations associated with above various symmetry-breaking phases have been investigated within the framework of mean-field approach, where new features, such as the softening of particle- and hole-like excitations, have been explored Dogra et al. (2016). The quantum fluctuation has not been taken into account within such mean-field calculations. However, it has been shown to play an important role in the strongly-interacting regime, leading to intriguing physical phenomena, such as the appearance of Higgs mode in the excitation spectrum of the Bose-Hubbard model Freericks and Monien (1994, 1996).

In this work, a systematic strong-coupling expansion method has been developed by performing two successive Hubbard-Stratonovich transformations Sengupta and Dupuis (2005) for the system of ultracold bosons in an optical lattice strongly coupled to a single mode of a high-finesse optical cavity Landig et al. (2016); Klinder et al. (2015). The effect of quantum fluctuations can thus been systematically explored in the strongly-interacting regime, which has not been studied previously Dogra et al. (2016); Chen et al. (2016). Interestingly, we find that the quantum fluctuations at the SS-to-CDW transition cause an unexpected effect of the modulation of density imbalance, neglected in the previous studies  Dogra et al. (2016). While being small, this modulation results in dramatic changes of the behavior of low-energy excitations. A free switching between the odd (type I) and even (type II) power energy-momentum dispersion NG modes can be achieved along the phase boundary. This prediction would pave a new way to manipulate NG modes in a cold atom based system.

Effective model — Same as the ETH experiment Landig et al. (2016), here we consider a Bose-Einstein condensate (BEC), such as ⁸⁷Rb atoms, subjected to an ultrahigh-finesse optical cavity, where three mutually orthogonal standing waves form a highly anisotropic $3$D optical lattice $V_{L}(\mathbf{r})=V_{y}\cos^{2}(q_{0}y)+V_{2D}[\cos^{2}(q_{0}x)+\cos^{2}(q_{0}z)])$ with the lattice depth $V_{y}\gg V_{2D}$ and the wave vectors of laser fields denoted by $q_{0}$. The BEC is thus split into a stack of $2$D layers, where the bosonic atoms are exposed to an additional potential $V_{c}(\mathbf{r})=\hbar\bar{\eta}(\hat{a}+\hat{a}^{{\dagger}})\cos(q_{0}x)\cos(q_{0}z)-\hbar(\Delta_{c}-U_{0}\cos^{2}(q_{0}z))\hat{a}\hat{a}^{{\dagger}}$ in the $xz$-plane formed by the coherent scattering between one free space lattice and one intracavity optical standing wave with the same wave-vector $q_{0}$ Landig et al. (2016); Klinder et al. (2015). $\bar{\eta}$ is the two-photon Rabi frequency determining the scattering rate. $\hat{a}$ and $\hat{a}^{{\dagger}}$ are the annihilation and creation operators for the cavity photon, respectively. $\Delta_{c}=\omega_{p}-\omega_{c}$ describes the discrepancy between pumping light (lattice standing-wave) frequency $\omega_{p}$ and the cavity resonance frequency $\omega_{c}$. $U_{0}$ captures the maximum light shift per atom resulting from the effect of the dispersive shift of the cavity resonance frequency Landig et al. (2016).

Since the coherent scattering of light between the lattice and the cavity mode creates a dynamical checkerboard superlattice for the atoms, the effective Hamiltonian describing the atomic dynamics dressed by the cavity field can be expressed as

	$\displaystyle H$	$\displaystyle=$	$\displaystyle\sum_{\mathbf{r},\mathbf{r}^{\prime}}-\mathrm{T}_{\sigma\sigma^{\prime}}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)\psi_{\sigma}^{\dagger}(\mathbf{r})\psi_{\sigma^{\prime}}\left(\mathbf{r}^{\prime}\right)-\sum_{\mathbf{r},\sigma}\mu_{\sigma}n_{\sigma}(\mathbf{r})$		(1)
		$\displaystyle+$	$\displaystyle\frac{U}{2}\sum_{\mathbf{r},\sigma}n_{\sigma}(\mathbf{r})(n_{\sigma}(\mathbf{r})-1)-\delta_{c}|\alpha|^{2},$

where $\psi_{\sigma}(\mathbf{r})$ is the bosonic atom field. $\sigma=e,o$ stands for even and odd sites, respectively, describing distinct sublattices of the checkerboard lattice. Here we choose the nearest neighbor two sites as one unit cell and $\mathbf{r}$ ($\mathbf{r}^{\prime}$) is the lattice index capturing the location of the unit cell. The expression of the hopping matrix $\mathrm{T}_{\sigma\sigma^{\prime}}$ is given in the supplementary material (SM). $U$ captures the strength of the repulsive interaction between bosonic atoms determined by the effective $s$-wave scattering length, which can be tuned by means of the Feshbach resonance and lattice depth. $n_{\sigma}(\mathbf{r})=\psi_{\sigma}^{\dagger}(\mathbf{r})\psi_{\sigma}(\mathbf{r})$ is the density operator. The onsite energies $\mu_{e}=\mu-\eta\left(\alpha^{*}+\alpha\right)$ and $\mu_{o}=\mu+\eta\left(\alpha^{*}+\alpha\right)$ are introduced for even and odd sites, respectively, where $\mu$ is the chemical potential and $\eta$ is the energy shift due to the two-photon Rabi process. $\delta_{c}=\hbar(\Delta_{c}-\delta)$ is the energy detuning between the cavity and the pumping lights, where $\delta$ is the dispersive shift of the cavity due to the BEC Landig et al. (2016). $\alpha$ is the mean value of the cavity field, which can be determined by the steady equation $i\hbar\partial_{t}\alpha=<[\hat{\alpha},{\tilde{\mathbf{H}}}]>-i\kappa\alpha=0$ with ${\tilde{\mathbf{H}}}=-\sum_{\mathbf{r},\mathbf{r}^{\prime}}\mathrm{T}_{\sigma\sigma^{\prime}}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)\psi_{\sigma}^{\dagger}(\mathbf{r})\psi_{\sigma^{\prime}}\left(\mathbf{r}^{\prime}\right)-\mu\sum_{\mathbf{r},\sigma}n_{\sigma}(\mathbf{r})+\frac{U}{2}\sum_{\mathbf{r},\sigma}n_{\sigma}(\mathbf{r})(n_{\sigma}(\mathbf{r})-1)+\eta(\hat{a}+\hat{a}^{{\dagger}})\sum_{\mathbf{r}}(n_{e}(\mathbf{r})-n_{o}(\mathbf{r}))-\hbar(\Delta_{c}-\delta)\hat{a}^{{\dagger}}\hat{a}$ in the regime of large decay rate $\kappa$ Baumann et al. (2010) and leads to the following relation

	$\displaystyle\alpha$		$\displaystyle=\frac{\eta\sum_{\mathbf{r}}\left[<n_{e}(\mathbf{r})>-<n_{o}(\mathbf{r})>\right]}{\delta_{c}+i\kappa}$		(2)
			$\displaystyle\overset{\left|\Delta_{c}\right|\gg\kappa,\left|\delta\right|}{\approx}\frac{N\eta\left(\left\langle n_{e}\right\rangle-\left\langle n_{o}\right\rangle\right)}{\hbar\Delta_{c}}$

where $<n_{e}(\mathbf{r})>$ and $<n_{o}(\mathbf{r})>$ are the average density of bosonic atoms at even and odd sites, respectively. $2N$ is the total lattice site.

Path integral approach — In the following, we develop a strong-coupling expansion of the model Hamiltonian in Eq. (1), which extends the treatment in Bose-Hubbard model Sengupta and Dupuis (2005), where both superfluid and Mott phases can be captured in the same time by this method. Through introducing the complex field $\psi_{\sigma}$, the thermodynamic properties of the system can be obtained from the partition function $Z$ as a functional integral with the action $S\left[\psi_{\sigma}^{*},\psi_{\sigma}\right]=\int_{0}^{\beta}d\tau\left\{\sum_{\sigma,\mathbf{r}}\psi_{\sigma}^{*}(\mathbf{r})\partial_{\tau}\psi_{\sigma}(\mathbf{r})+H\left[\psi^{*},\psi\right]\right\}$. Here, $\tau$ is a imaginary time and $\beta=1/k_{B}T$ is inverse temperature. We then perform a Hubbard-Stratonovich transformation through introducing the auxiliary field $\phi_{\sigma}$ to decouple the intersite hopping in the action and obtain

	$\displaystyle Z=\int D\left[\psi_{\sigma}^{*},\psi_{\sigma},\phi_{\sigma}^{*},\phi_{\sigma}\right]\exp\{-\int_{0}^{\beta}d\tau\sum_{\mathbf{r},\mathbf{r}^{\prime}}\phi_{\sigma}^{*}(\mathbf{r})T_{\sigma\sigma^{\prime}}^{-1}\phi_{\sigma^{\prime}}(\mathbf{r^{\prime}})$
	$\displaystyle+[\sum_{\sigma,\mathrm{r}}\int_{0}^{\beta}d\tau\phi_{\sigma}^{*}(\mathbf{r})\psi_{\sigma}(\mathbf{r})+c.c]-S_{0}\left[\psi_{\sigma}^{*},\psi_{\sigma}\right]\}$
	$\displaystyle=Z_{0}\int D\left[\phi_{\sigma}^{*},\phi_{\sigma}\right]e^{-\int_{0}^{\beta}d\tau\sum_{\mathbf{r},\mathbf{r}^{\prime}}\phi_{\sigma}^{*}(\mathbf{r})T_{\sigma\sigma^{\prime}}^{-1}\phi_{\sigma^{\prime}}(\mathbf{r^{\prime}})}$
	$\displaystyle\times\langle\exp[\sum_{\sigma,\mathrm{r}}\int_{0}^{\beta}d\tau\phi_{\sigma}^{*}(\mathbf{r})\psi_{\sigma}(\mathbf{r})+c.c]\rangle_{0}$
	$\displaystyle=Z_{0}\int D\left[\phi_{\sigma}^{*},\phi_{\sigma}\right]\exp\{-\int_{0}^{\beta}d\tau\sum_{\mathbf{r},\mathbf{r}^{\prime}}\phi_{\sigma}^{*}T_{\sigma\sigma^{\prime}}^{-1}\phi_{\sigma^{\prime}}+W\left[\phi_{\sigma}^{*},\phi_{\sigma}\right]\},$
			(3)

where $T_{\sigma\sigma^{\prime}}^{-1}$ represents the inverse hopping matrix. $S_{0}$ and $Z_{0}$ are the action and partition function in the limit of $t=0$. $\langle\cdots\rangle_{0}$ stands for averaging with $S_{0}$. $W\left[\phi_{\sigma}^{*},\phi_{\sigma}\right]=\ln\langle\exp[\sum_{\sigma,\mathrm{r}}\int_{0}^{\beta}d\tau\phi_{\sigma}^{*}(\mathbf{r})\psi_{\sigma}(\mathbf{r})+\mathrm{c.c.}]\rangle_{0}$ is the generation function linking to connected local Green’s functions $G^{Rc}$ through the relation [25] $W\left[\phi^{*},\phi\right]=\sum_{R=1}^{\infty}\frac{(-1)^{R}}{(R!)^{2}}\sum^{\prime}_{\sigma_{1}\cdots\sigma_{R}^{\prime}}G_{\left\{\sigma_{i},\sigma_{i}^{\prime}\right\}}^{Rc}\phi_{\sigma_{1}}^{*}\cdots\phi_{\sigma_{R}}^{*}\phi_{\sigma_{R}^{\prime}}\cdots\phi_{\sigma_{1}^{\prime}}$, where $\sum^{\prime}$ means that all the fields share the same value of the site index and ${\left\{\sigma_{i},\sigma_{i}^{\prime}\right\}}\equiv{\sigma_{1}\cdots\sigma_{R},\sigma_{1}^{\prime}\cdots\sigma^{\prime}_{R}}$. After doing a power expansion of $W\left[\phi^{*},\phi\right]$ to quartic order, we obtain that

	$\displaystyle S[\phi_{\sigma}^{*},\phi_{\sigma}]$	$\displaystyle=$	$\displaystyle\int_{0}^{\beta}d\tau\sum_{\mathbf{r},\mathbf{r}^{\prime}}\phi_{\sigma}^{*}(\mathbf{r})T_{\sigma\sigma^{\prime}}^{-1}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)\phi_{\sigma^{\prime}}\left(\mathbf{r}^{\prime}\right)-W\left[\phi_{\sigma}^{*},\phi_{\sigma}\right]$		(4)
		$\displaystyle=$	$\displaystyle\int_{0}^{\beta}d\tau\sum_{\mathbf{r},\mathbf{r}^{\prime}}\phi_{\sigma}^{*}(\mathbf{r})T_{\sigma\sigma^{\prime}}^{-1}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)\phi_{\sigma^{\prime}}\left(\mathbf{r}^{\prime}\right)$
		$\displaystyle+$	$\displaystyle\int d\tau_{1}d\tau_{2}\underset{\mathbf{r}}{\sum}G_{\sigma\sigma}(\mathbf{r,}\tau_{1}-\tau_{2}\mathbf{)}\phi_{\sigma}^{\ast}(\mathbf{r},\tau_{1})\phi_{\sigma}(\mathbf{r},\tau_{2})$
		$\displaystyle-$	$\displaystyle\frac{1}{2!}\int\underset{\alpha=1}{\overset{4}{\Pi}}d\tau_{\alpha}\underset{\mathbf{r}}{\sum}\chi_{\sigma\sigma^{\prime}}(\mathbf{r,}\tau_{1},\tau_{2},\tau_{3},\tau_{4})$
		$\displaystyle\times$	$\displaystyle\phi_{\sigma}^{\ast}(\tau_{1})\phi_{\sigma}(\tau_{2})\phi_{\sigma^{\prime}}^{\ast}(\tau_{3})\phi_{\sigma^{\prime}}(\tau_{4})+O(\phi^{6}),$

where $G$ is the local Green’s function and $\chi$ is the two-particle vertex in the local limit, i.e., $t=0$ (see details in SM). However, starting from the above action, it is inconvenient to calculate physical quantities, such as the excitation spectrum or the momentum distribution. These difficulties can be circumvented through performing another Hubbard-Stratonovich transition to decouple the hopping term Sengupta and Dupuis (2005). It is shown that the auxiliary field of this transformation has the same correlation functions as the original boson field (see details in SM). Therefore, we can use the same notation for both fields. The partition function can thus be expressed as

	$\displaystyle Z=Z_{0}\int D\left[\psi_{\sigma}^{*},\psi_{\sigma},\phi_{\sigma}^{*},\phi_{\sigma}\right]\exp\{\int_{0}^{\beta}d\tau\sum_{\mathbf{r},\mathbf{r}^{\prime}}\psi_{\sigma}^{*}T_{\sigma\sigma^{\prime}}\psi_{\sigma^{\prime}}$
	$\displaystyle-[\sum_{\sigma,\mathrm{r}}\int_{0}^{\beta}d\tau\psi_{\sigma}^{*}(\mathbf{r})\phi_{\sigma}(\mathbf{r})+\mathrm{c.c.}]+W\left[\phi_{\sigma}^{*},\phi_{\sigma}\right]\},$		(5)

Integrating out the $\phi_{\sigma}$ field in Eq. (5), the effective action can be obtained

	$\displaystyle S[\psi_{\sigma}^{*},\psi_{\sigma}]=\int_{0}^{\beta}d\tau\sum_{\mathbf{r},r^{\prime}}\psi_{\sigma}^{*}(\mathbf{r},\tau)\left[-T_{\sigma\sigma^{\prime}}\left(\mathbf{r}-\mathbf{r}^{\prime}\right)\right]\psi_{\sigma^{\prime}}\left(\mathbf{r}^{\prime},\tau\right)$
	$\displaystyle+\int_{0}^{\beta}d\tau_{1}d\tau_{2}\sum_{\mathbf{r}}\psi_{\sigma}^{*}\left(\mathbf{r},\tau_{1}\right)\left[-G_{\sigma\sigma}^{-1}\left(\tau_{1}-\tau_{2}\right)\right]\psi_{\sigma}\left(\mathbf{r},\tau_{2}\right)$
	$\displaystyle+\frac{1}{2}g_{\sigma\sigma^{\prime}}\int_{0}^{\beta}d\tau\sum_{\mathbf{r}}\left|\psi_{\sigma}(\mathbf{r},\tau)\right|^{2}\left|\psi_{\sigma^{\prime}}(\mathbf{r},\tau)\right|^{2},$		(6)

where $g_{\sigma\sigma^{\prime}}$ captures the amplitude of the boson-boson interaction determined by the local two-particle vertex in its static limit (see details in SM).

Phase diagram under the saddle-point approximation — Starting from the above effective action, we first perform a saddle-point approximation to determine the ground state of our proposed system. The saddle-point action derived from Eq. (6) can be expressed as

	$\displaystyle S_{saddle}$	$\displaystyle=$	$\displaystyle-\bar{G}_{ee}^{-1}{\left|{{{\bar{\psi}}_{e}}}\right|^{2}}-\bar{G}_{oo}^{-1}{\left|{{{\bar{\psi}}_{o}}}\right|^{2}}-4t{{\bar{\psi}}_{o}}^{*}{{\bar{\psi}}_{e}}-4t{{\bar{\psi}}_{e}}^{*}{{\bar{\psi}}_{o}}$		(7)
		$\displaystyle-$	$\displaystyle\frac{1}{4}\frac{{{{\bar{\chi}}_{ee}}}}{{\bar{G}_{ee}^{4}}}{\left|{{{\bar{\psi}}_{e}}}\right|^{4}}-\frac{1}{4}\frac{{{{\bar{\chi}}_{oo}}}}{{\bar{G}_{oo}^{4}}}{\left|{{{\bar{\psi}}_{o}}}\right|^{4}},$

where $\bar{G}_{\sigma\sigma}\equiv G_{\sigma\sigma}(i\omega=0)$ and $\bar{\chi}_{\sigma\sigma}\equiv\chi_{\sigma\sigma}(i\omega=0)$ with $G_{\sigma\sigma}(i\omega)$ and $\chi_{\sigma\sigma}(i\omega)$ being Fourier transforms of the single-particle Green’s function and the two-particle vertex, respectively. The saddle-point value $\bar{\psi}_{\sigma}$ can be obtained from minimizing $S_{saddle}$ and we obtain

	$\displaystyle n_{s,e}$	$\displaystyle=$	$\displaystyle\frac{\bar{G}_{ee}^{-1}+4tc^{-1}}{g_{ee}}\text{ \ if }\bar{G}_{ee}^{-1}+4tc^{-1}>0$
	$\displaystyle n_{s,e}$	$\displaystyle=$	$\displaystyle 0\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ otherwisde,}$		(8)
	$\displaystyle n_{s,o}$	$\displaystyle=$	$\displaystyle\frac{\bar{G}_{oo}^{-1}+4tc}{g_{oo}}\text{ \ \ \ \ \ \ if }\bar{G}_{oo}^{-1}+4tc>0$
	$\displaystyle n_{s,o}$	$\displaystyle=$	$\displaystyle 0\text{ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ otherwise},$		(9)

where $n_{s,o(e)}=|\bar{\psi}_{\sigma}|^{2}$ is the superfluid density of the pseudospin ${\sigma}$ component and $c=\bar{\psi}_{e}/\bar{\psi}_{o}$. By utilizing an iterative algorithm to solve a complete set of self-consistent equations (2), (8) and (9), superfluid order parameters $\bar{\psi}_{\sigma}$ for even and odd lattice sites and the mean value of cavity field $\alpha$ can be determined. As defined in Eq. (2), $\alpha$ characterizes the average imbalance between even and odd lattice sites, which can be used as the order parameter describing the CDW order.

We display the phase diagram as a function of the lattice depth $V_{2D}$ and detuning $\Delta_{c}$ in Fig. 1(a), to compare with the realistic experiments, such as the ETH experimental setup Landig et al. (2016). There are four different phases as shown in the phase diagram, which consists of a superfluid (SF) phase, a Mott insulator (MI) state, a charge-density-wave (CDW) state and a supersolid (SS) phase. The difference among these four phases can be captured by order parameters defined above, for example, as shown in Fig. 1(b). A SF phase is characterized by finite and equal superfluid order parameters and vanishing even-odd imbalance indicated by $\alpha=0$. It is quite different in the SS phase, where a finite even-odd imbalance and nonzero superfluid order parameters are observed. In both MI and CDW states, superfluid order parameters vanish. However, the presence of a finite even-odd imbalance $\alpha\neq 0$ distinguishes between MI and CDW states. The obtained phase diagram is consistent with other mean-field calculations Dogra et al. (2016); Chen et al. (2016). However, the path integral approach constructed here can provide a systematic way beyond the mean-field approach for investigating the effect of high-order fluctuations. The unexpected results driven by the fluctuations will be unveiled below.

Manipulating Goldstone modes at the SS-CDW transition — Since typically fluctuations will play a dominate role in determining physical properties at the critical region, to explore the unexpected results arising from the fluctuations, we construct the critical theory describing the phase transition region as shown in Fig. 1(a). By performing both a spatial and a temporal gradient expansion of the action in Eq. (6), as well as the cumulant expansion in powers of $\psi_{\sigma}$ Sengupta and Dupuis (2005); Fisher et al. (1989); Faccioli and Salasnich (2018), the effective action can be rewritten as

$\displaystyle\begin{aligned} S=&\int\nolimits_{0}^{\beta}d\tau\int d\mathbf{r\{}\left[\psi_{e}^{\ast}(\mathbf{r},\tau)r_{e}\psi_{e}(\mathbf{r},\tau)+\psi_{o}^{\ast}(\mathbf{r},\tau)r_{e}\psi_{o}(\mathbf{r},\tau)-\psi_{o}^{\ast}(\mathbf{r},\tau)4t\psi_{e}(\mathbf{r},\tau)-\psi_{e}^{\ast}(\mathbf{r},\tau)4t\psi_{o}(\mathbf{r},\tau)\right]\\ &+\left[K_{e1}\psi_{e}^{\ast}(\mathbf{r},\tau)\partial_{\tau}\psi_{e}(\mathbf{r},\tau)+K_{o1}\psi_{o}^{\ast}(\mathbf{r},\tau)\partial_{\tau}\psi_{o}(\mathbf{r},\tau)\right]+\left[K_{e2}\left|\partial_{\tau}\psi_{e}(\mathbf{r},\tau)\right|^{2}+K_{o2}\left|\partial_{\tau}\psi_{o}(\mathbf{r},\tau)\right|^{2}\right]\\ &+\left[K_{3}\nabla\psi_{e}^{\ast}(\mathbf{r},\tau)\nabla\psi_{o}(\mathbf{r},\tau)+K_{3}\nabla\psi_{o}^{\ast}(\mathbf{r},\tau)\nabla\psi_{e}(\mathbf{r},\tau)\right]+\frac{g_{ee}}{2}\left|\psi_{e}(\mathbf{r},\tau)\right|^{4}+\frac{g_{oo}}{2}\left|\psi_{o}(\mathbf{r},\tau)\right|^{4}\\ &+\frac{g_{eo}}{2}\left|\psi_{e}(\mathbf{r},\tau)\right|^{2}\left|\psi_{o}(\mathbf{r},\tau)\right|^{2}+\frac{g_{oe}}{2}\left|\psi_{o}(\mathbf{r},\tau)\right|^{2}\left|\psi_{e}(\mathbf{r},\tau)\right|^{2}\},\end{aligned}$

(10)

where $r_{e}=-{G}_{ee}^{-1}(i\omega=0)$, $r_{o}=-{G}_{oo}^{-1}(i\omega=0)$, $K_{e1}=\frac{\partial G_{ee}^{-1}\left(i\omega\right)}{\partial(i\omega)}|_{i\omega=0}$, $K_{o1}=\frac{\partial G_{oo}^{-1}\left(i\omega\right)}{\partial(i\omega)}|_{i\omega=0}$. $K_{e(o)2}$ and $K_{3}$ are the coefficients of the second-order temporal and spatial derivatives, respectively, which can be expressed in terms of the Green’s function (see SM for details).

Next, we introduce $\eta_{\sigma}(\mathbf{r},\tau)=\psi_{\sigma}(\mathbf{r},\tau)-\bar{\psi}_{\sigma}$ to describe the fluctuations. After expanding the action in Eq. (10) to the quadric order of fluctuation fields $\eta_{\sigma}$, we get

$\displaystyle S=\frac{1}{2}\underset{\mathbf{k},\omega}{\sum}\tilde{\eta}^{{\dagger}}({\mathbf{k},i\omega})M({\mathbf{k},i\omega})\tilde{\eta}({\mathbf{k},i\omega}),$

(11)

where $\tilde{\eta}^{{\dagger}}=[\eta_{e}^{\ast}(\mathbf{k},i\omega),\eta_{e}(\mathbf{-k},-i\omega),\eta_{o}^{\ast}(\mathbf{k},i\omega),$ $\eta_{o}(\mathbf{-k},-i\omega)]$ and the matrix elements of $M$ can be expressed as

$\displaystyle\begin{split}&M_{11}=r_{e}-i\omega K_{e1}+K_{e2}\omega^{2}+2g_{ee}\bar{\psi}_{e}^{2}+\frac{g_{eo}}{2}\bar{\psi}_{o}^{2}+\frac{g_{oe}}{2}\bar{\psi}_{o}^{2},\\ &M_{22}=r_{e}+i\omega K_{e1}+K_{e2}\omega^{2}+2g_{ee}\bar{\psi}_{e}^{2}+\frac{g_{eo}}{2}\bar{\psi}_{o}^{2}+\frac{g_{oe}}{2}\bar{\psi}_{o}^{2},\\ &M_{33}=r_{o}-i\omega K_{o1}+K_{o2}\omega^{2}+2g_{oo}\bar{\psi}_{o}^{2}+\frac{g_{eo}}{2}\bar{\psi}_{e}^{2}+\frac{g_{oe}}{2}\bar{\psi}_{e}^{2},\\ &M_{44}=r_{o}+i\omega K_{o1}+K_{o2}\omega^{2}+2g_{oo}\bar{\psi}_{o}^{2}+\frac{g_{eo}}{2}\bar{\psi}_{e}^{2}+\frac{g_{oe}}{2}\bar{\psi}_{e}^{2},\\ &M_{12}=M_{21}=g_{ee}\bar{\psi}_{e}^{2},\quad M_{34}=M_{43}=g_{oo}\bar{\psi}_{o}^{2},\\ &M_{13}=M_{24}=M_{31}=M_{42}\\ &=-4t+K_{3}k^{2}+\frac{g_{eo}}{2}\bar{\psi}_{e}\bar{\psi}_{o}+\frac{g_{oe}}{2}\bar{\psi}_{o}\bar{\psi}_{e},\\ &M_{14}=M_{23}=M_{32}=M_{41}=\frac{g_{eo}}{2}\bar{\psi}_{e}\bar{\psi}_{o}+\frac{g_{oe}}{2}\bar{\psi}_{o}\bar{\psi}_{e}.\end{split}$

Therefore, the single particle excitation spectrum at the critical region as shown in Fig. 1(a), can be obtained by solving the equation $det[M({\mathbf{k},i\omega})]=0$.

Here we focus on the transition from SS to CDW. In SS phase, there are two superfluid order parameters $\bar{\psi}_{\sigma}\equiv\sqrt{n_{s,\sigma}}e^{i\theta_{\sigma}}$ for even and odd lattice sites, where the amplitude and phase of the order parameter emerge as two independent degree of freedoms, instead of being conjugate to each other. The phase fluctuation leads to the Goldstone mode, while the amplitude fluctuation leads to the Higgs mode. When further considering the coupling between even and odd lattice sites (hopping term between them), the two Higgs modes arising from the amplitude fluctuation will couple to each other and split into two new Higgs modes (two higher modes obtained from $det[M({\mathbf{k},i\omega})]=0$). While the two Goldstone modes resulting from the phase fluctuation will also couple to each other and split into one gapless Goldstone mode associated with the overall phase $\theta_{o}+\theta_{e}$ plus a gapped mode linking to the relative phase $\theta_{o}-\theta_{e}$.

At the transition from SS to CDW, we find the new feature of the lowest Goldstone mode excitation. Previous studies show that the slight modulation of density imbalance induced by the superradiant cavity light has a negligible effect on the single particle excitation spectrum Dogra et al. (2016). Here we firstly show that at the transition from SS to CDW, even such negligible modulation of density imbalance will result in a dramatic effect on the excitations. Distinct from the mean-field approach Dogra et al. (2016); Chen et al. (2016), such as employing the Gutzwiller ansatz, the critical theory constructed in Eq. (10) can unveil the dominate role of fluctuations. It is found that tuning the system along the phase boundary between SS and CDW as shown in Fig 1(a) can freely switch between two types of Goldstone modes in its single particle excitation, i.e., type I and type II with odd and even power energy-momentum dispersion (as shown in Fig. 2), respectively.

To further understand the underlying physics, we analytically solve the equation $det[M({\mathbf{k},i\omega})]=0$ and find that in the long wave limit the lowest single particle excitation (positive branch) can be approximately expressed as (See details in SM)

$\displaystyle E^{+}_{Goldstone}\simeq\sqrt{I_{2}}k+\sqrt{I_{4}}k^{2}.$

(12)

Therefore, the switching from type I to type II Goldstone modes can be achieved by changing the coefficients $I_{2}$ and $I_{4}$. Around the transition from SS to CDW, it is shown that the quantum fluctuations make even the slight modulation of density imbalance play a dominant role in determining $I_{2}$ and $I_{4}$. As shown in Fig. 3, the slight modulation of density imbalance leads to two distinct regions: (I) $I_{2}\neq 0$, (II) $I_{2}=0$ and $I_{4}\neq 0$, which corresponds type I and type II Goldstone modes, respectively. To distinguish these two different types of Goldstone modes, we calculate the density of states (DOS) for the lowest single particle excitation as $N(E)=\frac{1}{N}\sum_{\mathbf{k}}\delta(E-E^{+}_{Goldstone})$, which is directly related to the Bragg spectroscopy signal Pippan et al. (2009). As shown in Fig. 2, with linear dispersion (type I), we find $N(E)\propto E$ when $E\rightarrow 0$. It is distinct from the quadratic dispersion (type II), where N(E) is a constant when $E\rightarrow 0$. The experimental advances in the momentum-resolved Bragg spectroscopy in optical lattices Ernst et al. (2010b) make the detection of this signal accessible.

Conclusion — By performing two successive Hubbard-Stratonovich transformations, we have developed an effective field theory to study a bosonic lattice gas inside a cavity in the strongly interacting regime. Through taking into account the quantum fluctuations, we firstly find that the slight modulation of density imbalance, neglected in the previous studies, leads to dramatic changes of the behavior of low-energy excitations. The switching between type I and type II Goldstone modes can be driven along SS-CDW transition. Our findings would bridge the cavity light and strongly interacting quantum matters and may open up a new direction in this field.

Acknowledgement — This work is supported by NSFC (Grant No. 12074305, 11774282, 11950410491), the National Key Research and Development Program of China (2018YFA0307600), Cyrus Tang Foundation Young Scholar Program and the Fundamental Research Funds for the Central Universities(H.W., S.L., M.A and B.L.) and by the AFOSR Grant No. FA9550-16-1-0006, the MURI-ARO Grant No. W911NF17-1-0323 through UC Santa Barbara, and the Shanghai Municipal Science and Technology Major Project through the Shanghai Research Center for Quantum Sciences (Grant No. 2019SHZDZX01) (W.V.L.). We also thank the HPC platform of Xi’An Jiaotong University, where our numerical calculations was performed.

Supplementary Material:
Manipulating Goldstone modes via the superradiant light in a bosonic lattice gas inside a cavity