Quench dynamics of Rydberg-dressed bosons on two-dimensional square lattices

TL;DR

This paper investigates the non-equilibrium dynamics of Rydberg-dressed bosons on a 2D lattice, revealing universal early-stage behavior and interaction-dependent long-term evolution, with implications for experimental observation via time-of-flight imaging.

Contribution

It introduces a quench protocol in an extended Bose-Hubbard model with long-range interactions and analyzes the resulting dynamics, highlighting universal early behavior and the role of density waves.

Findings

Universal early-stage superfluid order parameter dynamics

Long-range interactions significantly influence late-time evolution

Density wave excitations are prominent at slow quench rates

Abstract

We study the dynamics of bosonic atoms on a two-dimensional square lattice, where atomic interactions are long-ranged with either a box or soft-core shape. The latter can be realized through laser dressing ground-state atoms to electronically excited Rydberg states. When the range of interactions is equal or larger than the lattice constant, the system is governed by an extended Bose-Hubbard model. We propose a quench process by varying the atomic hopping linearly across phase boundaries of the Mott insulator-supersolid and supersolid-superfluid phases. Starting from a Mott insulating state, the dynamical evolution of the superfluid order parameter exhibits a universal behaviour at the early stage, largely independent of interactions. The dynamical evolution is significantly altered by strong, long-range interactions at later times. Particularly, we demonstrate that density wave excitation is important when the quench rate is small. Moreover, we show that the quench dynamics can be analyzed through time-of-flight images, i.e., measuring the momentum distribution and noise correlations.

Equations

$$H=-J(t)\sum _{⟨ij⟩}{\hat{b}}_i^{†}{\hat{b}}_j+\frac{U}{2}\sum _i{\hat{n}}_i({\hat{n}}_i-1)+\sum _{i\ne j}{V}_{ij}{\hat{n}}_i{\hat{n}}_j,$$

$$i\frac{\partial }{\partial t}{f}_m^i$$

$$S(\mathbf{k})=\frac{1}{M^2}\sum _{j,l}{e}^{i\mathbf{k}\cdot ({\mathbf{r}}_j-{\mathbf{r}}_l)}⟨{\hat{n}}_j{\hat{n}}_l⟩,$$

$$E(\mathbf{k})=\sqrt{\epsilon (\mathbf{k}{)}^2+{\rho }_0\epsilon (\mathbf{k})\left(U+\tilde{V}(\mathbf{k})\right)},$$

$$\Sigma (\mathbf{k})={\rho }_0\frac{U+\tilde{V}(\mathbf{k})}{2-\cos (k_{x}d)-\cos (k_{y}d)}.$$

$$\widehat{ϵ}={\left(\frac{a_Q{\tau }_0}{J_s/U}\right)}^{\frac{1}{1+z\nu }},\hat{t}=\frac{J_s/U}{a_Q}{\left(\frac{a_Q{\tau }_0}{J_s/U}\right)}^{\frac{1}{1+z\nu }}.$$

$$n(\mathbf{k})\approx |{\tilde{\varphi }}_{\mathbf{k}}{|}^2,$$

$$C(\mathbf{q})\approx {\left|{\tilde{n}}_{\mathbf{q}}\right|}^2,$$

$$n(\mathbf{k})=|{\tilde{\varphi }}_{\mathbf{k}}{|}^2+\frac{1}{M}\sum _i\left(n_i-|{\varphi }_i{|}^2\right).$$

$$C(\mathbf{k},{\mathbf{k}}^{\prime })$$

$$\Delta (\mathbf{q})F(n-|\varphi {|}^2;\mathbf{q})-4X({\varphi }^{\ast };(n-|\varphi {|}^2)\varphi )+$$

$$2X({\varphi }^{\ast };\gamma -\eta {\varphi }^{\ast })+\sqrt{M}\Delta (\mathbf{q})F({\varphi }^{\ast };\mathbf{k})\tilde{F}(\varphi ;{\mathbf{k}}^{\prime })$$

$$\displaystyle\phantom{{}=}+2\sqrt{M}\Big{[}Y(\phi,\phi,\eta^{*}-\phi^{*2};\mathbf{k},\mathbf{k}^{\prime})+$$

$$\displaystyle\phantom{{}=}\quad Y(n-|\phi|^{2},\phi,\phi^{*};\mathbf{q},\mathbf{k})\Big{]}+\frac{C_{0}}{M},$$

Full text

I Introduction

In the past decades, there has been a growing interest in the study of ultracold atoms, which is largely driven by the unprecedented level of control over external trapping potentials, internal states and interactions between atoms with electromagnetic fields Jaksch et al. (1998); Stenger et al. (1999); Theis et al. (2004); Chin et al. (2010); Schachenmayer et al. (2010); Johnson and Rolston (2010); Anderson et al. (2011); Viteau et al. (2011); Macrì and Pohl (2014). Various lattice models Gross and Bloch (2017), such as the Bose-Hubbard model Fisher et al. (1989), have been studied and realized experimentally Greiner et al. (2002a). This has opened opportunities to probe static properties Jaksch and Zoller (2005); Lewenstein et al. (2007); Bloch et al. (2008), such as Mott insulator-superfluid phase transition Fisher et al. (1989); van Oosten et al. (2001); Greiner et al. (2002a); Zwerger (2003), spin-orbit coupling Wall et al. (2016); Kolkowitz et al. (2017); Zhang and Liu (2018), supersolidity Sengupta et al. (2005); Scarola et al. (2006); Batrouni et al. (2006); Menotti et al. (2007); Yi et al. (2007); Iskin and Freericks (2009); Iskin (2011); Li et al. (2012); Landig et al. (2016); Li et al. (2018), entanglement Alba et al. (2013); Islam et al. (2015), topology Goldman et al. (2016); Lohse et al. (2016), etc. In cold atom systems, many parameters can be manipulated and monitored dynamically. This allows exploring non-equilibrium dynamics in addition to steady states. Theoretical and experimental works have investigated Landau-Zener transitions Wu and Niu (2003); Tomadin et al. (2008); Deng et al. (2015), Kibble-Zurek mechanism Zurek et al. (2005); Dziarmaga et al. (2012); del Campo and Zurek (2014); Shimizu et al. (2018a, b, c); Weiss et al. (2018), transport Scherg et al. (2018); Brown et al. (2019); Fujiwara et al. (2019), and excitations of Higgs and Goldstone modes Léonard et al. (2017); Di Liberto et al. (2018).

Recently, growing interest has been spent on investigating the dynamics of Bose-Hubbard models (BHMs) driven by external periodic  Eckardt et al. (2005); Gaul et al. (2009) or linear fields Green and Sondhi (2005). A recent review can be found in Refs. Polkovnikov et al. (2011); Kennett (2013). A particularly interesting topic is the universal dynamics found in BHMs when linearly changing the hopping strength. Due to the Kibble-Zurek mechanism (KZM) Zurek et al. (2005); Dziarmaga et al. (2012); del Campo and Zurek (2014), the dynamics is frozen around the phase boundary while adiabatic away from it. Many quantities, such as correlation lengths and topological defects, exhibit universal behaviours.

In BHMs, the competition between the hopping and two-body on-site interactions Freericks and Monien (1994); Elesin et al. (1994); Fölling et al. (2005) leads to a Mott insulator (MI) to superfluid (SF) phase transition at some critical $J/U$. When the interaction length is greater than the lattice constant $d$, we obtain an extended Bose-Hubbard model (eBHM). Its ground state can have non-uniform, periodic densities, such as the density wave (DW), supersolid (SS), Haldane insulator, etc. Kühner et al. (2000); Rossini and Fazio (2012). Such a situation has been examined extensively using atoms with weak magnetic or electric dipole moments Góral et al. (2000); Griesmaier et al. (2005); Lu et al. (2011); Aikawa et al. (2012); Büchler et al. (2007); Micheli et al. (2005); Kotochigova and Tiesinga (2006); Gorshkov et al. (2011); Baier et al. (2016), where the dominant interaction is between the nearest-neighbour sites. Drastically different dynamics is found when quenching the eBHMs Shimizu et al. (2018b, c).

In this work, we will go beyond the nearest-neighbour interaction regime and examine the dynamics of eBHMs with even longer-range interactions. This situation can be realized by using Rydberg-dressed atoms confined in a two-dimensional optical lattice Balewski et al. (2014); Macrì and Pohl (2014); Jau et al. (2016); Zeiher et al. (2016); Mukherjee et al. (2016); Kościk and Sowiński (2018, 2019). Rydberg atoms have long lifetimes $10\sim 100~{}\mu\text{s}$, growing with $n^{3}$ where the principal quantum number $n\gg 1$. Large polarizabilities ($\propto n^{7}$) give strong and long-range interactions (e.g. van der Waals interactions $\propto n^{11}$) between Rydberg atoms. However, atomic hopping typically occurs on a much slower time scale than the Rydberg lifetime. One could instead weakly couple ground-state atoms to Rydberg states with a far detuned laser Henkel et al. (2010); Honer et al. (2010); Johnson and Rolston (2010); Pupillo et al. (2010); Chougale and Nath (2016). This results in a soft-core shape potential between Rydberg dressed ground-state atoms. Its range $r_{c}$ could extend several lattice constants before decaying significantly (see Fig. 1a). The interaction strength is nearly a constant for distance $r<r_{c}$ and decays quickly when $r>r_{c}$. Such an interaction potential may be approximated by a box potential, i.e. atomic interaction is a constant when $r\leq r_{c}$ and zero otherwise (Fig. 1a). Box-type interactions have been used to study the extended Bose-Hubbard model, with a focus on nearest-neighbour interactions ($r_{c}=d$) Kühner et al. (2000); Iskin (2011); Rossini and Fazio (2012). As the two types of potential share similar profiles in momentum space (Fig. 1b), we will show that the respective dynamics shows common features.

This paper is organized as follows. In Sec. II we introduce the eBHM of the Rydberg dressed atoms. Ground-state phase diagrams are calculated using the Gutzwiller method. We show the relation between the roton instability and density modulation by analysing the Bogoliubov spectra. In Sec. III, we discuss universal dynamics for the box interaction in the MI-SF phase transition and dynamics in the SS and SF phases. In Sec. IV, we propose that the quench dynamics can be measured through time-of-flight (TOF) density distributions Greiner et al. (2002b); Li et al. (2012); Baier et al. (2016) and covariance Altman et al. (2004); Fölling et al. (2005); Rom et al. (2006); Jeltes et al. (2007); Toth et al. (2008); Hu et al. (2010) at different probing times. In Sec. V, dynamics for the soft-core interaction is discussed. It is found that the dynamics of the eBHM with the soft-core interaction is largely similar to the box interaction. We conclude in Sec. VI.

II The extended Bose-Hubbard Model and ground-state phase

The Hamiltonian of the Rydberg dressed atoms in the 2D square lattice is given by an eBHM,

[TABLE]

where $\langle ij\rangle$ stands for the nearest-neighbour sites, and $J(t)$ is the time-dependent hopping. $U$ is the on-site interaction, including contributions from s-wave scattering and the long-range interaction. In this work, two types of long-range interactions are considered. The soft-core interaction has a form $V_{ij}=2V_{0}[1+(r_{ij}/r_{c})^{6}]^{-1}$, where $r_{ij}$ is the distance between site $i$ and $j$, $r_{c}$ is the soft-core radius, and $V_{0}$ is the interaction strength. In the case of a box potential, the interaction is given by $V_{ij}=V_{0}\Theta(r_{c}-r_{ij})$, where the Heaviside function $\Theta(r)$ defines the box length $r_{c}$. Here we will use $r_{c}$ to denote the interaction range for both box and soft-core interactions.

We employ the Gutzwiller approach to calculate the ground state and dynamics of the Hamiltonian Krauth et al. (1992); Seibold and Lorenzana (2001); Schiró and Fabrizio (2010); von Oelsen et al. (2011). The Gutzwiller approach is a mean-field method and predicts qualitatively accurate phase boundaries. Decoupling the many-body wave function into local wave functions $|\Psi_{N}\rangle\approx\Pi_{i}|\Psi_{i}\rangle$, where the $i$-th site wave function is expanded using Fock state basis $|m\rangle$ ($m=0,1,\cdots$) as $|\Psi_{i}\rangle=\sum_{m}f_{m}^{i}(t)|i;m\rangle$ with $|i;m\rangle$ to represent $m$ atoms in $i$-th site and $f_{m}^{i}(t)$ to be the corresponding probability amplitude. The equation of motion of coefficients $f_{m}^{i}$ is

[TABLE]

where $\phi_{i}=\langle\hat{b}_{i}\rangle$ is the SF order parameter.

We find the ground state by solving Eq. (2) in the imaginary time via the Nelder-Mead algorithm. This is done with $32\times 32$ sites and the maximal occupation in each site is 12. In the MI phase, particle numbers at each site are integer and SF order parameter $\phi_{i}=0$. The SF phase is determined by a non-zero, uniform SF order parameter $\phi_{i}=\phi$. The static structure factor is used to identify density-modulated phases,

[TABLE]

where $\mathbf{k}\equiv(k_{x},k_{y})$ is the reciprocal lattice vector (in terms of $1/d$), and $M$ is the total number of lattices Bandyopadhyay et al. (2019). In DW states the structure factor $S(\pi,\pi)\neq 0$ and SF order parameter $\phi_{i}=0$. In SS phases, both SF order parameter and the structure factor are non-zero.

Ground-state phase diagrams for the box interaction have been examined by similar methods (see, e.g. Ref. Iskin and Freericks (2009)). It was shown that the ground state exhibits MI, DW, SS and SF phases by varying hopping, on-site and long-range interactions. We calculate the phase diagram with average particle number in each site $n_{i}=1$. One example for the box interaction with $r_{c}=d$ is shown in Fig. 2a. When the tunnelling and $V_{0}$ are weak, the ground state is a MI. Starting from the MI, the system undergoes a MI-SF phase transition when $J$ increases. With larger $V_{0}$, the ground state enters DW phases when the hopping is small. We find a DW-SS and then SS-SF transition by increasing hopping $J$ Bandyopadhyay et al. (2019). More cases and discussions of the phase diagram can be found in  Iskin and Freericks (2009). In the case of soft-core interactions, the phase diagram has similar distributions of phases. In Fig. 2a, we give one example for $r_{c}=2.5d$. In this case, the SS phase region becomes larger in the parameter space due to the longer-range interaction.

The emergence of the SS phase can be further analyzed using the Bogoliubov theory. In the SF region, the ground state is a homogeneous condensate, while SS leads to spatial modulations in the condensate. When this happens, the Bogoliubov spectrum of the homogeneous SF phase Rey et al. (2003); Macrì et al. (2013),

[TABLE]

becomes unstable. Here \varepsilon(\mathbf{k})=2J\big{[}2-\cos(k_{x}d)-\cos(k_{y}d)\big{]} is the kinetic energy term, $\rho_{0}$ is the condensate density and $\tilde{V}(\mathbf{k})=\frac{1}{M}\sum_{i\neq j}V_{ij}\exp(i\mathbf{k}\cdot\mathbf{r}_{ij})$, is the Fourier transformation of the long-range interaction. In free space, $\tilde{V}(\mathbf{k})$ is calculated analytically, which yields $\tilde{V}(\mathbf{k})=V_{0}\frac{2\pi r_{c}}{|\mathbf{k}|}J_{1}(|\mathbf{k}|r_{c})$ for the box interaction and $\tilde{V}(\mathbf{k})=V_{0}\frac{4\pi r_{c}^{2}}{3}\sum_{m=0,\pm 1}e^{i2m\pi/3}K_{0}(e^{im\pi/3}|\mathbf{k}|r_{c})$ for the soft-core interaction, where $J_{n}(x)$ is the first kind of Bessel function and $K_{n}(x)$ is the second kind of modified Bessel function Hsueh et al. (2012). In the two-dimensional lattice, the Fourier transform of the interaction potential is numerically calculated and one example is presented in Fig. 1b. It is found that $\tilde{V}(\mathbf{k})$ is negative for a large range of $kd$.

When the long-range interaction is strong, the negative components of $\tilde{V}(\mathbf{k})$ cause a roton instability. Here the Bogoliubov spectra $E(\mathbf{k})$ becomes imaginary. The critical hopping (i.e. the lower bound) to observe the roton instability is $J_{l}=-\min[\Sigma(\mathbf{k})]$, where

[TABLE]

A typical $\Sigma(\mathbf{k})$ for the soft-core interaction with $r_{c}=d$ and $V_{0}=0.6U$ is shown in Fig. 2c. When $J=0.5U$, the roton instabilities are triggered and found in the white area in the Bogoliubov spectrum (Fig. 2d). Through numerically fitting $J_{l}$ versus $V_{0}$, its slope for the soft-core interaction becomes larger than that of the box interaction when increasing $r_{c}$ (Fig. 2e). On the other hand, momentum $\mathbf{k}_{\rm rot}$ corresponding to the onset of roton instabilities varies with $r_{c}$. From the numerical data, one finds that $|\mathbf{k}_{\rm rot}|\approx 3\pi/2r_{c}$ when $r_{c}\gg d$ (Fig. 2f).

III Quench Dynamical with box interaction

When $t<0$, the long-range interaction is not present yet. We prepare the system in the MI state with mean particle number $n_{i}=\langle b_{i}^{\dagger}b_{i}\rangle=1$. When $t>0$, the box interaction is turned on instantaneously and the hopping is increased linearly, i.e., $J(t)=a_{Q}t$ with $a_{Q}$ being the quench rate. The time sequence is shown in Fig.1c. Response of the system is captured by average quantities such as superfluid fraction $\rho_{s}=\sum_{\langle i,j\rangle}{\rm Re}[\phi_{i}^{*}\phi_{j}]/zM$ Roth and Burnett (2003); Damski et al. (2003), density variance $\sigma_{n}=\sum_{i}\sqrt{\langle{\hat{n}_{i}}^{2}\rangle-\langle{\hat{n}_{i}}\rangle^{2}}/M$ and vortex nucleation. Here $z=4$ is the coordination number of the 2D square lattice and $M$ is the total number of lattice points. In the MI phase ($J(t)/U\ll 1$), both $\rho_{s}$ and $\sigma_{n}$ vanish.

The dynamics of the eBHM is obtained by numerically solving Eq. (2) using the fourth-order-Runge-Kutta algorithm. In the initial state $f_{1}^{i}\approx e^{i\theta_{i}}$ where $\theta_{i}$ is a random phase that uniformly distributes from [math] to $2\pi$. We then enforce particle number fluctuations in the order of $10^{-3}$ to $m\neq 0$ states and the normalization condition $\sum_{m}|f_{m}^{j}|^{2}=1$. The lattice size is chosen to be $128\times 128$ for box interactions and $48\times 48$ for soft-core interactions with a maximal occupation 7 for each site. The dynamics remains largely the same if we increase the number of sites. We cast 40 trajectories for different initial states and evaluate physical quantities by calculating their average values.

III.1 Superfluid fraction, density variance and vortex density

When the box interaction is weak, the dynamics is largely similar to that of the BHM, as shown in Fig. 3a,b. Initially both $\rho_{s}$ and $\sigma_{n}$ remain small with increasing $J(t)$. During this stage, atomic tunnelling is negligible until $J(t)$ increases to $\sim 0.17U$. This forms the early stage of evolution. At the second stage $J(t)>0.17U$, $\rho_{s}$ and $\sigma_{n}$ increase rapidly. Here the long-range interaction acts as a weak perturbation on top of the strong on-site interaction (see lower panel of Fig. 3a,b, where the difference between interacting and noninteracting case is plotted). The value $J/U\approx 0.17$ is related to the Kibble-Zurek mechanism, which will be discussed in the next section.

Dynamics changes drastically at the two stages when the long-range interaction becomes stronger, as illustrated in Fig. 3c,d. At early times, SF density $\rho_{s}$ is independent of the long-range interaction. However, $\rho_{s}$ decrease with increasing $r_{c}$ at later times. The density variance $\sigma_{n}$ is more sensitive to $r_{c}$. Different curves depart from each other when $J(t)>0.17U$. From a mean-field level, we can understand these curves by noting that the density-density interaction at a given site becomes stronger due to the long-range interaction (i.e., with multiple sites). Eventually the local SF order parameter is reduced with larger soft-core radius (Fig. 3c). At the same time, the density fluctuations are enhanced (Fig. 3d), as density waves are excited during the quench (see examples in Fig. 5).

When $J(t)$ is increased, topological defects (vortices) can be created in the many-body state. Numerically, number densities of vortices $n_{v}$ are evaluated according to $n_{v}=1/(2\pi M)\sum_{i}|\arg(\phi_{i+\hat{x}}^{*}\phi_{i})+\arg(\phi_{i+\hat{x}+\hat{y}}^{*}\phi_{i+\hat{x}})+\arg(\phi_{i+\hat{y}}^{*}\phi_{i+\hat{x}+\hat{y}})+\arg(\phi_{i}^{*}\phi_{i+\hat{y}})|$ with $\hat{x}$ and $\hat{y}$ to be the unit vector along $x$ and $y$ axis. The development of the topological defects in the vicinity of the phase transition point is related to quantities such as the healing length, density fluctuation, etc., as found in BHMs Dziarmaga et al. (2012) and eBHMs with dipolar interactions Yi et al. (2007); Baier et al. (2016). Due to the random phases in the initial state, there are vortices even before the SF order parameter develops. We thus will have phase winding but no currents. To exclude this situation, we define a superfluid vortex density $\rho_{s}n_{v}$, which takes into account contributions from both the SF density and vortex density.

For intermediate quench rate $a_{Q}/U=0.01$, the superfluid vortex density increases rapidly around $J(t)/U=0.17$, as shown in Fig. 4a. In this case, there will be some vortices even when $J(t)$ is in the SF phase region (Fig. 4b). When the quench rate is further decreased to $a_{Q}/U=10^{-3}$, the vortex density increases quickly when $J(t)>0.071U$. This comes from the fact that low-energy modes are excited during the slow quench, which create many vortices. Subsequently, the SF vortex density decreases with increasing $J(t)$. When $J(t)$ is large, vortex and anti-vortex pairs recombine at a higher rate, which speeds up the relaxation of the system to a homogeneous SF with nearly homogeneous phases (Fig. 4c).

III.2 Kibble-Zurek dynamics

The key feature of the KZM in a BHM is that dynamics is divided into frozen and adiabatic region across the phase transition. It is convenient to define a distance parameter $\epsilon=J(t)/J_{s}-1$, which depends on the critical hopping $J_{s}$ of the MI-SF transition. Dynamics is qualitatively different before and after a frozen parameter $\hat{\epsilon}$. As the instantaneous relaxation time $\tau(t)$ diverges around the phase transition, the dynamics is frozen to the initial state if $|\epsilon|<\hat{\epsilon}$. However, dynamical evolution becomes adiabatic when $|\epsilon|>\hat{\epsilon}$, where many dynamical quantities change significantly. The relaxation time can be estimated as $\tau=|\epsilon/\dot{\epsilon}|$. On the other hand, Landau’s mean-field theory predicts that $\tau=\tau_{0}|\epsilon|^{-z\nu}$, where $\nu$ is the critical exponent about correlation length, $\tau_{0}$ is a coefficient of the relaxation time for a BHM, and $z$ is the dynamical exponent. The frozen parameter and frozen (KZM) time are Dziarmaga et al. (2012),

[TABLE]

In practice, it is difficult to determine the KZM time from numerical calculations. In this work, we estimate the time when ${\rm d}^{2}\rho_{s}/{\rm d}t^{2}$ is maximized (indicated by arrows in Fig. 3a,c). In Shimizu, et al’s work the KZM time is chosen to be the time when $|\phi(t_{0}+\hat{t})|=2|\phi(t_{0})|$ Shimizu et al. (2018a). This makes minor differences in $\tau_{0}$, but the scaling exponents remain the same.

To find the universal behaviour in $\rho_{s}$, we change the quench rate from low to high, and find the respective KZM time. The corresponding frozen parameter $\hat{\epsilon}$ is shown in Fig. 5a. It is found that $\hat{t}$ ($\hat{\epsilon}$) is largely independent to types of interactions. The data can be described by a single set of fitting parameters $J_{s}/U=0.0449$, $z\nu=0.442$ and $U\tau_{0}=21.0$. The resulting $J_{s}$ is close to the decoupling approach result $J_{s}/U=0.0429$ van Oosten et al. (2001) (Fig. 2a,b).

At the KZM time, the vortex density $n_{v}\propto\xi^{\sf def}/\xi^{\sf dim}\propto a_{Q}^{({\sf dim}-{\sf def})\nu/(1+z\nu)}$, where $\xi$ is the correlation length, and ${\sf dim}=2$ is the dimension of the lattice and ${\sf def}$ is that of topological defects. Here we argue that ${\sf def}=1$ because the vortices are always created in pairs (e.g., Fig. 4c). From the numerical simulation, the scaling exponent $({\sf dim}-{\sf def})\nu/(1+z\nu)=0.359$ (see Fig. 5b). Additionally, using the fitting results of $\hat{\epsilon}$, we obtain $\nu=0.518$ and $z=0.854$. Here one can make a comparison with mean-field results where $\nu=0.5$, $z=2$, and 3D $XY$ model, which is frequently compared with 2D BHM, where $\nu=2/3$ and $z=1$ Fisher et al. (1989); Elstner and Monien (1999); Hohenadler et al. (2011). In Ref. Shimizu et al. (2018a), it was found that the vortex density has anomalous behaviour in slow quench regime ($a_{Q}/U<10^{-2}$). In our simulations, $n_{v}$ shows larger fluctuation as $a_{Q}$ decreases due to the finite size effect. The average values still agree with the universal scaling.

The breakdown of the universal dynamics is more apparent for slow quenches, where the tunnelling $J(\hat{t})$ at the Kibble-Zurek time $\hat{t}$ is smaller than $J_{l}$. For fast quenches, however, the system will not fully respond to the roton mode, where dynamics is still universal. To demonstrate this, we calculate the density variance $\sigma_{n}$, shown in Fig. 5c-e. The variance reaches a maximal value around $J(t)\sim J_{l}$ when the quench rate is low. In this region, the system suffers significantly from the roton instability, exhibiting apparent density patterns (Fig. 5e3).

IV Time-of-flight analysis

The phases shown above can be probed dynamically through analyzing the momentum distribution and noise correlation Altman et al. (2004); Jeltes et al. (2007); Toth et al. (2008); Hu et al. (2010). By releasing the optical lattice at different times, interference patterns shown in the time-of-flight images encode the phase information. In the deep SF regime, the momentum distribution $n(\mathbf{k})=\langle\hat{b}^{\dagger}_{\mathbf{k}}\hat{b}_{\mathbf{k}}\rangle$, $\hat{b}_{\mathbf{k}}$ being the bosonic operator in momentum space, is approximately given by

[TABLE]

where $\tilde{\phi}_{\mathbf{k}}=1/\sqrt{M}\sum_{i}\phi_{i}\exp[i\mathbf{k}\cdot\mathbf{r}_{i}]$ is the Fourier component of the SF order parameter $\phi_{i}$.

Non-uniform density structures can be characterized e.g. by the noise correlation. This is done by calculating the covariance $C(\mathbf{k},\mathbf{k^{\prime}})=\langle\hat{b}^{\dagger}_{\mathbf{k}}\hat{b}_{\mathbf{k}}\hat{b}^{\dagger}_{\mathbf{k}^{\prime}}\hat{b}_{\mathbf{k}^{\prime}}\rangle-\langle\hat{b}^{\dagger}_{\mathbf{k}}\hat{b}_{\mathbf{k}}\rangle\langle\hat{b}^{\dagger}_{\mathbf{k}^{\prime}}\hat{b}_{\mathbf{k}^{\prime}}\rangle$. The covariance can be obtained through HBT-type interference measurements. In the MI regime, the covariance is simplified to be ,

[TABLE]

where $\tilde{n}_{\mathbf{q}}$ is the Fourier component of occupation $n_{i}$ and the relative momentum $\mathbf{q}=\mathbf{k}^{\prime}-\mathbf{k}$. The analytical expressions for $n(\mathbf{k})$ and $C(\mathbf{k},\mathbf{k}^{\prime})$ for general situations are presented in the Appendix A. Similar results can be found in Refs. Altman et al. (2004); Fölling et al. (2005); Toth et al. (2008).

When system dynamics is frozen, $\tilde{\phi}_{\mathbf{k}}$ and $n(\mathbf{k})$ are almost uniform except for a small peak area centred at $\mathbf{k}=0$ (see Fig. 6). Such a feature is hardly visible in the distribution of $\tilde{n}_{\mathbf{q}}$ and $C(\mathbf{q})$. Note that all these distributions are flat in the initial state. At a later time and for weak long-range interactions, $n(k)$ has a higher peak around $k=0$, signifying the appearance of the SF component. Here only the first Brillouin zone is plotted. When the soft-core radius is small, widths of the momentum distribution decreases slowly with increasing interaction strengths (Fig. 7a1). Here the roton instability is found at a large momentum $\sqrt{2}\pi/r_{c}$ (factor $\sqrt{2}$ is due to the 2D square lattice, see Fig. 8 for details). However, the roton excitation is weak and only perturbs the momentum distribution. When the radius is large, $n(\mathbf{k})$ is affected by rotons (see Fig. 7b3). The appearance of the roton minima causes a flat dispersion relation (Fig. 7a3).

To reveal details of the data shown in Fig. 7, we calculate $\tilde{\phi}_{\mathbf{k}}$ and $\tilde{n}_{\mathbf{q}}$, i.e. Fourier transformation of the order parameter and spatial density. Distributions of these two quantities together with $n(\mathbf{k})$ and $C(\mathbf{q})$ are shown in Fig. 8. In the SF regime, the distributions are more visible in $\tilde{\phi}_{\mathbf{k}}$ and $n(\mathbf{k})$ where a peak is found around $|\mathbf{k}|=0$ (the first row in Fig. 8). For stronger interactions, peaks at different momentum are found (the second and third rows in Fig. 8). These peaks are more profound in the Fourier transformation of the density and covariance. For example, four peaks are found at $|k_{x}|d=|k_{y}|d=\pi$ in the second row. This is because these positions are determined by the soft-core radius Li et al. (2012), given by $\pi d/r_{c}$. When $r_{c}=d$, we thus find the peak positions at $|k_{x}|d=|k_{y}|d=\pi$. Increasing the soft-core radius, more and more peaks emerge. An example with $r_{c}=3d$ is shown in Fig. 8a3-d3. Here the longer-range interactions excite density waves with different characteristic wave-lengths. These extra length scales cause peaks in the momentum distribution.

V Quench dynamics with soft-core interaction

The more realistic soft-core interaction has a similar shape as the box potential at short distances, while decaying quickly with increasing distances. With the same initial state, the dynamics is qualitatively the same as that of the box interaction. This is demonstrated with an example where $r_{c}/d=2.5$. The variables $\rho_{s}$ and $\sigma_{n}$ are shown in Fig. 9a,b. When the soft-core interaction is weak, the dynamics is again similar to the non-interacting case. But as $V_{0}$ grows, $\rho_{s}$ will decrease and $\sigma_{n}$ increases. The rough border between weak and strong interactions can be estimated by the roton instability, which leads to $J\approx 0.29U$. The evolution of the vortex density with varying quench rate is shown in Fig. 9c. Compared to Fig. 4, the notable difference is that decay of $n_{v}$ becomes relatively faster as $a_{Q}$ decreases. In the SF region, the vortex is almost entirely eliminated (Fig. 9d).

Although dynamics between the box and soft-core interaction is similar, the tail of the soft-core interaction alters the time-of-flight results. As shown in Fig. 10a,c, the strength that facilitates density structures is similar to Fig. 7b2,b3, and the peak positions of $C(\mathbf{q})$ also agree with the roton instability analysis. However, the patterns of covariance are apparently different. For the box interactions, $C(\mathbf{q})$ has a clearly octupole geometry. In the soft-core case, more wave numbers are excited in the dynamics such that the peak has less contrast at momentum $|\mathbf{q}|=|\mathbf{k}_{\rm rot}|$ (see Fig. 10d).

VI Conclusion

In this paper, we studied dynamical properties of a two-dimensional eBHM by quenching the atomic hopping. When quenching the hopping from a MI to a SF phase, two stages are found in the dynamical evolution. The dynamics is universal and frozen initially due to the Kibble-Zurek mechanism. After that, the SF order parameter rises quickly with time (hopping strength). For weaker interactions, we observe universal dynamics even after the onset of the SF order parameter. When the long-range interaction is strong, the universal dynamics disappears due to non-negligible SS components. The system eventually enters the SF regime when the hopping is large enough. We also proposed TOF experiments to measure the quench dynamics and determine effects induced by the long-range interaction.

In the future, it is worth exploring situations that are relevant to current cold atom experiments. One could investigate roles played by dimensionality (1D, 3D) and lattice structures (triangular, honeycomb, etc) in the dynamics of the Rydberg dressed gases. Another interesting research topic is to explore the quench dynamics induced by time-dependent long-range interactions. This can be realized experimentally by (adiabatically or abruptly) turning on the Rydberg dressing laser Zeiher et al. (2016). Starting from a SF state, this permits us to understand, e.g. dynamics of quantum depletion in the presence of long-range interactions. Moreover, it has been shown that ground-state phase transitions in spin models can be probed in the quench dynamics Heyl et al. (2013); Titum et al. (2019). An interesting question here is whether one could probe the phase transitions (i.e. SS-SF phase transition) in the extended Bose-Hubbard model through monitoring dynamical quantities in a quench process.

Acknowledgements.

We thank Huaizhi Wu and Mark Fromhold for many useful comments. This work is supported by the UKIERI-UGC Thematic Partnership No. IND/CONT/G/16-17/73, EPSRC Grants No. EP/M014266/1 and No. EP/R04340X/1, Indo-French Centre for the Promotion of Advanced Research, and National Natural Science Foundation of China under Grants No. 11304386 and No. 11774428. The numerical simulation is carried out at the University of Nottingham High Performance Computing Facility. We also acknowledge the support of the University of Nottingham through a RPA grant and China Scholarship Council.

Appendix A Time-of-flight interference and noise correlation

In the time-of-flight experiments, the atomic motion is nearly governed by ballistic expansion. Therefore, the field operator at space-time $\{t;\mathbf{x}\}$ yields $\hat{b}_{\mathbf{k}}=\tilde{w}(\mathbf{k})\sum_{i}\hat{b}_{i}\exp[i\mathbf{k}\cdot\mathbf{r}_{i}]$ with $\mathbf{k}=m\mathbf{x}/(\hbar t)$. Here $\tilde{w}(\mathbf{k})$ is the Fourier component of Wannier function of the lowest band of the optical lattice. We make an assumption that the optical lattice is deep enough such that $|\tilde{w}(\mathbf{k})|$ is approximately a constant. Under Gutzwiller approximation, the expectation value of density operator $n(\mathbf{k})=\langle\hat{n}_{\mathbf{k}}\rangle=\langle\hat{b}^{\dagger}_{\mathbf{k}}\hat{b}_{\mathbf{k}}\rangle$ reads

[TABLE]

Compared to Eq. 7, the second term in Eq. (9) is a constant, which preserves the total number of atoms.

Similarly, the expression for the covariance $C(\mathbf{k},\mathbf{k^{\prime}})=\langle\hat{n}_{\mathbf{k}}\hat{n}_{\mathbf{k^{\prime}}}\rangle-\langle\hat{n}_{\mathbf{k}}\rangle\langle\hat{n}_{\mathbf{k^{\prime}}}\rangle$ is found to be

[TABLE]

where $\mathbf{q}=\mathbf{k}^{\prime}-\mathbf{k}$, $\mathbf{Q}=\mathbf{k}^{\prime}+\mathbf{k}$, and $\eta_{i}=\langle\hat{b}_{i}\hat{b}_{i}\rangle$, $\gamma_{i}=\langle\hat{b}^{\dagger}_{i}\hat{b}_{i}\hat{b}_{i}\rangle$, $\lambda_{i}=\langle\hat{b}^{\dagger}_{i}\hat{b}^{\dagger}_{i}\hat{b}_{i}\hat{b}_{i}\rangle$. $F$ and $\tilde{F}$ represents for the Fourier and inverse Fourier transformation, and $\Delta(\mathbf{q})=1/\sqrt{M}\sum_{i}\exp[i\mathbf{q}\cdot\mathbf{r}_{i}]$, $X(f,g)={\rm Re}[F(f;\mathbf{k})\tilde{F}(g;\mathbf{k})+F(f;\mathbf{k}^{\prime})\tilde{F}(g;\mathbf{k}^{\prime})]$, $Y(f,g,h;\mathbf{k}_{1},\mathbf{k}_{2})={\rm Re}[F(f;\mathbf{k}_{1})F(g;\mathbf{k}_{2})\tilde{F}(h;\mathbf{k}_{1}+\mathbf{k}_{2})]$, $C_{0}=\sum_{i}[8n_{i}|\phi_{i}|^{2}-2n_{i}^{2}-|\eta_{i}|^{2}+\lambda_{i}-6|\phi_{i}|^{4}-4{\rm Re}(\gamma^{*}_{i}\phi_{i}-\eta^{*}\phi_{i}^{2})]$.