Non-equilibrium dynamics of ultracold lattice bosons inside a cavity

The study of dynamics of quantum many-body system is one of the most exciting frontiers in modern condensed matter physics Polkovnikov et al. (2011); Aoki et al. (2014). The so-called Kibble-Zurek mechanism (KZM) has been used to understand certain universal features of the dynamics across a continuous phase transition both in classical Kibble (1976, 1980); Zurek (1985, 1993, 1996) and quantum Damski (2005); Zurek et al. (2005); Polkovnikov (2005); Dziarmaga (2005); Mukherjee et al. (2007); Sen et al. (2008); Dziarmaga (2010); Polkovnikov et al. (2011); Aoki et al. (2014); Cincio et al. (2009) realms. Based on conventional scaling laws near the criticality, KZM predicts a universal scaling for resulting density of defects when the system crosses a second-order phase transition by a linear quench. Recently, tremendous amount of efforts in both theoretical and experimental studies have been triggered to explore such quench dynamics, for instance in liquid crystals Bowick et al. (1994); Digal et al. (1999), quantum optical systems Xu et al. (2014), superconducting films Maniv et al. (2003); Golubchik et al. (2010), trapped ions Chandran et al. (2012); Francuz et al. (2016); Bermudez et al. (2009, 2010); Dziarmaga and Zurek (2014); Gardas et al. (2017), as well as ultracold gases Navon et al. (2015); Chen et al. (2011); Braun et al. (2015); Shimizu et al. (2018a, b, c).

In particular, ultracold gases in optical lattices provide a versatile tool for simulating and studying dynamics of quantum many-body physics Jaksch et al. (1998); Greiner et al. (2002); Bloch et al. (2008); Cirac and Zoller (2012); 201 (2012), including the early studies of single-component bosonic lattice gases Shimizu et al. (2018a, b, c), Rydberg-dressed atoms Zhou et al. (2020); Keesling et al. (2019) and dipolar bosons Sable et al. (2021) in optical lattices. Lots of interesting properties of these systems, such as verifying the Kibble-Zurek (KZ) scaling law between the defect formation and quench rate Chen et al. (2011); Braun et al. (2015); Anguez et al. (2016); Clark et al. (2016), have already been well studied. At the same time, there has also been a great interest in exploring the deviation from the KZ scaling in quench dynamics across the critical regime of the phase transition. Various new effects, such as new statistics of the resulting defects del Campo (2018); Gomez-Ruiz et al. (2020), fantastic evolution of correlation functions Roychowdhury et al. (2021) and unusual saturation of the defect density Donadello et al. (2016); Ko et al. (2019); Goo et al. (2021); del Campo et al. (2010); Liu et al. (2018); Chesler et al. (2015) have been unveiled.

In this work, motivated by recent progresses in experimental investigation of the bosonic lattice gas inside a cavity Landig et al. (2016); Klinder et al. (2015), where the effect of cavity-mediated interactions result in the observation of a rich equilibrium phase diagram including Mott insulator (MI), superfluid (SF), supersolid (SS) and charge-density-wave (CDW) phases, we study its non-equilibrium quench dynamics. This set-up offers new possibilities for exploring various non-equilibrium dynamics via varying the parameters of the Hamiltonian crossing different quantum phase transitions (QPT). And here we focus on studying the quench dynamics across the continuous phase transition between the CDW and SS phases. It is found that there is a threshold of the quench rate. Below that threshold, various physical quantities, such as the correlation length and vortex density, show a good agreement with the KZ scaling law. While above that threshold, the saturation of the defect density has been found, indicating a deviation of the KZ scaling at fast quench.

This paper is organized as follows. In Section II, through employing the static Gutzwiller (GW) method, the equilibrium phase diagram is obtained, which is consistent with other mean-field calculations. In Section III, the protocol of quench is introduced and the non-equilibrium quench dynamics across the continuous phase transition between the CDW and SS phases have been studied by the time-dependent Gutzwiller (tGW) method. The scaling laws of the correlation length and the number of defects have been investigated. It is shown that there is a threshold of the quench rate separating two distinct scaling regimes. One is satisfied with the KZ scaling law, the other is characterized with the saturation of defect density, showing a deviation from the KZ scaling. We further show that such distinct scaling behaviors can be understood through comparing the relaxation time and the quench rate.

Let us consider load a Bose-Einstein condensate (BEC) of ${}^{87}\mathrm{Rb}$ into a highly anisotropic $3$D optical lattice coupled to an ultrahigh-finesse optical cavity, being similar as the ETH experimental setup Landig et al. (2016). Since the coherent scattering of light between the lattice and cavity modes creates a dynamical checkerboard superlattice for the atoms Baumann et al. (2010); Mottl et al. (2012); Wang et al. (2022), the effective Hamiltonian describing the atomic dynamics dressed by the cavity field can be expressed as

where $\hat{b}_{i}$ and $\hat{b}_{i}^{{\dagger}}$ are the annihilation and creation operators for bosonic atoms at the lattice site ${\mathbf{r}}_{i}$. ${e}$ and ${o}$ refer to the even and odd sites of the lattice, respectively. $\hat{n}_{i}=\hat{b}_{i}^{{\dagger}}\hat{b}_{i}$ is the on-site particle number operator. $N$ is the total number of lattice sites. $J$ captures the tunneling amplitude between nearest neighbors and $\mu$ is the chemical potential. The strength of the on-site repulsion is labeled as $U$, which can be tuned through the Feshbach resonance technics. $U_{L}=-N\hbar M_{0}^{2}\eta^{2}/\Delta_{c}$ describes the strength of the cavity-mediated interaction between atoms, where $M_{0}$ is the overlap between the Wannier function and the lattice potential $(\propto\cos\frac{2\pi x}{\lambda}\cos\frac{2\pi y}{\lambda})$ and $\eta$ is the two-photon Rabi frequency of the scattering process. $\Delta_{c}=\omega_{p}-\omega_{c}$ represents the dispersive shift of pumping frequency $\omega_{p}$ and cavity resonance frequency $\omega_{c}$.

The cavity-mediated interaction ($U_{L}$ term in Eq. (1)) favors an overall even-odd sites imbalance, which can be characterized by the defined charge-density-wave (CDW) order parameter as $\theta=2\langle(\sum\limits_{{i\in e}}{{{\hat{n}}_{{i}}-}}\sum\limits_{{i\in o}}{{{\hat{n}}_{{i}}}})\rangle/N$, where $\langle...\rangle$ means the expectation value in the ground state. The superfluid order parameters for the even and odd lattice sites are introduced as $\phi_{e/o}=\left\langle b_{e/o}\right\rangle$, respectively. Then, we employ the static Gutzwiller (GW) method to study the equilibrium zero-temperature phase diagram of the system. Consider starting with the GW ansatz

where $i\in e,o$ and $f_{n_{e}}^{(e)}$ , $f_{n_{o}}^{(o)}$ are the variational coefficients. $n_{max}$ is the maximum number of particles at each site. Then, the above defined order parameters can be rewritten as

Here, we take $n_{max}=7$ for both even and odd sites, which is proved to be large enough for the study of our proposed system. The average number density of atoms is determined through the relation $\left\langle n_{i}\right\rangle=\underset{n_{i}=0}{\overset{n_{max}}{\sum}}\left|f_{n_{i}}^{(i)}\right|^{2}n_{i}$. By minimizing the expectation ground state energy $\langle{{\Psi_{gw}}}|H|{{\Psi_{gw}}}\rangle$, the variational coefficients in $|{{\Psi_{gw}}}\rangle$ can be determined (see Appendix A for details). Therefore, order parameters defined above can be obtained through the GW method. The equilibrium zero-temperature phase diagram can thus be obtained, such as shown in Fig. 1(a). The difference among distinct phases can be captured by the order parameters, for example, as shown in Fig. 1(b). In the supersolid (SS) phase, a finite even-odd imbalance and nonzero superfluid order parameters coexist, which is distinct from a superfluid (SF) phase where the even-odd imbalance is vanished and there is finite and equal superfluid order parameters. In the charge-density-wave (CDW) state, superfluid order parameters vanish. And the presence of a finite even-odd imbalance distinguishes the CDW from a Mott insulator (MI) state. Our phase diagram is in good agreement with other existing studies Chen et al. (2016); Dogra et al. (2016); Niederle et al. (2016); Sundar and Mueller (2016); Flottat et al. (2017); Himbert et al. (2019); Wang et al. (2022).

In this section, we will study the dynamics across the phase transition from CDW to SS as shown in Fig. 1(a). The quench protocol is constructed through tuning the tunneling amplitude $J$ by the following relation

where $J_{i}$ and $J_{f}$ are the initial and final values of the quench tunneling amplitude, which satisfy the relations $J_{i}<J_{c1}$ and $J_{c1}<J_{f}<J_{c2}$ with $J_{c1}$ and $J_{c2}$ being the equilibrium phase transition points of CDW to SS and SS to SF as shown in Fig. 1(a), respectively. Here $J_{f}+J_{i}=2J_{c1}$ means that the system approaches the quantum critical point between CDW and SS at $t=0$. $\tau_{Q}$ is the quench time, which captures the quench rate. We then employ the time-dependent Gutzwiller (tGW) method for investigating the real-time dynamics of our proposed system via varying $J(t)$ as in Eq.(4). In the practical calculation, we use the interaction strength $U$ as the unit of energy and $\hbar/U$ is used as the unit of time $t$. In the tGW approximation, the Hamiltonian in Eq. (1) is devised into a single-site Hamiltonian and can be approximated as $H_{GW}$ (see details in Appendix B). Then, we introduce the tGW wave function as

The coefficients $f_{n_{i}}^{(i)}(t)$ can be determined by solving the following Schr$\ddot{o}$dinger equation $i\hbar\partial_{t}\left|\Psi_{tgw}(t)\right\rangle=H_{GW}(t)\left|\Psi_{tgw}(t)\right\rangle$ for various initial states corresponding to the CDW phase and $H_{GW}(t)$ comes from replacing $J$ by $J(t)$ in $H_{GW}$. Therefore, the coefficients of tGW wave function can be determined through the following relation

Then, we use the fourth-order Runge-Kutta method to study the time evolution of the system and calculate various physical quantities, such as the SF order parameter $\Phi$, correlation length ${\xi}$ and vortex number $N_{v}$, which are defined as follows

where $\Omega_{i}=\frac{1}{4}[\sin(\theta_{i+\hat{e}_{x}}-\theta_{i})+\sin(\theta_{i+\hat{e}_{x}+\hat{e}_{y}}-\theta_{i+\hat{e}_{x}})+\sin(\theta_{i+\hat{e}_{y}}-\theta_{i+\hat{e}_{x}+\hat{e}_{y}})+\sin(\theta_{i}-\theta_{i+\hat{e}_{y}})]$ with $\theta_{i}$ being the phase of superfluid order parameter $\phi_{i}$ and $\hat{e}_{x}(\hat{e}_{y})$ being the unit vector along the $x(y)$ direction. The system size is chosen as $100\times 100$ in most of our numerical calculations and we have verified that this is a sufficiently large size where the size-dependent effect disappears.

As shown in Fig. 2, the typical behavior of SF order parameter $\Phi$ as the function of time $t$ is demonstrated. At $t=0$, the system crosses the quantum critical point between CDW and SS at $J_{c1}$ and the relaxation time of the system diverges. The dynamics is thus frozen and the state of the system cannot follow the change in the hopping amplitude. This scenario persists till the transition time $\hat{t}$, which is determined by the relation $|\Phi(\hat{t})|=2|\Phi({t=0})|$. Then, $\Phi$ develops very rapidly. After that rapid increase, $\Phi$ starts to fluctuate and such a oscillatory trend sets the average value tend towards the steady state value. We also study the defects formation during the course of the quench dynamics by investigating the vortex number $N_{v}$ defined above, as an indicator of the defect density. When the quench dynamics in the CDW region, the vortex density is high, which is caused by the effect of phase fluctuations introduced during the initial state preparation. Then, at $t=0$ the system crosses the quantum critical point and enters into the SS regime where spontaneous symmetry breaking occurs. The SF order parameters acquire finite values and small domains are gradually formed (shown in Fig. 3). However, such domains are formed where there is a phase coherence within, but not between the domains. As the quench is continued further, due to the annihilation of vortex-antivortex pairs, domains are merged and the size of domain becomes larger. The number of that thus decreases. After approaching the steady state, phase coherence is established in the system and the vortex number approaches zero.

Next, we will investigate the behavior of quench dynamics with different quench rate $(\propto 1/\tau_{Q})$. To study the scaling laws for the above defined quantities, such as correlation length ${\xi}$ and vortex number $N_{v}$, with respect to the quench time, we monitor ${\xi}$ and $N_{v}$ at the transition time $\hat{t}$, because $\hat{t}$ captures the freeze out time indicating the moment that the system catches the quench speed beginning to evolve quickly and can thus be linked to the universal scaling law connected with the relaxation time. As shown in Fig. 4, there are two different regimes separating distinct scaling behaviors for both ${\xi}$ and $N_{v}$ as a function of $\tau_{Q}$. For slow quench, such as $\tau_{Q}$ from $\tau_{Q}=10$ to $\tau_{Q}=400$ (Fig. 4), it is shown that both ${\xi}$ and $N_{v}$ satisfy a fairly good scaling law as ${\xi}\propto{\tau_{Q}}^{b}$ and ${N_{v}}\propto{\tau_{Q}}^{-d}$ with $b\approx 0.33$ and $d\approx 0.64$, respectively. The relation $d\approx 2b$ holds. Therefore, the KZ hypothesis predicted power law scaling is satisfied in the slow quench regime. Moreover, the KZM also predicts that the transition time $\hat{t}$ should have the scaling law as $\hat{t}\sim\tau_{Q}^{\nu z/(1+\nu z)}$ with $\nu$ and $z$ being the critical exponent of the equilibrium correlation length and the dynamical critical exponent, respectively. Through extracting the exponent from the scaling of $\hat{t}$, such as shown in Fig. 4(c) and utilizing the relation $b=\frac{\nu}{1+\nu z}$ and $d=\frac{2\nu}{1+\nu z}$, $\nu$ and $z$ can be estimated. We get $z\approx 0.98$ and $\nu\approx 0.49$, where $z$ is fairly close to that expected value $z=1$ from the 3D XY model, but $\nu$ dose not coincide with $\nu=0.672$ predicated in the 3D XY model Burovski et al. (2006), indicating that the spatial inhomogeneity effects the critical behavior Shimizu et al. (2018c). While considering fast quench, as shown in Fig. 4, it is found that the vortex defect density demonstrates a plateau whose value is a constant independent of the quench rate and ${\xi}(\hat{t})$ also has the similar behavior. It indicates that there is a violation of the KZ scaling in the fast quench regime.

In the following, we will make a detailed analysis of the above distinct scaling behaviors and determine the critical quench rate ($\propto 1/\tau_{Qc}$) separating two different quench regimes. To understand the quench dynamics here, we first study the relaxation time of the system. Since the transition time $\hat{t}$ captures the moment that the system catches the quench speed and begin to evolve quickly, it can be linked to the relaxation time $\tau$ through the following relation

Therefore, through solving the above equation, the relaxation time can be obtained. As shown in Fig. 5, we display the relaxation time as a function of the hopping amplitude $J$. Through comparing the relaxation time and the quench rate, distinct scaling behaviors during the dynamic procedure can be understood. Let us consider fixing $J_{f}$ at a certain value at $J_{f}=\bar{J}$. Then, various quench procedures with different quench time $\tau_{Q}$ can be distinguished by a critical quench time $\tau_{Qc}$ determined by the relation $J(\hat{t})=\bar{J}$, for instance, as indicated by the blue-dashed line in Fig. 5. Therefore, $\tau_{Qc}$ can be determined by the slope of that blue-dashed line decides a threshold of the quench rate separating two distinct scaling regime. Let us first consider the fast quench above that threshold. For example, as indicated by the red-solid line in Fig. 5, the linear quench ends before the time that such line (stands for $J(t)$) meets the $\tau(J)$ line. Therefore, the transition time $\hat{t}$ for any fast quench above the threshold of the quench rate will be the same and be fixed at $\hat{t}=\tau(\bar{J})$, being independent of the values of $\tau_{Q}$. The domains of defect for different quench rates will thus have the same average size, which naturally explains the plateau of defect density appears in rapidly quenched dynamics. While considering the slow quench below that threshold, for instance, as indicated by the green-solid line in Fig. 5, the $J(t)$ line intersects the $\tau(J)$ line before the end of linear quench. It means that the system cannot follow the change of $J(t)$ and the dynamics is frozen, which will persist till the transition time marking the impulse period of the quench. After the evolution time is equal to the relaxation time, the system evolves promptly. Therefore, the adiabatic-impulse approximation is valid and the KZ scaling is satisfied. Furthermore, the critical quench time $\tau_{Qc}$ can also be determined by the intersection point of the un-smooth transition from the plateau to KZM power law scaling of the correlation length or vortex density, for instance as shown by the dashed-line in the inset of Fig. 5. The obtained critical quench time $\tau_{Qc}$ is consistent with the way through comparing the relaxation time and the quench rate as mentioned above.

In conclusion, we have examined the quench dynamics crossing a continuous phase transition between CDW and SS phases of a bosonic lattice gas with cavity-mediated interactions. By using the time-dependent Gutzwiller method, we find that various physical quantities, such as the correlation length and vortex density, show two distinct scaling regimes with respect to the quench time. There is a threshold of quench rate. Below that, the scaling behavior is in good agreement with the KZ scaling law. When above that threshold, the scaling behavior is characterized with the saturation of various quantities, indicating a deviation from the KZ scaling. Furthermore, our proposal can be viably simulated experimentally in ultracold gases, since such a system has already been implemented in experiments. It thus paves an alternative way for investigating the dynamics of the system combined with the cavity light and interacting quantum matters.