Tensor-network study of correlation-spreading dynamics in the two-dimensional Bose-Hubbard model

Recent developments in analog quantum simulators based on cold atoms and trapped ions call for cross-validating the accuracy of quantum-simulation experiments with use of quantitative numerical methods; however, it is particularly challenging for dynamics of systems with more than one spatial dimension. Here we demonstrate that a tensor-network method running on classical computers is useful for this purpose. We specifically analyze real-time dynamics of the two-dimensional Bose-Hubbard model after a sudden quench starting from the Mott insulator by means of the tensor-network method based on infinite projected entangled pair states. Calculated single-particle correlation functions are found to be in good agreement with a recent experiment. By estimating the phase and group velocities from the single-particle and density-density correlation functions, we predict how these velocities vary in the moderate interaction region, which serves as a quantitative benchmark for future experiments and numerical simulations.


I. INTRODUCTION
State-of-art experimental platforms of cold atoms and trapped ions as analog quantum simulators have offered unique opportunities for studying far-from-equilibrium dynamics of isolated quantum many-body systems. Thanks to their high controllability and long coherence time, these platforms have already addressed a variety of intriguing phenomena that are in general difficult to simulate with classical computers, such as correlation spreading 1-3 and relaxation 4-6 after a quantum quench, many-body localization in a disorder potential [7][8][9] , and quantum scar states 10,11 . Nevertheless, accurate numerical methods using classical computers are highly demanded at the current stage of the studies of quantum many-body dynamics, since the classical computation still has complementary advantages over the quantum simulation in that it is free of noise and much more accessible owing to its wide dissemination. In this sense, it is important to cross-check the validity of quantum-simulation experiments and some numerical methods by comparing them with each other.
In particular, direct comparisons between experimental and numerical outputs have been made for dynamical spreading of two-point spatial correlations of the Bose-Hubbard model 1,3,[12][13][14][15] , which can be realized experimentally with ultracold bosons in optical lattices 16 . The correlation spreading has attracted much theoretical interest [12][13][14][15][17][18][19][20][21][22][23][24][25][26][27] in the sense that it is closely related to fundamental phenomena, including the propagation of quantum information and the thermalization. In one spatial dimension, quasi-exact numerical methods based on matrix product states (MPSs) have been used to validate the performance of the quantum simulators 1,3,12 . In two dimensions (2D), by contrast, accurate numerical simulations are challenging. Indeed, the comparisons with respect to a single-particle correlation have shown that a few types of the truncated Wigner approximation (TWA) fail to capture the real-time evolution accurately enough to extract the prop- * rkaneko@phys.kindai.ac.jp † danshita@phys.kindai.ac.jp agation velocity of the correlation 3,14 . Moreover, while the propagation velocities obtained by a two-particle irreducible strong-coupling (2PISC) approach quantitatively agree with the experimental value, and the approach is applicable to much weaker interaction than in the experiment, it does not necessarily provide the exact value of the correlation itself 15 .
In this paper, we present quantitative numerical analyses of the correlation-spreading dynamics of the 2D Bose-Hubbard model starting from a Mott insulating initial state with unit filling. To this end, we employ the tensor-network method based on the infinite projected entangled pair state (iPEPS) [28][29][30][31][32][33][34][35] or the tensor product state [36][37][38][39][40] , which is an extension of MPS to 2D systems [see Fig.1(a)]. The iPEPS studies on real-time dynamics of isolated [41][42][43][44][45][46][47][48][49] and open 41,42,[50][51][52] quantum systems in 2D have begun very recently. Previous simulations suggest that iPEPS can represent relatively low-entangled states in short-time dynamics for simple spin 1/2 systems 41-44 and some itinerant electron systems 46 . This observation may be valid for real-time dynamics in Bose-Hubbard systems; however, little is known about it until now. We find that the single-particle correlation computed with iPEPS, as well as the estimated propagation velocity of the correlation front, agrees very well with the experimental result 3 , demonstrating that iPEPS can be useful for actual quantum-simulation experiments. We also conduct numerical simulations in a moderate interaction region, which has not been addressed by the previous experiments 1, 3 . From the real-time evolution of the single-particle and density-density correlations, we show that the phase and group velocities approach each other when the interaction decreases.

A. Model
We consider the Bose-Hubbard model on a square lattice 53,54 . The Hamiltonian is given aŝ whereˆ † andˆare the creation and annihilation operators at site ,ˆ=ˆ †ˆis the number operator, is the strength of the hopping between nearest-neighbor sites, is the strength of the onsite interaction, and is the chemical potential. The notation indicates that sites and are nearest neighbors. For simplicity, we ignore the effects of the trap potential and the Gaussian envelopes of optical lattice lasers, which do not affect short-time dynamics. We set the lattice spacing lat to be unity. The ground state at the commensurate filling is the Mott insulating (superfluid) state for ( ). Hereafter, we will consider a sudden quench and a quench with a short time [see Fig. 1(b) and Supplementary Note 1 for details].

B. Quench starting from the Mott insulator: Comparison with the exact diagonalization and the experiment
Let us first focus on the case of a sudden quench. We compare our results of iPEPS with those of the exact diagonalization (ED) method and obtain consistent results in a short time. In the ED simulations using the QuSpin library 55,56 , we choose the system sizes × up to 5 × 4 and use the periodic-periodic boundary condition. We examine to what extent the energy is conserved in the iPEPS simulations. The FIG. 2. Single-particle correlation functions sp r ( ) in the case of a sudden quench. Comparison is made between the infinite projected entangled pair state algorithm (iPEPS, blue lines with symbols) and the exact diagonalization method (ED, gray lines). The unit of time is taken as the inverse of the strength of the hopping . is the strength of the interaction.
is the bond dimension. s is the system size. The correlations at distances (a) |r| = 1 and (b) |r| = √ 2 are shown. Both results overlap in a short time.
grand potential density ˆ at = 0 starting from the Mott insulator ⊗ | = 1 should remain constant. They well converge for the bond dimensions ≥ 6 and remain nearly constant up to ∼ 0.4ℏ/ (see Supplementary Note 2 for the time dependence of the grand potential density). We also investigate how the single-particle correlations converge with increasing bond dimensions. The equal-time single-particle correlation function at a distance r = ( , ) for the system size s is defined as Here , denotes the summation over ( , ) that satisfies | − | = and | − | = . In the iPEPS simulations, 1/ s × , is replaced by 1/2 × = , with and being sublattice sites because of the translational invariance. As shown in Next, we compare the correlations of iPEPS with those of the experiment 3 for a finite quench time. In the experiment, a quench to the Mott insulating region has been investigated so far. Figures 3(a-c) show the time evolution of correlations at distances |r| = 1, 2, and 3, respectively. Qualitative behavior is essentially equivalent to the case of the sudden quench, although the correlation function shifts to an earlier time. For |r| = 1, both data show a peak at ∼ 0.12ℏ/ final . Similarly, the first-peak times are consistent with each other for |r| = 2 and 3, and they become longer with increasing distances. When the energy is approximately conserved (namely, for 0.4ℏ/ final , see the time dependence of the grand potential density in Supplementary Note 4), the intensities of correlations also overlap very well. They are also consistent with those obtained by TWA 3,13,14 , while the iPEPS simulations can deal with a slightly longer time and capture the correlation peaks more clearly (see also Supplementary Note 5 for a detailed comparison with the TWA results). To see how well they match more quantitatively, we also compare the firstpeak position of iPEPS with that of the experiment 3 as shown in Fig. 3(d). Both iPEPS and experimental results agree very well.

C. Estimates of group and phase velocities in the moderate interaction region
Having confirmed the applicability of iPEPS simulations to real-time evolution of the Bose-Hubbard model, we study how information propagates by a sudden quench in the moderate interaction region. There are two kinds of velocity that are relevant to the correlation spreading. One is the group velocity gr , which corresponds to the propagation of the envelope of the wave packet and is a suitable quantity to characterize the spreading of correlations. In non-relativistic quantum manybody systems, gr is bounded above, and the upper bound is known as the Lieb-Robinson bound [23][24][25][26][27]57,58 . Notice that the Lieb-Robinson bound for the Bose-Hubbard model has not been rigorously derived with a few exceptions for limited situations 18,[23][24][25][26][27] . The phase velocity ph is the other characteristic quantity, which corresponds to the propagation of the first peak of the wave packet, and does not have to obey the Lieb-Robinson bound.
Although the exact Lieb-Robinson bound is not known for the Bose-Hubbard model, there are some values that can be used as a guide. As discussed in previous studies 1,3,12 in the weak interaction region, the single-particle dispersion up to constant is approximately given as (k) ∼ −2 cos ( = , in 2D), which is equivalent to the dispersion of free particles. The velocity of the correlation spreading would be well characterized by the group velocity of the singleparticle excitation. The largest velocity of a single quasiparticle (along the horizontal or vertical direction ) is described by the maximal slope of the dispersion and is given by = max [ | (k)|/ ]/ℏ = 2 /ℏ. Because both doublon and holon quasiparticles propagate with the group velocity , the front of the correlation function moves at the speed of front , which should be smaller than max = 2 = 4 /ℏ. Therefore, this speed max can be regarded as the Lieb-Robinson-bound-like value. Likewise, in the strong interaction region, the doublon and holon dispersions up to constant are approximately given as ( ) (k) ∼ −4 cos and (ℎ) (−k) ∼ −2 cos , respectively. Because the doublons and holons propagate with respective velocities 4 /ℏ and 2 /ℏ, front should be smaller than these two sum max = 6 /ℏ. Although we know the approximate limit values, the intermediate interaction region is yet to be explored.
To estimate the group velocity from the single-particle correlations, long-time simulations are required in general. However, it is challenging in the iPEPS simulations. To circumvent the difficulty, we estimate the group velocity by the density- . The black symbol corresponds to the first peak in the correlation function obtained by cubic spline interpolation of data points. The propagation velocities along the horizontal or vertical axis are extracted by the data at |r| = 1, 2, and 3. The velocity estimated from the density-density correlation functions is slower than that from the single-particle correlation functions. density correlation. It is known that the propagation velocity of the first peak of this correlation agrees very well with the group velocity 1,12 . The equal-time density-density correlation function at a distance r = ( , ) for the system size s is defined as where · · · c denotes a connected correlation function. In our simulations, ˆ( )ˆ( ) c = ˆ( )ˆ( ) − 1 because ˆ( ) = 1 for all sites and time steps. As in sp r ( ), the summation is replaced by that within sublattice sites in the iPEPS simulations. The parity-parity correlation closely related to the density-density one can be measured in experiments by using the quantum-gas microscope techniques 1 .
We extract the propagation velocities from the first peak in both correlations for |r| = 1, 2, and 3. For simplicity, we consider the sudden quench hereafter. When the interaction becomes weaker, we have confirmed that the energy is conserved in a longer time frame; typically, 0.9ℏ/ for / ∼ 5 (see Supplementary Note 6 for the time dependence of the grand potential density in the weaker interaction region). All the correlation peaks for |r| ≤ 3 appear in this time frame (see Fig. 4). The first peak of the single-particle correlation appears at ∼ 0.35ℏ/ for |r| = 1, while it appears at ∼ 0.65ℏ/ for |r| = 3. By contrast, the first peak FIG. 5. Propagation velocities as functions of the ratio between the interaction ( ) and hopping ( ) strengths. The group ( gr , red circles) and phase ( ph , blue squares) velocities are estimated from the density-density and single-particle correlation functions for the bond dimensions = 8 and = 9 using the infinite projected entangled pair state (iPEPS) algorithm. The data for = 8 and = 9 overlap within the error bars. The velocities and their error bars are obtained by extrapolation of the distance dependence of the peak time. The results obtained by the two-particle irreducible strong-coupling (2PISC) approach 15 (triangles), the truncated Wigner approximation (TWA) 13 (diamonds), and the experiment 3 (a star) are shown. Both velocities gradually merge with decreasing interaction. of the density-density correlation appears at ∼ 0.25ℏ/ for |r| = 1, while it appears at ∼ 0.7ℏ/ for |r| = 3. It takes a long time for propagation in the latter case. (See also the correlations for other interaction parameters given in Supplementary Note 7. Extraction of propagation velocities in the intermediate and strong interaction regions is summarized in Supplementary Notes 8 and 9, respectively.) The first-peak time is almost a linear function of the distance, and the system exhibits the light-cone-like spreading of correlations (see the time dependence of distance summarized in Supplementary Note 8).
We summarize the interaction dependence of the group and phase velocities in Fig. 5. In the weak interaction region, the estimated group velocities are gr ∼ 4 /ℏ. They are similar to those obtained by the TWA at filling factor = 10 13 . They are also consistent with the group velocity gr ( = 0) = 4 /ℏ of a single particle 13 . In the strong interaction region, the estimated group velocity gr ∼ (8 ± 2) /ℏ at / = 19.6 coincides with that obtained by the 2PISC approach 15 within the error bar of extrapolation. It is also comparable to the group velocity ] of a quasi-particle in the large limit 1,3,12 . Similarly, the estimated phase velocity agrees very well with the results of the 2PISC approach 15 and the experiment 3 . In the intermediate region, the estimated group velocity is closer to the single-particle group velocity in the superfluid region, whereas it is comparable to the 2PISC result near and above the critical point c / ∼ 16.7 [59][60][61] .
In all parameter regions, no anomalies appear in the propagation velocities. As for the real-time dynamics after a sudden quench, there is no sign of the superfluid-Mott insulator quantum phase transition. This is because non-universal highenergy excitations come into play during the time evolution. The quantum phase transition at zero temperature does not have to affect the time-evolved states.
Both group and phase velocities gradually converge to the same value as / is decreased. This phenomenon can be understood in terms of the separation of the energy scales. When the interaction is much stronger than the hopping , the correlation function oscillates rapidly as a function of time 12,22 . The correlation function exhibits the envelope of the wave packet. The time scale of the period of oscillation is ∼ 1/ , which determines the phase velocity ph ∼ . On the other hand, the time scale of the period of the envelope is ∼ 1/ , which determines the group velocity gr ∼ . Hence, the group and phase velocities differ as long as . When the interaction becomes comparable to the hopping , they start to coincide by slowing down the vibration. Note that this phenomenon occurs irrespective of the presence or absence of phase transitions.

III. CONCLUSIONS
We have studied real-time dynamics of the 2D Bose-Hubbard model after a sudden quench starting from the Mott insulator with unit filling. We have employed the 2D tensornetwork method based on the iPEPSs, which are the 2D extension of the well-known MPSs in one dimension. Calculated single-particle correlation functions reproduce the recent experimental results very well. The iPEPS algorithm can simulate real-time dynamics long enough for extracting the propagation velocities from correlations. This fact suggests that, for the quench dynamics starting from the Mott insulator in the 2D Bose-Hubbard model, time-evolved states are not so highly entangled before and even slightly after the time at which the correlation front is reached. This finding raises questions about our understanding of how quantum states get entangled with real-time evolution.
We have also estimated the group and phase velocities in the moderate interaction region, in which the 2PISC approach and the TWA are not applicable. The estimated group velocities are continuously connected without singularity in the middle. Our findings would be useful in the future analog quantum simulation and in the future examination of the rigorous Lieb-Robinson bound of Bose-Hubbard systems. The ability of the tensor-network method that accurately calculates the real-time dynamics of 2D quantum many-body systems opens up the possibility of applying it to other quantum-simulation platforms, such as Rydberg atoms, trapped ions, and superconducting circuits.

A. Real-time evolution by infinite projected entangled pair states
We prepare iPEPS with a two-site unit cell [see Fig. 1(a)]. The symbols and phys denote the virtual bond dimension and the dimension of the local Hilbert space, respectively. The former improves the accuracy of the wave function, whereas the latter corresponds to the maximum particle number max as phys = max + 1. Although max can take infinity in Bose-Hubbard systems, it is practically bounded above in the presence of interaction 62,63 . We can choose finite phys in the simulations of real-time dynamics. In the case of a sudden quench to the Mott insulating region ( / > c / ∼ 16.7 59-61 ), we set the dimension of the local Hilbert space as phys = 3 because the number of particles deviates only slightly from unity 62,63 . For / < c / , we choose phys = 5 so that the wave functions can further take into account the effect of particle fluctuations. When is close to zero (at / = 2 in our simulations), we use slightly larger phys = 7 (see Supplementary Note 10 for the details of the choice of the dimensions of the local Hilbert space). The initial Mott insulating state ⊗ | = 1 can be represented with the bond dimension = 1. As for static properties, the Bose-Hubbard model was investigated by finite PEPS or iPEPS, and the phase transition between the Mott insulating and superfluid phases was reproduced [64][65][66][67][68][69][70][71][72] .
The applying the two-site gate −ˆ/ℏ to neighboring tensors, we approximate the local tensors by the singular value decomposition in such a way that the virtual bond dimension of iPEPS remains . In the actual simulations, the second-order Suzuki-Trotter decomposition is used for this simple update algorithm 32,76 , and the tensor-network library TeNeS 77-79 is adopted. The wave functions are optimized up to the bond dimension = 9. Qualitative behavior of correlation functions is found to be nearly the same for ≥ 6. When extracting the propagation velocities, we mainly use the data for = 8 and = 9 to ensure sufficient convergence of physical quantities. We do not preserve the U(1) symmetry during the calculation. Even without respecting the symmetry, we have numerically found that at these values of , the number of particles is nearly conserved during the real-time evolution starting from the Mott insulator.
To compare our results obtained by iPEPS with the experiment 3 , we consider a quench with a short time Q = 0.1ms 13,14 [see Fig. 1(b)]. For − Q < < 0, both and are controlled. The wave function is updated as |Ψ( + ) ∼ −ˆ( )/ℏ |Ψ( ) with the time-dependent Hamiltonianˆ( ) in this region. For > 0, both parameters are fixed. We take final = ( = 0) ∼ 0.0612ℏ/ Q as the unit of energy. The discrete time step for the real-time evolution is set to be /(ℏ/ final ) = Q /(ℏ/ final )/15 ∼ 0.00408 for all . To compare the iPEPS results with the exact real-time dynamics in finite-size systems, we also consider a sudden parameter change and set the time step as /(ℏ/ ) = 0.005. We have checked that the simulations with doubled and halved do not change the results significantly.

V. DATA AVAILABILITY
The data obtained by the iPEPS and ED simulations in this paper are available at https://doi.org/10.5281/zenodo.6085592.
The experimental data 3 and the data obtained by the TWA simulations 13 in this paper are available from the authors upon request.

VI. CODE AVAILABILITY
The codes in this paper are available from the authors upon request.
We briefly review the detailed setup in the experiment [S1] and the time dependence of the hopping and the interaction. In the experiment, ultracold Bose gas of 174 Yb atoms is confined in the optical lattice with the lattice spacing lat = 266nm.
At the time = − Q = −0.1ms, the Mott insulator at = 1 is prepared in the optical lattice with the depth 0 = 15 R . The interaction and the hopping are ( 0 = 15 R ) ∼ 0.648 R and ( 0 = 15 R ) ∼ 0.00652 R , respectively, and their ratio is / ∼ 99.4. Here R /(2 ℏ) = 4021.18Hz is the recoil energy of the system.
For − Q < < 0, the lattice depth is decreased rapidly to 0 = 9 R . The change of the depth is nearly a linear function of time, namely, 0 ( )/ R = 9 − 6 / Q . Finally, the interaction and the hopping become ( 0 = 9 R ) ∼ 0.474 R and ( 0 = 9 R ) ∼ 0.0242 R , respectively. Their ratio becomes / ∼ 19.6, which corresponds to the Mott insulating region close to the critical point c / ∼ 16.7 [S2-S4]. In the iPEPS simulations, we take final = ( 0 = 9 R ) as the unit of energy. The corresponding unit of timescale is unit = ℏ/ final ∼ 1.63ms, and the length of the quench time is expressed as Q = 0.1ms ∼ 0.0612 unit .
For > 0, the hopping and the interaction are fixed. The single-particle correlation functions are obtained by measuring the time-of-flight interference pattern.
All the time dependence of the hopping ( ) and the interaction ( ) are summarized in Fig. S1.

Supplementary Note 2: FURTHER COMPARISON WITH THE EXACT DIAGONALIZATION FOR A SUDDEN QUENCH TO STRONGER INTERACTION
To examine to what extent the energy is conserved in the case of a sudden quench at / = 19.6, we calculate the time dependence of the grand potential density at = 0 obtained by iPEPS in Fig. S2. It is nearly constant ( MI / = −˜/ with˜= / = −0.371 [S2-S4]) for 0.4ℏ/ . The data almost converge for the bond dimensions ≥ 6 in this short time.
We also compare the single-particle correlation functions at a far distance |r| = 2 obtained by iPEPS with those obtained by ED in Fig. S3. In contrast to the ED data at distances |r| = 1 and |r| = √ 2, they do not converge for the system sizes s ≤ 20. This is because a distance |r| = 2 reaches a half the length of the lattice, and consequently, the boundary effect is significant. The data of iPEPS simulations, in which FIG. S1. Time dependence of the depth 0 of the optical lattice, the hopping / R , the interaction / R , and the ratio / . The recoil energy is R /(2 ℏ) = 4021.18Hz. the quantities at the thermodynamic limit are obtained directly, converge very well for ≥ 6. The first peaks in the ED simulations gradually converge to those in the iPEPS simulations as the system size increases. This observation suggests that the iPEPS result represents the correlation functions in the thermodynamic limit fairly well.

Supplementary Note 3: FURTHER COMPARISON WITH THE EXACT DIAGONALIZATION FOR A SUDDEN QUENCH TO WEAKER INTERACTION
To demonstrate the applicability of the iPEPS method for a relatively weaker interaction region, we also compare the single-particle correlation functions obtained by the iPEPS and ED methods. As examples, we show the correlation functions for / = 10 and 4.
For / = 10, the ED results show relatively large size dependence for s ≤ 16 (see Fig. S4). The iPEPS results are well converged for ≥ 7. The data obtained by iPEPS and those by ED for s = 16 overlap very well within a short time frame ( /ℏ 0.3).
On the other hand, for / = 4, the ED results are nearly converged for s ≥ 9 (see Fig. S5). Likewise, the iPEPS results are nearly converged for ≥ 8 when /ℏ 0.6. As in the case of / = 10, the data obtained by iPEPS for = 9 and those by ED overlap very well for /ℏ 0.8.

Supplementary Note 4: FURTHER COMPARISON WITH THE EXPERIMENT FOR A FINITE-TIME QUENCH
As in the sudden quench case, we show the time dependence of the grand potential density at = 0 obtained by iPEPS in the case of a finite-time quench at / final = 19.6 in Fig. S6.  Fig. S7). The intensities along the horizontal direction are slightly smaller than those along the vertical direction in the experiment. Nevertheless, the intensity obtained by iPEPS is almost consistent with that in the experiment. The first-peak times for |r| = 1, 2, and 3 obtained by iPEPS also agree very well with the experimental results.
We also compare the single-particle correlation function at a distance |r| = √ 2 in Fig. S8. The first-peak times are nearly the same between the experiment and iPEPS simulations. However, the intensity obtained by iPEPS is nearly 1.5 times larger than that in the experiment, and the difference is at more than a few-sigma level. Note that the intensity of correlations at a distance |r| = √ 2 by iPEPS overlaps perfectly with those by ED (see the main text) and by TWA (see the next section) within the time shorter than ∼ 0.2ℏ/ final . At the moment, we do not know the origin of the difference between the experimental and theoretical results.

Supplementary Note 5: COMPARISON WITH THE TRUNCATED WIGNER APPROXIMATION FOR A FINITE-TIME QUENCH
We compare the single-particle correlation functions of iPEPS and TWA [S5] in the case of a finite-time quench in  Figure S10 shows the time dependence of the grand potential density at each / . When / ∼ 10, the energy gradually increases with increasing time. This tendency is similar to what we have found at / = 19.6. By contrast, when / ∼ 2, the energy decreases with increasing time. In the intermediate region, especially for / ∼ 5, the behavior of energy going up and going down cancel each other out. Remarkably, it is conserved for 0.9ℏ/ , much longer than the time at / = 19.6. This indicates that the real-time evolution can be simulated in a longer time in the moderate interaction region.

Supplementary Note 7: SINGLE-PARTICLE AND DENSITY-DENSITY CORRELATION FUNCTIONS FOR SEVERAL INTERACTION PARAMETERS
To investigate the propagation velocities, we use the singleparticle and density-density correlation functions. Here we demonstrate how they behave for small and large interaction regions. When the interaction is small ( / = 2), it is difficult to follow the correlation peaks in the single-particle correlation functions within accessible simulation time, as shown in Fig. S11(a). The first peak at |r| = 1 is broad, while peaks do not appear at |r| = √ 2, 2, and 3 for 0.9ℏ/ . This behavior persists until / ∼ 4. Because of the lack of data points of correlation peaks, we can extract the phase velocities only FIG. S10. Time dependence of the grand potential density in the unit of hopping energy for a sudden quench at each / . The energy densities at = 0 are given as dashed lines. When / ∼ 5, the energy density is nearly conserved for a rather longer time 0.9ℏ/ .
for / 5. On the other hand, we observe clear first peaks in the density-density correlation functions [see Fig. S11(b)]. The first-peak times are nearly consistent with those obtained by TWA at high filling = 10 [S6].
When the interaction is strong ( / = 19.6), we success- fully capture all the first peaks up to |r| ≤ 3 within the reliable simulation time 0.4ℏ/ (see Fig. S12). The maximum intensity in the single-particle correlation function (max | sp |r | ( )/ bond | ∼ 0.35) is comparable to that for the small interaction [compare Figs. S11(a) and S12(a)]. By contrast, the maximum intensity of the density-density correlation function (max | dd |r | ( )/ bond | ∼ 0.03) is much smaller than that for the small interaction [compare Figs. S11(b) and S12(b)]. The particle fluctuation is strongly suppressed in the large region, and the peak intensity decays rapidly as a function of a distance.
FIG. S13. Peak-time dependence of distance obtained from the singleparticle correlation functions for = 8 and = 9.
FIG. S14. Peak-time dependence of distance obtained from the density-density correlation functions for = 8 and = 9.

Supplementary Note 8: EXTRACTION OF PROPAGATION VELOCITIES IN THE MODERATE INTERACTION REGION
We show the peak-time dependence of distance obtained from the single-particle correlation functions sp |r | ( ) in Fig. S13. For all the values of / that we have taken in our analyses, the first-peak times of the correlations at |r| = 3 are well below the time when the energy density is no longer conserved (see Fig. S10). At / = 5, the data points for the bond dimensions = 8 and = 9 deviate a little. On the other hand, for / ≥ 10, they are nearly the same and are well converged. The phase velocity along the horizontal or vertical axis is estimated by fitting data points at |r| = 1, 2, and 3. The slopes are found to be nearly the same for = 8 and = 9.
We also show the peak-time dependence of distance obtained from the density-density correlation functions dd |r | ( ) in Fig. S14. In this case, the first-peak times of the correlations at |r| = 3 is comparable to or slightly less than the time when the energy density is no longer conserved (see Fig. S10). Therefore, data points up to |r| = 3 reach near the limit of what we can do in our iPEPS simulations. For all / , the data points for = 8 and = 9 are similar. The group velocity along the horizontal or vertical axis is estimated by fitting data points at |r| = 1, 2, and 3. FIG. S16. Estimated phase velocity from the single-particle correlation functions for the bond dimensions = 8 and = 9.
Note that the peak at |r| = √ 2 does not have to be on the line connecting data points at |r| = 1, 2, and 3. As discussed in Refs. [S1, S7-S9], the correlation spreading in 2D can be highly anisotropic in general. As for the low-energy physics, dynamics of the system is dominated by long wavelength excitations, and the correlation front propagates according to the Euclidean metrics. However, when the system is quenched, various wavelengths are mixed. In this situation, short wavelength excitations also participate, and the Manhattan metrics are more relevant to the propagation of the correlation front. Reflecting this fact, the first peaks of the correlation functions for |r| = 2 and |r| = √ 2, having the same Manhattan distance, appear nearly simultaneously (see Fig. S14).
Estimated group and phase velocities for 2 ≤ / ≤ 22 are summarized in the main text.

Supplementary Note 9: EXTRACTION OF PROPAGATION VELOCITIES FOR MUCH STRONGER INTERACTION
As for much stronger interaction, we successfully estimate the phase velocity using the correlation functions obtained by FIG. S17. Density-density correlation functions per bond at / = 30 that are the same as those in Fig. S15, but the vertical axes are magnified. the iPEPS simulations. On the other hand, it is much harder to extract the group velocity. In this section, we summarize our results for much stronger interaction.
As in the moderate interaction region, we extract the phase velocity from the single-particle correlation functions. They behave similarly between the moderate and stronger interaction regions although each first-peak time of the single-particle correlation shifts to an earlier time as interaction increases (compare, for example, Fig. S12 and Fig. S15). From the first-peak-time dependence of distance, we extract the phase velocity as shown in Fig. S16. The phase velocity is dominated by the energy scale of interaction and is nearly a linear function of / . The estimated value ph ∼ 30 /ℏ at / ∼ 40 is comparable to the value obtained by the 2PISC method in Ref. [S7].
On the other hand, within the accessible simulation time, we are not able to extract the group velocity for very large / . The group and phase velocities start to deviate even when we use the density-density correlation functions to extract the group velocity. To accurately estimate the group velocity, which is slower than the phase velocity, we have to carefully follow the envelope of the correlation function for slightly longer time. To examine how the envelope should look like, we plot the typical behavior of the density-density correlation functions for / = 30 in Fig. S17. When we focus on the correlation function at |r| = 3, the first-peak time is located at /ℏ ∼ 0.2. On the other hand, because this peak | dd |r |=3 ( /ℏ ∼ 0.2)| seems to be smaller than the second peak at /ℏ ∼ 0.5, it is likely that the peak of the true envelope of the correlation function is at the later time than /ℏ ∼ 0.2. However, because the energy is nearly conserved up to a short time /ℏ ∼ 0.4, it is not possible to estimate the peak of the envelope, which should be determined by connecting at least three peak points.