## Abstract

An optical lattice quantum simulator is an ideal experimental platform to investigate nonequilibrium dynamics of a quantum many-body system, which is, in general, hard to simulate with classical computers. Here, we use our quantum simulator of the Bose-Hubbard model to study dynamics far from equilibrium after a quantum quench. We successfully confirm the energy conservation law in the one- and three-dimensional systems and extract the propagation velocity of the single-particle correlation in the one- and two-dimensional systems. We corroborate the validity of our quantum simulator through quantitative comparisons between the experiments and the exact numerical calculations in one dimension. In the computationally hard cases of two or three dimensions, by using the quantum-simulation results as references, we examine the performance of a numerical method, namely, the truncated Wigner approximation, revealing its usefulness and limitation. This work constitutes an exemplary case for the usage of analog quantum simulators.

## INTRODUCTION

Rapid advances in analog quantum simulation using highly controllable systems with long coherence time, such as ultracold gases in optical lattices [for example, see (*1*, *2*)], Rydberg atoms in an optical tweezer array [for example, see (*3*–*5*)], and trapped ions [for example, see (*6*, *7*)], have notably expanded possibilities for studying dynamics of quantum many-body systems. One of the recent targets of optical lattice quantum simulators has been the investigation of the nonequilibrium dynamics arising after a quantum quench (*8*–*17*), where a parameter of the system is varied rapidly and substantially. In the case of one dimension (1D) for a short time scale, quantum quench dynamics can be exactly computed with classical computers by means of the matrix product state (MPS) method [for example, see (*18*, *19*)]. In pioneering works of quantum-simulation research, the outputs of experiments were directly compared with those of exact numerical simulations with classical computers to examine the performance of the quantum simulators (*10*–*12*).

A two-point spatial correlation as a function of the distance of the two points has been the intense theoretical interest (*20*–*30*), and in fact, in 1D systems, it has been shown that access to such a correlation function allows for exploring the dynamical spreading of quantum information, which is of great interest in connection with the Lieb-Robinson (LR) bound (*11*, *31*). An exact computation of the spatiotemporal evolution of these two-point correlations is, however, generally intractable for a long time scale or in higher dimensions. While a more recent work has used outputs from a quantum simulator built with ultracold fermions in a Floquet-engineered optical lattice in 3D as a reference for examining the performance of an approximate numerical method, namely, the nonequilibrium dynamical mean-field theory (*32*), a direct comparison with quantitative theoretical approaches in the quench dynamics in higher dimensions is still lacking.

Here, we investigate the energy redistribution dynamics and the spatiotemporal evolution of the single-particle correlation function, which is one of the simplest two-point spatial correlations, in quantum quench dynamics starting with a Mott insulating state by using an optical lattice quantum simulator of the Bose-Hubbard model (BHM) in 2D and 3D as well as 1D. The observation of the redistribution of the kinetic and interaction energies turns out to be the confirmation of the energy conservation in the quench dynamics of a Bose-Hubbard quantum simulator. Furthermore, we successfully observe the correlation spreading after a rapid quench from a Mott insulating state toward the quantum critical region in 2D as well as toward the Mott region in 1D. We compare the measured propagation velocity of the correlation front, which is defined from the first peak in the time evolution of the correlation function at each distance, with the LR-like bound set by the maximum velocity of the quasiparticles. In the 2D case, we find that the former velocity exceeds the latter one. This happens because the single-particle correlation spreads with two typical velocities, namely, the group velocity and the phase velocity, as was pointed out in the recent theoretical work (*30*), and the measured velocity corresponds to the phase one. Since the first peak propagating with the phase velocity decays rapidly with the distance, our observation is not contradicting the existence of the LR bound implying that any correlation functions outside the LR light cone must be exponentially suppressed.

In addition to these experimental findings, to corroborate the quantitative performance of our quantum simulator, we present a thorough comparison between the quantum-simulation results and state-of-the-art theoretical calculations. We use the exact MPS method in the 1D case, finding excellent agreement with the observations. As for the case of the quench toward a deep superfluid region in 3D, the time evolution of the kinetic and interaction energies is directly compared with numerical results obtained using the truncated Wigner approximation (TWA) on the basis of the Gross-Pitaevskii mean-field theory (*29*). The good agreement between the experiment and the theory establishes the predictive power of the TWA for this type of quench. In contrast, in the case of the quench toward the quantum critical region in 2D, the TWA fails to capture quantitatively the experimental results, although it captures some qualitative features. This indicates that our quantum simulation goes beyond current classical computation and the data serve as a useful reference for pushing out its boundary.

## RESULTS

### Investigating nonequilibrium dynamics of the BHM

We consider a system of ultracold bosonic atoms confined in an optical lattice. When an optical lattice potential is deep, the system is quantitatively described by the BHM (*33*, *34*)*j*, *J* is the tunneling-matrix element between nearest-neighbor sites, *U* is the on-site interaction energy, μ is the chemical potential, and *V _{j}* is the local potential offset at the site

*j*, which originates from the trap potential and the Gaussian envelopes of optical lattice lasers. Σ

_{〈j, l〉}represents the summation over all neighboring sites. The position of the site

*j*is denoted by

*x*

^{1},

*x*

^{2}, and

*x*

^{3}mean

*x*,

*y*, and

*z*, respectively.

*e*_{α}represents the unit vector in the

*x*

^{α}-direction, and

*D*is the spatial dimension.

When the atom number per site, namely, the filling factor *U*/*J* is varied, the BHM exhibits a second-order quantum phase transition between the Mott insulator and the superfluid. The system favors the superfluid phase for a relatively small *U*/*J*, while it does the Mott insulator phase for a relatively large *U*/*J*. For the unit filling case (*U*/*J*)* _{c}*= 3.4 (1D), 16.7 (2D), and 29.3 (3D), respectively [for review, see (

*35*)].

Our analog quantum simulator of the BHM is built with an ultracold Bose gas of ^{174}Yb atoms confined in a 3D optical lattice. We use this ^{174}Yb-atom-BHM quantum simulator to analyze dynamics after a quench of the ratio *U*/*J* starting with a Mott insulator state with unit filling. We convert a ^{174}Yb Bose-Einstein condensate (BEC) in a weakly confining harmonic trap into the initial Mott insulator state by slowly ramping up the optical lattice depth up to *s* ≡ *V*_{0}/*E*_{R} = 15 for all the three directions, where *V*_{0} is the depth of the optical lattice and *E*_{R} is the recoil energy of the optical lattice laser whose wavelength is 532 nm. See Methods for the preparation. The prepared state is deep in a Mott insulator regime (*U*/*J* = 100) and is well approximated as a product of local Fock states

To realize a quench of *U*/*J*, we rapidly ramp down the lattice depth for some directions toward a final value. For instance, in the case of the 1D quench, we ramp down the lattice depth only for the *x* direction while in the 3D case, we do it for all the three directions. The ramp-down speed is set to be 100 *E*_{R}/ms. By using the band-mapping techniques, we check that there is no discernible amount of the atoms in excited bands with this quench speed. We use the numerical values of *U* and *J* calculated as functions of lattice depth reported in (*35*) (see section S1).

After the quench process, we keep the lattice depth constant and let the system evolve. To obtain the single-particle correlation function at a certain distance _{x}α ≥ 0, in the unit of a lattice spacing *d* (=266 nm) (*36*)*K*_{Δ}. The kinetic energy of the BHM is equal to the sum of −*JK*_{Δ} at Δ = 1, where Δ = ∣ **Δ**∣ (*36*).

Moreover, we measure the onsite-interaction energy of the BHM, *36*, *37*) for the 3D case and photoassociation spectroscopy (*38*) for the 1D case where we confirm that there are almost no multiple occupancies larger than two (see Methods for details). Compared to the methods based on the quantum gas microscope techniques (*11*), our methods are rather efficient, especially in higher dimensions, for our current purposes of obtaining the ensemble average of the two-point correlation functions and the Hubbard energies, because less repetitions are needed thanks to the much larger number of atoms.

The experimental procedure and setup and typical high-resolution spectra are summarized in sections S1 and S2. It is worth noting that the dynamical evolution of the phase correlation, which is similar to the single-particle correlation, has been measured for weakly interacting Bose gases in 1D optical lattices by means of Talbot interferometry in (*39*). In contrast, the present work investigates the single-particle correlation in strongly correlated regimes in higher dimensions.

### Experimental confirmation of our methods: Dynamics of the 1D BHM after a sudden quench

First, we investigate the behaviors of atoms after a sudden quench in 1D. The results for the dynamical redistribution of the Hubbard energies and the spatiotemporal evolution of the atom correlations are shown in Figs. 1 and 2, respectively. Specifically, we ramp the lattice depth in the *x* direction down to *s* = 5 implying *U*/*J* = 6.8, where the ground state is a Mott insulator state close to the quantum critical point. Figure 1 shows the time evolution of the kinetic energy, the onsite-interaction energy, and the sum of the two energies. On a short time scale, while the sum of the two remains almost constant, the kinetic energy decreases and the interaction energy increases. After making a small overshoot, each energy ends up with an almost steady value, i.e., the energies are redistributed. These behaviors are expected for an isolated system but have never been observed experimentally before.

Figure 2 (A to D) shows the time evolution of the single-particle correlation function for several values of Δ. As the time evolves, the correlations first grow and each of them has the first (local) maximum at a certain time. We extract the peak time for each Δ by numerically fitting to the experimental data, which is plotted against Δ in Fig. 2G. The peak time increases linearly with the distance, i.e., the correlation exhibits a light-cone-like propagation. From the peak time versus Δ, we extract the propagation velocity as *v* = 5.5(7)*Jd*/ℏ.

The maximum velocity of a particle-hole excitation is given as (*11*, *21*)*v*_{max} corresponds to the sum of the maximum velocities of the doublon and the holon, which are respectively given by *J*/*U*. As long as *U*/*J* ≫ 1, Eq. 4 is valid regardless of the spatial dimension. At *U*/*J* = 6.8, *v*_{max} = 5.8*Jd*/ℏ such that the condition *v* < *v*_{max} is satisfied, as expected. A similar propagation behavior has been also observed in the case of the density-density correlation (*11*). In contrast, we will see later that *v* > *v*_{max} in the 2D case. We will explain that this observation is still compatible with the LR-like bound.

While these observations reveal important features of the nonequilibrium dynamics of BHM, this 1D study is also important from another aspect. Since our BHM quantum simulator is analog, it is imperative to examine its accuracy through a direct comparison with exact numerical calculations in 1D before applying it to the cases of higher dimensions, in which exact computation on classical computers is currently unavailable. In Figs. 1 and 2, we compare the experimental results in 1D with the exact numerical ones at zero temperature obtained with the MPS method. For details of the MPS calculations, see section S4. We see that the experimental observations are in good agreement with the exact numerical calculations with no fitting parameters.

### Dynamics of BHM after a sudden quench in higher dimensions

Having corroborated the quantitative validity of our BHM quantum simulator by the comparison between the theory and experiment in 1D, we now discuss the main result of this work, i.e., the quench dynamics in higher dimensions. Figure 3 shows the energy redistribution dynamics for the 3D case after the ramp-down of the lattice depth to *s* = 5 (*U*/*J* = 3.4), where the ground state is deep in the superfluid phase. The general tendency of the time evolution is similar to the 1D case: The two energies are redistributed on a time scale smaller than ℏ/*J* and the sum of the two remains almost constant within the displayed time window *t* ≲ ℏ/*J*.

We next investigate the dynamical spreading of the single-particle correlation after a quantum quench in 2D. The final lattice depth in this case is *s* = 9 implying *U*/*J* = 19.6, where the ground state is a Mott insulator phase near the quantum critical point. Figure 4 (A to D) shows the spatial distribution of the single-particle correlation at several hold times after a quench. We clearly observe that the correlation first grows between nearest-neighbor sites and then it propagates for larger distances at later times. More directly, Fig. 4 (E to G) shows the time evolution of the single-particle correlation function *K*_{Δ}*x*_{,} _{Δ}*y* for several values of (Δ* _{x}*, Δ

*), where Δ*

_{y}*(Δ*

_{x}*) denotes the distance in the*

_{y}*x*(

*y*) direction in units of the lattice spacing

*d*. The delay in the growth of the correlation for longer distance is clearly observed along the directions of

*x*,

*y*, and

*x*+

*y*, in Fig. 4 (E, F, and G, respectively). In the same manner as the 1D case, we extract the position of the first peak in the time evolution of the correlation at each distance, which is plotted against the Euclidean distance

*v*= 13.7(2.1)

*Jd*/ℏ (peak) and

*v*= 10.2(1.4)

*Jd*/ℏ (trough). According to Eq. 3, the maximum velocity of the particle-hole excitation is

*v*

_{max}= 8.4

*Jd*/ℏ, which is slower than the observed propagation velocity in Fig. 4.

## DISCUSSION AND OUTLOOK

The observation that *v* > *v*_{max} in the 2D case shown in Fig. 4 can be interpreted along the line explained in (*30*). The correlation spatially propagates as a wave packet, whose width spreads in time. This means that the velocity of the first peak in the time evolution of the correlation function at each distance, namely, the phase velocity, is faster than that of the center of the wave packet, namely, the group velocity. Moreover, the first peak decays rather rapidly as the distance becomes larger. The velocity extracted from the experimental data in the way described above corresponds to the phase velocity, while the meaningful propagation velocity, which should be compared with the LR-like bound, does to the group velocity. In the case of the final lattice depth *s* = 9 in 2D, we cannot accurately extract the group velocity because of the unclear separation of the two velocities. Instead, In section S5, we show an example, in which the phase velocity is well separated from the group velocity in the 1D case with large *U*/*J*. There, we also see that *v* > *v*_{max}. Hence, the behavior that *v* > *v*_{max} is not unique to the 2D case but can emerge regardless of the spatial dimension as long as *U*/*J* after the quench is sufficiently large. Notice that we observed *v* < *v*_{max} in the case of the final lattice depth *s* = 5 (*U*/*J* = 6.8) in 1D because the phase velocity is approximately equal to the group velocity at *U*/*J* = 6.8 (*30*).

Next, we discuss the usefulness and limitation of some numerical methods based on the quantum simulation results. Since there is no exact computation method applicable to the 2D and 3D cases, it is meaningful to examine the accuracy of some approximate methods by using the quantum simulation results as a quantitative reference. In (*32*), the time-dependent dynamical mean-field theory has been examined by comparison with quantum simulation results for real-time dynamics of the Fermi-Hubbard model. This method is not suited for computing the nonlocal spatial correlations analyzed in the present work because it ignores the momentum dependence of the correlation functions. Instead, we choose the TWA approximation, which is supposed to accurately capture semiclassical dynamics of the BHM at least on a short time scale [see (*29*) and references therein]. In Fig. 3, where the energy redistribution dynamics in 3D is depicted, we also show the numerical calculations as solid lines obtained with the TWA (*29*). In the TWA calculations, we take the Mott insulator state of Eq. 2 as the initial state and set the system size to be 30^{3} sites. We ignore the trapping potential because it is irrelevant to the dynamics within the time window *t* ≲ ℏ/*J* as was discussed in the 1D case. The TWA results are in good agreement with the experimental observations. More details of the TWA calculations are described in (*29*).

Let us turn our attention to the correlation spreading in 2D shown in Fig. 4. The solid lines in Fig. 4 (E to G) represent the results obtained by using the TWA. The TWA agrees with the experiment on a very short time scale (*t* < 0.1ℏ/*J*). Moreover, the peak positions and the values at a relatively long time (*t* > 1ℏ/*J*) for a short distance, say Δ = 1, are reasonably captured. However, it fails to capture some important properties of the correlation dynamics, such as the locations of the correlation troughs and the almost converged value of the correlation for Δ > 1. This disagreement is consistent with the general fact that the TWA is less accurate when *tJ*/ℏ is larger. This failure of the TWA indicates that one needs to push out the boundary of currently available numerical techniques for quantitative description of the quantum simulation results. One possible candidate is to extend the SU(*N*) TWA (*40*, *41*) for analyzing the BHM with unit filling.

In both of the 1D and 2D quench cases, we observed that the peaks propagated linearly with a constant velocity (see Figs. 2G and 4, H and I). However, extrapolations to *t* = 0 have nonzero offsets. In addition, our numerical results also support the existence of the offsets. The offsets reflect the difference between the speed for the creation of a particle-hole pair and that for its propagation. The former speed determines the time giving the first peak at Δ = 1, while the latter does those at Δ > 1. It is noted that the dependence of the propagation velocity on distance in the case of the 1D quench was already numerically discussed in (*11*).

Our quantum simulation platform for studying nonequilibrium dynamics can be straightforwardly applied to other quantum many-body systems such as the Fermi-Hubbard model [with SU(*N*) symmetry (*42*, *43*)], the Bose-Fermi Hubbard model, and the spinful BHM. In addition, it is interesting to extend our work to a study of quench dynamics on a quantum system with controlled dissipation, which has recently attracted much interest (*44*).

## METHODS

### Preparation of the initial Fock state

Details of our experimental setup are described in (*36*). We first prepare a BEC of ^{174}Yb atoms confined in an optical far-off resonant trap (FORT) whose wavelength is 532 nm. The trap frequencies of the FORT are given by (ω_{x′}, ω_{y′}, ω* _{z}*) = 2π × (28,130,160) Hz, where the

*x*

^{′}and

*y*

^{′}axes were tilted from the

*x*and

*y*axes, to which two of the optical lattices are directed, by 45

^{∘}. Then, we slowly ramp up the optical lattice depth for all the three directions from

*s*= 0 to 5 in 100 ms and from 5 to 15 in another 100 ms. A typical number of atoms is chosen to be

*N*= 1.3 × 10

^{4}such that the filling factor is unity.

### Lattice quench

We perform the quench by sudden decrease of the optical lattice with depth of *s E*_{R} in 0.01(15 − *s*) ms. See also section S1. The excitation of the atoms into higher bands is negligible with this procedure. For the cases of the 1D and 2D quench, we decrease the lattice depth along the one direction of *x* and two directions of *x* and *y*, respectively. It is to be noted that when lattice depth is 10.6 *E*_{R}, *U*/*J* is equal to 29.34, which is the critical lattice depth for the superfluid-Mott transition at

### Measurement of the ensemble average of the nonlocal atom correlation

Here, we briefly describe a method for obtaining the ensemble average of the nonlocal atom correlation *K*_{Δ} of Eq. 3. Details are described in (*36*). The atomic-density distribution *n*(**r**) after the TOF *t* is given by*w*_{0}(**r**) and **k** = *m***r**/ℏ*t*.

When *t* is long enough, the structure factor *S*(**k**) is expressed as*K*_{Δ} can be easily obtained by Fourier transformation. For real experiments, two factors should be taken into account: an interaction effect and the finite-TOF effect.

A careful estimation of our experimental conditions (*36*) shows that the ratio of the interaction energy *Un*(*n* − 1)/2 to the kinetic energy ℏω* _{L}* is mostly far lower than 1, justifying our ignorance of the interaction effect during TOF. The finite-TOF effect is small but not negligible so that we determined the nonlocal atom correlation by extrapolation based on the theoretical model described in (

*36*).

### Measurement of ensemble average of the interaction energy

To measure the ensemble average of the interaction energy *36*). First, we increase the optical lattice depth quickly to freeze the hopping of atoms. The ramp-up time is smaller than the hopping time but large enough to prevent the atoms from being excited into the higher band of the optical lattice. For example, the ramp-up time is 0.1 ms from 5*E*_{R} to 15*E*_{R}.

Subsequently, we perform a site-occupancy-resolved spectroscopy. We use two methods: the high-resolution spectroscopy using the optical transition between the ^{1}S_{0} and ^{3}P_{2} (*m _{J}* = 0) electronic states of Yb atoms and the photoassociation spectroscopy. The excellent resolution of the spectroscopy using the

^{1}S

_{0}and

^{3}P

_{2}(

*m*= 0) transition allows us to distinguish different site occupancies, owing to quite different two-body interactions of

_{J}*U*/

_{eg}*h*= −8.5 kHz and

*U*/

_{gg}*h*= 3.2 kHz at 15

*E*

_{R}. From the area of the spectra, we obtain the total number

*N*of

_{n}*n*-occupied sites. The interaction energy is obtained as

*U*= (1/2)Σ

*(*

_{n}N_{n}n*n*− 1). The correlation factors induced by occupancy-dependent Rabi frequencies and a loss of atoms in the

^{3}P

_{2}(

*m*= 0) state during the spectroscopy were studied in our previous work (

_{J}*36*).

Another method that we use is the photoassociation, which is a process to create one molecule from two atoms by light. The created molecule rapidly escapes from the trap so that we can measure the total number of doubly occupied sites as the loss of the atoms. It is noted that the method is invalid in the case of triple and higher occupancies. For example, photoassociation in a triply occupied site induces only two-atom loss and one atom remains, which is the same result as the case of a doubly occupied site, concerning the loss of atoms.

In our experiment, we use the high-resolution spectroscopy using the optical transition between the ^{1}S_{0} and ^{3}P_{2} (*m _{J}* = 0) for the 3D quench experiment. In contrast, we use the photoassociation for the 1D quench experiment, where we additionally check the absence of triply occupied sites or higher by means of the high-resolution spectroscopy.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/40/eaba9255/DC1

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank J. Sakamoto for experimental assistance and T. Kuno, Y. Watanabe, and T. Sagawa for comments and discussions. The MPS-based calculations in Figs. 1 and 2 and fig. S6 were performed using the ITensor Library (http://itensor.org). The TWA simulation in Fig. 4 was carried out at the Yukawa Institute Computer Facility.

**Funding:**This work was supported by the Grant-in-Aid for Scientific Research of the Ministry of Education, Culture Sports, Science, and Technology / Japan Society for the Promotion of Science (MEXT/JSPS KAKENHI) nos. JP25220711, JP16H00801, JP17H06138, JP18H05405, JP18K03492, and JP18H05228; the Impulsing Paradigm Change through Disruptive Technologies (ImPACT) program; Japan Science and Technology Agency CREST (no. JPMJCR1673); and MEXT Quantum Leap Flagship Program (Q-LEAP) grant no. JPMXS0118069021.

**Author contributions:**Y.Takasu and Y.Takahashi conceived and planned the experiments. Y.Takasu, T.Y., H.A., and Y.F. carried out the experiments. Y.Takasu analyzed the experimental data. K.N., S.G., and I.D. developed the theoretical framework and planned and carried out the numerical simulations. Y.Takasu, I.D., and Y.Takahashi wrote the paper with input from all authors.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the corresponding author.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).