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. Quantitative numerical simulations on classical computers are integral to analog quantum simulator experiments. This work studies the propagation of correlations in the out of equilibrium dynamics of the Bose-Hubbard model in two dimensions via tensor-network methods, comparing to predictions in known limits and a recent experiment, also covering the intermediate regime of moderate interactions.

S tate-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][2][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 quantumsimulation 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 propagation 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. 1a). 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][42][43][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 quantumsimulation 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.

Results
Model. We consider the Bose-Hubbard model on a square lattice 53,54 . The Hamiltonian is given aŝ whereâ y i andâ i are the creation and annihilation operators at site i,n i ¼â y iâi is the number operator, J is the strength of the hopping between nearest-neighbor sites, U is the strength of the onsite interaction, and μ is the chemical potential. The notation 〈ij〉 indicates that sites i and j 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 d lat to be unity. The ground state at the commensurate filling is the Mott insulating (superfluid) state for U ≫ J (U ≪ J). Hereafter, we will consider a sudden quench and a quench with a short time (see Fig. 1b and Supplementary Note 1 for details).
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 L x × L y 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 grand potential density hĤi at T = 0 starting from the Mott insulator i n i ¼ 1 should remain constant. They well converge for the bond dimensions D ≥ 6 and remain nearly constant up to t~0.4ℏ/J (see Supplementary Note 2 for the time dependence of the grand potential density). We also investigate how the for −τ Q < t < 0 with τ Q being a finite quench time. In the case of a sudden quench, we discard the region −τ Q < t < 0.
single-particle correlations converge with increasing bond dimensions. The equal-time single-particle correlation function at a distance r = (x, y) for the system size N s is defined as Here ∑ 0 i;j denotes the summation over (i, j) that satisfies |x j − x i | = x and |y j − y i | = y. In the iPEPS simulations, 1=N s ∑ 0 i;j is replaced by 1=2 ∑ i ¼ A;B ∑ 0 j with A and B being sublattice sites because of the translational invariance. As shown in Fig. 2, C sp jrj ¼ 1 ðtÞ :¼ ½C sp r ¼ ð1;0Þ ðtÞ þ C sp r ¼ ð0;1Þ ðtÞ=2 exhibits a peak at t~0.15ℏ/J in both results, and they overlap in this short time. For t ≳ 0.15ℏ/J, the correlation functions of ED start to exhibit a significant finite-size effect, whereas those of iPEPS converge for D ≥ 6. We observe similar behavior for C 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. Figure 3a-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 t~0.12ℏ/J 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 t ≲ 0.4ℏ/J 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 first-peak position of iPEPS with that of the experiment 3 as shown in Fig. 3d. Both iPEPS and experimental results agree very well.
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 v 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 many-body systems, v 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 v 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 ϵ U ( J ðkÞ $ À2J∑ α cos k α (α = x, y 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 single-particle 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 v ¼ max k α ½djϵ U ( J ðkÞj=dk α =_ ¼ 2J=_. Because both doublon and holon quasiparticles propagate with the group velocity v, the front of the correlation function moves at the speed of v front , which should be smaller than v max ¼ 2v ¼ 4J=_. Therefore, this speed v 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 ϵ ðdÞ U ) J ðkÞ $ À4J ∑ α cos k α and ϵ ðhÞ U ) J ðÀkÞ $ À2J ∑ α cos k α , respectively. Because the doublons and holons propagate with respective velocities 4J/ℏ and 2J/ℏ, v front should be smaller than these two sum v max ¼ 6J=_. 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-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 = (x, y) for the system size N s is defined as where 〈⋯〉 c denotes a connected correlation function. In our simulations, hn i ðtÞn j ðtÞi c ¼ hn i ðtÞn j ðtÞi À 1 because hn i ðtÞi ¼ 1 for all sites and time steps. As in C sp r ðtÞ, 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, t ≲ 0.9ℏ/J for U/J~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 t~0.35ℏ/J for |r| = 1, while it appears at t~0.65ℏ/J for |r| = 3. By contrast, the first peak of the density-density correlation appears at t~0.25ℏ/J for |r| = 1, while it appears at t~0.7ℏ/J 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. Supplementary  Notes 8, 9, respectively.) The first-peak time is almost a linear function of the distance, and the system exhibits the light-conelike spreading of correlations (see the time dependence of distance summarized in Supplementary Note 8).

Extraction of propagation velocities in the intermediate and strong interaction regions is summarized in
We summarize the interaction dependence of the group and phase velocities in Fig. 5. In the weak interaction region, the estimated group velocities are v gr~4 J/ℏ. They are similar to those obtained by the TWA at filling factor ν = 10 13 . They are also consistent with the group velocity v gr (U = 0) = 4J/ℏ of a single particle 13 . In the strong interaction region, the estimated group velocity v gr~( 8 ± 2)J/ℏ at U/J = 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 v gr ðU ) JÞ ¼ 6J=_ ½1 þ OðJ 2 =U 2 Þ of a quasiparticle in the large U 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 U c /J~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. ). 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.
Both group and phase velocities gradually converge to the same value as U/J is decreased. This phenomenon can be understood in terms of the separation of the energy scales. When the interaction U is much stronger than the hopping J, 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/U, which determines the phase velocity v ph~U . On the other hand, the time scale of the period of the envelope is~1/J, which determines the group velocity v gr~J . Hence, the group and phase velocities differ as long as U ≫ J. When the interaction U becomes comparable to the hopping J, they start to coincide by slowing down the vibration. Note that this phenomenon occurs irrespective of the presence or absence of phase transitions.

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 tensor-network 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, timeevolved 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.

Methods
Real-time evolution by infinite projected entangled pair states. We prepare iPEPS with a two-site unit cell (see Fig. 1a). The symbols D and D 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 n max as D phys ¼ n max þ 1. Although n max can take infinity in Bose-Hubbard systems, it is practically bounded above in the presence of interaction 62,63 . We can choose finite D phys in the simulations of real-time dynamics. In the case of a sudden quench to the Mott insulating region (U/J > U c /J~16.7 59-61 ), we set the dimension of the local Hilbert space as D phys = 3 because the number of particles deviates only slightly from unity 62,63 . For U/J < U c /J, we choose D phys = 5 so that the wave functions can further take into account the effect of particle fluctuations. When U is close to zero (at U/J = 2 in our simulations), we use slightly larger D 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 i n i ¼ 1 can be represented with the bond dimension D = 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 wave function at each time ΨðtÞ ¼ e ÀitĤ=_ Ψð0Þ is obtained by realtime evolving iPEPS [41][42][43] . The real-time evolution operator in a small time step dt can be approximated by the Suzuki-Trotter decomposition 73-75 as e ÀidtĤ=_ $ Q hiji e ÀidtĤ ij =_ , whereĤ ij ¼ ÀJðâ y iâj þâ y jâi Þ þ U½n i ðn i À 1Þ þ n j ðn j À 1Þ=ð2zÞ À μðn i þn j Þ=z with the coordination number z = 4 is the local Hamiltonian satisfyingĤ ¼ ∑ hijiĤij . After applying the two-site gate e ÀidtĤ ij =_ 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 D. 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 D = 9. Qualitative behavior of correlation functions is found to be nearly the same for D ≥ 6. When extracting the propagation velocities, we mainly use the data for D = 8 and D = 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 D, 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.1 ms 13,14 (see Fig. 1b). For −τ Q < t < 0, both J and U are controlled. The wave function is updated as Ψðt þ dtÞ $ e ÀidtĤðtÞ=_ ΨðtÞ with the time-dependent HamiltonianĤðtÞ in this region. For t > 0, both parameters are fixed. We take J final = J(t = 0)~0.0612ℏ/τ Q as the unit of energy. The discrete time step for the real-time evolution is set to be dt/(ℏ/J final ) = τ Q / (ℏ/J final )/15~0.00408 for all t. To compare the iPEPS results with the exact realtime dynamics in finite-size systems, we also consider a sudden parameter change and set the time step as dt/(ℏ/J) = 0.005. We have checked that the simulations with doubled and halved dt do not change the results significantly.

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. Fig. 5 Propagation velocities as functions of the ratio between the interaction (U) and hopping (J) strengths. The group (v gr , red circles) and phase (v ph , blue squares) velocities are estimated from the density-density and single-particle correlation functions for the bond dimensions D = 8 and D = 9 using the infinite projected entangled pair state (iPEPS) algorithm. The data for D = 8 and D = 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.