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

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-superﬂuid phases. Starting from a Mott insulator state, the dynamical evolution of the superﬂuid order parameter exhibits a universal behaviour at the early stage, largely independent of interactions. The dynamical evolution is signiﬁcantly 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-ﬂight images, i.e. measuring the momentum distribution and noise correlations.

Recently, growing interest has been spent on investigating the dynamics of Bose-Hubbard models (BHMs) driven by external periodic [50,51] or linear fields [52]. A recent review can be found in Refs. [53,54]. A particularly interesting topic is the universal dynamics found in BHMs when linearly changing the hopping strength. Due to the Kibble-Zurek mechanism (KZM) [38][39][40], 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 [55][56][57] leads to a Mott insulator (MI) to superfluid (SF) phase transition at some critical J/U . When the interaction length is greater The tunneling J(t)/U = aQt (orange) is increased linearly during t ∈ (0, t f ). This is done by reducing lattice depth (grey). At the same time, long-range interactions (green dotted) are turned on. When t > t f , atoms are released from the optical lattice to initialize the time-of-flight experiment.
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 [9,[69][70][71][72][73][74]. Rydberg atoms have long lifetimes 10 ∼ 100 µs, growing with n 3 where the principal quantum number n 1. Large polarizabilities (∝ n 7 ) give strong and long-range interactions (e.g. van der Waals interactions ∝ 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 [6,[75][76][77][78]. 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). For distance r < r c , the interaction strength is nearly a constant, while decays quickly if r > r c . Such an interaction potential may be approximated by a box potential, i.e. atomic interaction is a constant when r ≤ 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-neighbor interactions (r c = d) [27,58,59]. As the two types of potential share similar profiles in momentum space (Fig. 1b), we will show that the respective dynamics show common features.
This paper is organized as follows. In Sec. II we introduce the eBHM of the Rydberg dressed atoms. Groundstate 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 [28,68,79] and covariance [57,[80][81][82][83][84] 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 softcore 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, where ij 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 softcore 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 where the Heaviside function Θ(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 [85][86][87][88]. The Gutzwiller approach is a mean field method and predicts qualitatively accurate phase boundaries. Decoupling the many-body wave function into local wave functions |Ψ N ≈ Π i |Ψ i , where the i-th site wave function is expanded using Fock state basis |m (m = 0, 1, · · · ) as |Ψ i = m f i m (t)|i; m with |i; m to represent m atoms in i-th site and f i m (t) to be be the corresponding probability amplitude. The equation of motion of coefficients where φ i = b i 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 × 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 φ i = 0. The SF phase is determined by a non-zero, uniform SF order parameter φ i = φ. The static structure factor is used to identify density-modulated phases, where k ≡ (k x , k y ) is the reciprocal lattice vector (in terms of 1/d), and M is the total number of lattices [89]. In DW states the structure factor S(π, π) = 0 and SF order parameter φ 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. [26]). It was shown that the ground state exhibits MI, DW, SS and SF phases by varying hopping, on-site and longrange 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 [89]. More cases and discussions of the phase diagram can be found in [26]. 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 [90,91], is the first kind of Bessel function and K n (x) is the second kind of modified Bessel function [92]. 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Ṽ (k) is negative for a large range of kd.
When the long-range interaction is strong, the negative components ofṼ (k) cause a roton instability. Here the Bogoliubov spectra E(k) becomes imaginary. The critical hopping (i.e. the lower bound) to observe the roton instability is A typical Σ(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 k rot corresponding to the onset of roton instabilities varies with r c . From the numerical data, one finds that |k rot | ≈ 3π/2r c when r c 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 When the long-range interaction is weak, a MI-SF transition is found by increasing J. Starting from the DW phase (at stronger V0), the ground state undergoes a DW-SS and then a SS-SF transition when increasing J. The dashed line is calculated from the roton instability analysis, which is close to the SS-SF phase boundary. (c) Σ(k) as defined in Eq. (5). The momentum krot is found at the minimal Σ(k) when the spectra become complex. (d) Bogoliubov spectra for J = 0.5U . The white areas indicate roton instability. In (c) and (d) we consider the soft-core interaction with rc = 2d and V0 = 0.6U . (e) Critical tunneling J l for box-type (black solid) and softcore (red dashed) interactions. The interaction lengths are rc = {d, 2d, 3d}, respectively. (f) Momentum |krot| at the roton minimum, in term of its inverse, for box (square) and soft-core (dot) interactions. Here the interaction strength is V0 = 2U .
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 are captured by average quantities such as superfluid fraction [93,94], density variance σ n = i n 2 i − n i 2 /M and vortex nucleation. Here z = 4 is the coordination number of the 2D square lattice, M is the total number of lattice points. In the MI phase (J(t)/U 1), both ρ s and σ n vanish. The dynamics of the eBHM is obtained by numerically solving Eq. (2) using the 4th order Runge-Kutta algorithm. In the initial state f i 1 ≈ e iθi where θ i is a random phase that uniformly distributes from 0 to 2π. We then enforce particle number fluctuations in the order of 10 −3 to m = 0 states and the normalization condition m |f j m | 2 = 1. The lattice size is chosen to be 128 × 128 for box interactions and 48 × 48 for soft-core interactions with a maximal occupation 7. 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.

A. 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 ρ s and σ n remain small with increasing J(t). During this stage, atomic tunnelling is negligible until J(t) increases to ∼ 0.17U . This forms the early stage of evolution. At the second stage J(t) > 0.17U , ρ s and σ 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 non-interacting cases are plotted). The value J/U ≈ 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 ρ s is independent of the long-range interaction. However, ρ s decrease with increasing r c at later times. The density variance σ 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 The development of the topological defects in the vicinity of phase transition point is related to quantities such as the healing length, density fluctuation, etc as found in BHMs [39] and eBHMs with dipolar interactions [25,68]. 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 ρ s n v , which takes into account of 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).

B. 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 = 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 parame-terˆ . As the instantaneous relaxation time τ (t) diverges around the phase transition, i.e. dynamics is frozen to the initial state if | | <ˆ . However, dynamical evolution becomes adiabatic when | | >ˆ , where many dynamical quantities change significantly. The relaxation time can be estimated as τ = | /˙ |. On the other hand, Landau's mean field theory predicts that τ = τ 0 | | −zν , where ν is the critical exponent about correlation length, τ 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 [39], In practice, it is difficult to determine the KZM time from numerical calculations. In this work, we estimate the time when d 2 ρ s /dt 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 |φ(t 0 +t)| = 2|φ(t 0 )| [41]. This makes minor differences in τ 0 , but the scaling exponents remain the same.
To find the universal behaviour in ρ s , we change the quench rate from low to high, and find the respective KZM time. The corresponding frozen parameterˆ is shown in Fig. 5a. It is found thatt (ˆ ) 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ν = 0.442 and U τ 0 = 21.0. The resulting J s is close to the decoupling approach result J s /U = 0.0429 [16] (Fig. 2a,b).
At the KZM time, the vortex density n v ∝ ξ def /ξ dim ∝ , where ξ is the correlation length, and dim = 2 is the dimension of the lattice, def is that of topological defects. Here we argue that def = 1 because the vortices are always created in pairs (e.g., Fig. 4c). From the numerical simulation, the scaling exponent (dim − def)ν/(1 + zν) = 0.359 (see Fig. 5b). Additionally using the fitting results ofˆ , we obtain ν = 0.518, z = 0.854. Here one can make a comparison with mean-field results where ν = 0.5, z = 2, and 3D XY model, which is frequently compared with 2D BHM, where ν = 2/3 and z = 1 [11,95,96]. In Ref. [41], it was found that the vortex density has anomalous behavior 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(t) at the Kibble-Zurek timet 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 σ n , shown in Fig. 5c-e. The variance reaches a maximal value around J(t) ∼ 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 [80,[82][83][84]. By releasing the optical lattice at different times, interference patterns shown in the time-of-flight images encode the phase information. In By integrating the angular direction, we obtain radial distribution of the momentum density n(k) (a1-a4) and covariance C(q) (b1-b4) at different interaction strengths. The probe time is U t = 28 with quench rate aQ/U = 0.01. We consider box interactions with radius rc/d = {1, 2, 3, 4} from left to right, respectively. Dotted lines indicate the momentum where the covariance has maximal values. See Fig. 8 for snapshots of momentum distributions.
the deep SF regime, the momentum distribution n(k) = b † kb k ,b k being the bosonic operator in momentum space, is approximately given by is the Fourier component of the SF order parameter φ i .
Non-uniform density structures can be characterized e.g. by the noise correlation. This is done by calculating the covariance C(k, The covariance can be obtained through HBT-type interference measurements. In the MI regime, the covariance is simplified to be , whereñ q is the Fourier component of occupation n i and the relative momentum q = k − k. The analytical expressions for n(k) and C(k, k ) for general situations are presented in the Appendix A. Similar results can be found in Refs. [57,80,83]. When system dynamics is frozen,φ k and n(k) are almost uniform except for a small peak area centered at k = 0 (see Fig. 6). Such feature is hardly visible in the distribution ofñ q and C(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 √ 2π/r c (factor √ 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(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φ k andñ q , i.e. Fourier transformation of the order parameter and spatial density. Distributions of these two quantities together with n(k) and C(q) are shown in Fig. 8. In the SF regime, the distributions are more visible inφ k and n(k) where a peak is found around |k| = 0 (the 1st row in Fig. 8). For stronger interactions, peaks at different momentum are found (the 2nd and 3rd 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 = π in the second row. This is because these positions are determined by the soft-core radius [28], given by πd/r c . When r c = d, we thus find the peak positions at |k x |d = |k y |d = π. 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 decays 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 ρ s and σ n are shown in Fig. 9a,b. When the soft-core interaction is weak, the dynamics is again similar to the noninteracting case. But as V 0 grows, ρ s will decrease and σ n increases. The rough border between weak and strong interactions can be estimated by the roton instability, which leads to J ≈ 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(q) also agree with the roton instability analysis. However, the patterns of covariance are apparently different. For the box interactions, C(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 |q| = |k rot | (see Fig. 10d). Radial distribution (c) and two dimensional distribution (d) of the covariance C(q). The soft-core radius rc = 3d, quench rate aQ/U = 0.01, probe time U t = 28, and interaction strength V0/U = 0.36. The dotted vertical line marks the momentum corresponding to the roton minimum.

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 [71]. 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 [97,98]. 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.