Tensor network simulation of the quantum Kibble-Zurek quench from the Mott to superfluid phase in the two-dimensional Bose-Hubbard model

TL;DR

This paper uses tensor network simulations to study the Kibble-Zurek mechanism in the 2D Bose-Hubbard model, demonstrating that even short quenches can verify predicted scaling laws for correlations and excitations.

Contribution

It provides the first tensor network simulation of the Kibble-Zurek quench in the 2D Bose-Hubbard model, confirming the scaling laws with short ramp times.

Findings

Correlation length and excitation energy follow KZ power laws.

Single particle correlation functions collapse onto a universal scaling function.

Short quenches suffice to test KZ scaling predictions.

Abstract

Quantum simulations of the Bose-Hubbard model (BHM) at commensurate filling can follow spreading of correlations after a sudden quench for times long enough to estimate their propagation velocities. In this work we perform tensor network simulation of the quantum Kibble-Zurek (KZ) ramp from the Mott towards the superfluid phase in the square lattice BHM and demonstrate that even relatively short ramp/quench times allow one to test the power laws predicted by the KZ mechanism (KZM). They can be verified for the correlation length and the excitation energy but the most reliable test is based on the KZM scaling hypothesis for the single particle correlation function: the correlation functions for different quench times evaluated at the same scaled time collapse to the same scaling function of the scaled distance. The scaling of the space and time variables is done according to the KZ power laws.

Tables

Table 1. Table 1: KZ ramp - correlation growth. The speed at which the single particle correlations are spreading in the KZ regime estimated in Fig. 6 for decreasing values of threshold h ℎ h in Eq. ( 16 ). The brackets enclose maximal error bars of the last digit. Its upper speed limit is 6 ​ ( 2 ) ​ J 6 2 𝐽 6(2)J according to Refs. Kaneko and Danshita, 2022 ; Takasu  et al. , 2020 .

h	v
0.3	3.2(7)J
0.2	3.8(8)J
0.1	5.0(11)J
0.05	5.6(17)J

Equations

$$ϵ(t)=\frac{t}{{\tau }_Q}.$$

$$\hat{t}\propto {\tau }_Q^{z\nu /(1+z\nu )}.$$

$$\frac{d\xi }{dt}=\frac{dϵ}{dt}\frac{d\xi }{dϵ}\propto {\tau }_Q^{-1}\frac{1}{|ϵ{|}^{\nu +1}},$$

$$\hat{\xi }\propto {\tau }_Q^{\nu /(1+z\nu )},$$

$$\hat{t}\propto {\hat{\xi }}^z.$$

$$\hat{\xi}^{\Delta_{\cal O}}\mbox{$\langle\psi(t)|$}{\cal O}_{r}\mbox{$|\psi(t)\rangle$}=F_{\cal O}\left(t/\hat{\xi}^{z},R/\hat{\xi}\right),$$

$$\delta =||L-R||,$$

$$\Delta =\sum _i{\delta }_i,$$

$$H$$

$$|11111111....⟩$$

$$C^{sp}_{R}=\frac{1}{2}\mbox{$\langle\psi(t)|$}b_{i}^{\dagger}b_{j}+b_{j}^{\dagger}b_{i}\mbox{$|\psi(t)\rangle$}.$$

$$J=J_c\left[1+ϵ(t)\right],$$

$$\epsilon(t)=\left\{\begin{array}[]{ll}\frac{t}{\tau_{Q}}-\frac{4}{27}\frac{t^{3}}{\tau_{Q}^{3}}&{\rm,when}~{}~{}t<0\\ \frac{t}{\tau_{Q}}&{\rm,when}~{}~{}t\geq 0\end{array}\right.$$

$${\hat{\xi }}^{2{\Delta }_{sp}}{C}_R^{sp}(t=0)=f_C\left(R/\hat{\xi }\right).$$

$${\hat{\xi }}^{2\Delta }{C}_R^{sp}(t)=F_C\left(t/\hat{t},R/\hat{\xi }\right),$$

$$F_C\left(t/\hat{t},R/\hat{\xi }\right)=h$$

$$Q(t)=\lim_{N\to\infty}N^{-1}\left[\mbox{$\langle\psi(t)|$}H(t)\mbox{$|\psi(t)\rangle$}-E_{\rm GS}(t)\right].$$

$${\hat{\xi }}^{z+d}Q(t)=F_Q\left(t/\hat{t}\right),$$

$$T_{\mathrm{KZ}}\propto {\hat{\xi }}^{-z}\propto {\tau }_Q^{-z\nu /(1+z\nu )}$$

Full text

I Quantum Kibble-Zurek mechanism

The Kibble-Zurek mechanism (KZM) originated from a scenario for topological defect formation in cosmological phase transitions driven by expanding and cooling Universe Kibble (1976); *K-b; *K-c. Kibble considered independent selection of broken symmetry vacua in causally disconnected regions. The result is a mosaic of broken symmetry domains, whose size is limited by the causal horizon, leading to topologically nontrivial configurations. However, the speed of light is not relevant for laboratory experiments in condensed matter systems where, instead, a dynamical theory for the continuous phase transitions Zurek (1985); *Z-b; *Z-c; del Campo and Zurek (2014) predicts the scaling of the defects density as a function of the quench rate employing equilibrium critical exponents. It has been verified by numerous simulations Laguna and Zurek (1997); Yates and Zurek (1998); Dziarmaga et al. (1999); Antunes et al. (1999); Bettencourt et al. (2000); Zurek et al. (2000); Uhlmann et al. (2007); *KZnum-h; *KZnum-i; Witkowska et al. (2011); Das et al. (2012); Sonner et al. (2015); Chesler et al. (2015); Liu et al. (2020) and condensed matter experiments Chung et al. (1991); Bowick et al. (1994); Ruutu et al. (1996); Bäuerle et al. (1996); Carmi et al. (2000); Monaco et al. (2002); Maniv et al. (2003); Sadler et al. (2006a); Weiler et al. (2008); Monaco et al. (2009); Golubchik et al. (2010); Chiara et al. (2010); Mielenz et al. (2013); Ulm et al. (2013); Pyka et al. (2013); Chae et al. (2012); Lin et al. (2014); Griffin et al. (2012); Donadello et al. (2014); Deutschländer et al. (2015); Chomaz et al. (2015); Yukalov et al. (2015); Navon et al. (2015); Liu et al. (2018); Rysti et al. (2021). Topological defects play central role in these studies as they survive inevitable dissipation.

Their role was played down in the quantum KZM (QKZM) that considers quenches across quantum critical points in isolated quantum systems Damski (2005); Zurek et al. (2005); Polkovnikov (2005); Dziarmaga (2005, 2010); Polkovnikov et al. (2011); Schützhold et al. (2006); Saito et al. (2007); Mukherjee et al. (2007); Cucchietti et al. (2007); Cincio et al. (2007); Polkovnikov and Gritsev (2008); Sengupta et al. (2008); Sen et al. (2008); Dziarmaga et al. (2008); Damski and Zurek (2010); De Grandi et al. (2010); Pollmann et al. (2010); Damski et al. (2011); Zurek (2013); Sharma et al. (2015); Dutta and Dutta (2017); Jaschke et al. (2017); Białończyk and Damski (2018); del Campo (2018); Puebla et al. (2019); Sinha et al. (2019); Rams et al. (2019); Mathey and Diehl (2020); Białończyk and Damski (2020a); Sadhukhan et al. (2020); Revathy and Divakaran (2020); Rossini and Vicari (2020); Hódsági and Kormos (2020); Białończyk and Damski (2020b); Roychowdhury et al. (2021); Sadhukhan et al. (2020); Schmitt et al. (2022); Nowak and Dziarmaga (2021); Dziarmaga and Rams (2022); Dziarmaga et al. (2022). It was tested by experiments Sadler et al. (2006b); Anquez et al. (2016); Baumann et al. (2011); Clark et al. (2016); Chen et al. (2011); Braun et al. (2015a); Gardas et al. (2018); Meldgin et al. (2016); Keesling et al. (2019); Bando et al. (2020); Weinberg et al. (2020); King et al. (2022a); Semeghini et al. (2021); Satzinger et al. (2021). Recent developments in Rydberg atoms’ quantum simulators  Ebadi et al. (2021); Scholl et al. (2021); Semeghini et al. (2021); Satzinger et al. (2021) and coherent D-Wave King et al. (2022a, b) open possibility to study the QKZM in two and three spatial dimensions and/or to employ it as a test of quantumness of the simulator Nowak and Dziarmaga (2021); King et al. (2022a); Dziarmaga and Rams (2022); Schmitt et al. (2022); Dziarmaga et al. (2022).

The QKZM can be described in brief as follows. A smooth ramp crossing the critical point at time $t=0$ can be linearized in its vicinity as

[TABLE]

Here $\epsilon$ is a dimensionless parameter in a Hamiltonian, that measures distance from the quantum critical point, and $\tau_{Q}$ is called a quench time. Initially, the system is prepared in its ground state far from the critical point. At first the evolution adiabatically follows the ground state of the changing Hamiltonian until the adiabaticity fails near time $-\hat{t}$ when the energy gap becomes comparable to the ramp rate: $\Delta\propto|\epsilon|^{z\nu}\propto|\dot{\epsilon}/\epsilon|=1/|t|$. This KZM timescale is

[TABLE]

Here $z$ and $\nu$ are the dynamical and the correlation length critical exponents, respectively.

From a causality point of view Zurek (1985); *Z-b; *Z-c; Sadhukhan et al. (2020), which is most straightforward when the dynamical exponent $z=1$ and the excitations have a definite speed of sound at the critical point, the correlation length initially grows as $\xi\propto|\epsilon|^{-\nu}$ in step with the correlation length in the adiabatic ground state, that would eventually diverge at the critical point, but near $-\hat{t}$ its diverging growth rate,

[TABLE]

exceeds the speed limit at which correlations can spread near the critical point. The following growth is limited by $2c$ Sadhukhan et al. (2020), where $c$ is the relevant speed of sound at the critical point. The correlation length at $-\hat{t}$,

[TABLE]

defines the characteristic KZ length. Despite the following growth between $-\hat{t}$ and [math], the correlation range when crossing the critical point is also proportional to $\hat{\xi}$ although usually a few times longer Sadhukhan et al. (2020). The causality picture can be generalized to $z\neq 1$ where $c$ has to be replaced by a relevant speed of excitations that depends on $\tau_{Q}$ Sadhukhan et al. (2020).

The two KZ scales are interrelated by

[TABLE]

Accordingly, in the KZM regime after $-\hat{t}$, observables are expected to satisfy the KZM dynamical scaling hypothesis Kolodrubetz et al. (2012); Chandran et al. (2012); Francuz et al. (2016) with $\hat{\xi}$ being the unique scale. For, say, a two-point observable ${\cal O}_{r}$, where $r$ is a distance between the two points, it reads

[TABLE]

where $|\psi(t)\rangle$ is the state during the quench, $\Delta_{\cal O}$ is the scaling dimension, and $F_{\cal O}$ is a non-universal scaling function.

In this paper we consider the QKZM in the 2D Bose-Hubbard model (BHM) on an infinite square lattice. We assume commensurate filling of one particle per site with a well defined Mott-superfluid quantum phase transition. A sudden quench from deep in the Mott phase to the superfluid side of the transition was studied both experimentally Takasu et al. (2020) and numerically Kaneko and Danshita (2022). After the quench the system was allowed to evolve with the final Hamiltonian for a time long enough to estimate the speed at which correlations were spreading – the central phenomenon in the causal interpretation of the QKZM. The aim of the present paper is to demonstrate numerically that these evolution times would also be long enough to verify the KZM scaling hypothesis.

The experimental set up Takasu et al. (2020), where the initial state is a Mott state with the commensurate $n=1$ particle per site, provides an opportunity to go beyond the previous experimental test Chen et al. (2011) where the initial atomic cloud had non-uniform occupation numbers in the range $n=1..3$. $n\approx 3$ in the center of the trap may be just large enough to explain why the measured power laws for relatively fast quenches were consistent with the QKZM but with the mean-field values of the critical exponents. Another attempt was made in Ref. Braun et al., 2015b but a limited range of available parameters made the experimental results inconclusive, though in good agreement with numerical simulations of the experimental set-up. On the numerical front a more tractable 1D version was considered Dziarmaga and Zurek (2014); Gardas et al. (2017) where the Kosterlitz-Thouless nature of the transition makes $\hat{\xi}$ only logarithmically dependent on $\tau_{Q}$ and, therefore, a clear-cut test of the KZM would require quench times ranging over many orders of magnitude. In contrast, the 2D transition is sharper, the KZM power laws are steeper and their experimental verification should be unambiguous. However, numerical simulation of the non-integrable 2D model is more demanding as the applicability of the numerically exact tensor-network DMRG-like methods becomes severely limited in 2D and one may be forced to resort to the mean-field Gutzwiller ansatzShimizu et al. (2018) instead. In this work we overcome the limitations of the quasi-1D DMRG by employing a genuine 2D tensor network.

II 2D tensor network algorithm

Typical quantum many body states can be represented efficiently by tensor networks Verstraete et al. (2008); Orús (2014). These include the matrix product states (MPS) in one dimension (1D) Fannes et al. (1992), the projected entangled pair state (PEPS) in 2D Nishio et al. (2004); Verstraete and Cirac (2004), or the multi-scale entanglement renormalization ansatz (MERA) Vidal (2007, 2008); Evenbly and Vidal (2014a, b) incarnating the real space renormalization group. Recently an infinite PEPS ansatz (iPEPS) was employed to simulate unitary time evolution on infinite latticesCzarnik et al. (2019); Hubig and Cirac (2019); Hubig et al. (2020); Abendschein and Capponi (2008); Kshetrimayum et al. (2020, 2021); Dziarmaga (2021, 2022a); Schmitt et al. (2022); Kaneko and Danshita (2022, 2023). The simulations include spreading of correlations after a sudden quench in the Bose-Hubbard model (BHM)Kaneko and Danshita (2022) and the transverse field Ising model Kaneko and Danshita (2023) as well as the KZ ramp in the latter Schmitt et al. (2022). In this work we perform simulations of the KZ ramp in the BHM that seem timely in view of the new opportunities opened by the recent experiment  Takasu et al. (2020).

We apply the neighbourhood tensor update (NTU) algorithm Dziarmaga (2021) that was previously used to simulate the many body localization Dziarmaga (2022a) and the KZ ramp in the Ising model Schmitt et al. (2022). The evolution operator is Suzuki-Trotter decomposed Trotter (1959); Suzuki (1966, 1976) into a product of nearest neighbor(NN) Trotter gates. As each Trotter gate increases the bond dimension along its NN bond, it has to be truncated back to its original value to prevent its exponential growth with time. The truncation has to be done in a way that minimizes an error afflicted to the quantum state. There are several numerical error measures, each of them implying a different algorithm: the simple update (SU) Hubig et al. (2020); Kshetrimayum et al. (2020), the full update (FU) Phien et al. (2015); Czarnik et al. (2019), the neighbourhood tensor update (NTU) Dziarmaga (2021); Schmitt et al. (2022); Dziarmaga (2022a), or gradient tensor update (GTU) Dziarmaga (2022b). The NTU error measure is explained in Fig. 1. This is the efficient and stable algorithm to be employed here.

In each Trotter gate the Frobenius norm of the difference between the left ($L$) and right ($R$) hand sides of Fig. 1b is minimized. The norm,

[TABLE]

is what we call an NTU error. For small enough time step it should become proportional to $dt$. $\delta$ is an estimate for an error inflicted on local observables by the bond dimension truncation. Accumulating Trotter errors can eventually derail the time evolution. In the worst case scenario the errors are additive. This motivates an integrated NTU error Sinha et al. (2022),

[TABLE]

where the sum is over all performed Trotter gates. For a second order Suzuki-Trotter decomposition on a bipartite square lattice, where each time step is a sequence of $8$ NN Trotter gates, which is 4 gates per site, $4\Delta$ estimates an error of a typical local observable. The observables are calculated with the help of the corner transfer matrix renormalization group Corboz et al. (2014); Corboz (2016).

III Bose Hubbard model

The Hamiltonian on an infinite square lattice is

[TABLE]

Here $b_{i}^{\dagger}$ and $b_{i}$, respectively, creates and annihilates a boson on site $i$, $n_{i}=b^{\dagger}_{i}b_{i}$ is the number operator, $J$ is the strength of the hopping between nearest-neighbor sites, and $U$ is on-site repulsion strength. $\langle i,j\rangle$ denotes summation over nearest-neighbor (NN) pairs in the hopping energy (every pair contributes to the sum only once). For the commensurate filling of $n=1$ particles per site the continuous Mott-superfluid quantum phase transition is located at $U/J=16.7$ Elstner and Monien (1999); Capogrosso-Sansone et al. (2008); Krutitsky (2016). The dynamical exponent $z=1$ and the correlation length exponent $\nu=0.67$, hence $\hat{\xi}\propto\tau_{Q}^{0.40}$.

In an optical lattice both $J$ and $U$ depend on the recoil energy. Deep in the tight binding regime the dependence of $J$ is roughly exponential while that of $U$ is relatively weak (if not negligible). In a tensor network simulation the dimension of the local Hilbert space has to be truncated to a finite physical dimension $d$, i.e., to occupation numbers $0,...,d-1$. This is self-consistent on the Mott side of the transition, including the critical point, thanks to limited variance of occupation numbers $n_{i}$.

IV Sudden quench revisited

As a benchmark, but also to make contact with Ref. Kaneko and Danshita, 2022, we begin with a sudden quench from deep in the Mott insulator phase to the superfluid. In this section we define the energy scale by setting $U=1$. The initial Hamiltonian has zero tunnelling, $J=0$, and the initial ground state is a Fock state:

[TABLE]

with one particle per site. This is a product state that can be represented by an initial iPEPS with bond dimension $1$. Then non-zero tunnelling is suddenly switched on at $t=0$. As in Ref. Kaneko and Danshita, 2022 we consider $1/J=19.6$, i.e., a quench withing the Mott phase. This quench has been performed experimentally in Ref. Takasu et al., 2020 although with a somewhat smoother ramp.

After the quench we follow time evolution of the single particle correlation function

[TABLE]

Here $r$ is a distance between sites $i$ and $j$. Figure 2 shows time evolution of the NN correlator, $C^{sp}_{1}$, up to $Jt=0.5$. Acceptable convergence in this time window requires bond dimension at least $D=11...14$. If we were looking just at $C^{sp}_{1}(Jt)$ then it might appear converged already for $D=6...9$ but closer inspection of the corresponding NTU error in the bottom panel of Fig. 2 reveals that the NTU error does not improve in this range of $D$ as if adding more bond dimension did not improve expressive power of the iPEPS ansatz for this problem. Hidden symmetries may require increasing $D$ not by $1$ but by $2$ or more in order to accommodate not just one more virtual state but a whole multiplet before the expressive power is improved Hasik and Becca (2019). The error begins to improve again from $D=10$ and already $D=11$ brings it down to an acceptable level. At the same time, the curves $C^{sp}_{1}(Jt)$ appear converging again but this time with an acceptable level of the integrated NTU error.

This test shows that a combination of the more $D$-efficient NTU algorithm, than the simple update used in Ref. Kaneko and Danshita, 2022, and higher bond dimensions can significantly increase simulable evolution time. The result encourages us to step beyond the sudden quench and attempt smooth KZ ramps that, by their very nature, take longer times.

V Kibble-Zurek ramp

The Kibble-Zurek quench also begins from the product state (10) but the hopping rate is increased by a smooth ramp instead of the sudden jump. Near the critical point the ramp can be approximated by a linear slope. It is convenient to parameterize the ramp as

[TABLE]

where $J_{c}$ is the critical point and $\epsilon(t)$ is varied from $-1$ to $\infty$ either as a straight linear ramp $\epsilon(t)=t/\tau_{Q}$ or, for instance,

[TABLE]

The former is just linear while the latter can be considered approximately linear in the neighborhood of the critical point at $t=0$, where $\epsilon(t)\approx t/\tau_{Q}$, provided that quench time $\tau_{Q}$ is long enough for $\hat{t}$ in (2) to fall withing the regime of validity of the linearization. The additional qubic term in (13) was added to make its first derivative equal to zero at the beginning of the ramp when $t=-3\tau_{Q}/2$. This smoothing prevents extra initial excitations that would be created by the abrupt beginning of the linear ramp and might overshadow the KZM excitations created near the critical point. They do not pose a problem for long enough $\tau_{Q}$ when their energy, proportional to $\tau_{Q}^{-2}$, becomes negligible compared to the KZM excitation energy that is proportional to $\hat{\xi}^{-3}\propto\tau_{Q}^{-1.2}$, but extra bond dimension would be necessary from the very beginning of the tensor network simulation in order to accommodate their extra entanglement. In principle the extra entanglement is not a problem for a quantum simulator/experiment but the relative suppression of the abrupt excitation still requires longer ramp times that are limited by dissipation. In either case there are good reasons to begin the ramp smoothly.

Furthermore, as the on-site repulsion strength, $U$, depends on the recoil energy relatively weakly — when compared to the hopping rate — here we conveniently assume that it is constant and choose the unit of energy such that $U=1$. Even if we allowed $U$ to be time-dependent it could be linearized near the critical point and the only effect of the time dependence would be effective multiplication of $\tau_{Q}$ by a constant factor. This factor would not affect the KZM scaling hypothesis.

In our simulations the tunneling rate is smoothly ramped up to the critical point at $J_{c}=1/16.7$ with a time step $dt-0.1$ that is short enough for the second order Suzuki-Trotter scheme to be accurate. As our aim is to verify the KZM power laws, quench times are incremented geometrically as $\tau_{Q}=0.1\cdot 2^{m/2}$, where $m$ is a non-negative integer up to $16$. Longer $\tau_{Q}$ require larger bond dimensions, up to $D=14$, as they allow for longer KZM correlation length $\hat{\xi}$ to build up. The accuracy/convergence was monitored with the NTU error as for the sudden quench. We present results obtained with the physical dimension $d=3$. Selective test with $d=4$ show that $d=3$ is accurate enough in consistency with small variance of occupation numbers in our simulations.

Our main focus is the single particle correlation function. It is the most sensitive probe of the KZM as it quantifies just how the long range order builds up when the system is driven across the Mott-superfluid transition. In particular, according to the general KZM scaling hypothesis (6), when the ramp is crossing the critical point at $t=0$ the correlator should satisfy:

[TABLE]

Here $f_{C}$ is a non-universal scaling function, $\Delta$ is an anomalous dimension, and $\hat{\xi}\propto\tau_{Q}^{\nu/(1+z\nu)}$ is the KZ correlation length. The correlator at the critical point is plotted in Fig. 3. The top panel shows raw data for $C^{sp}_{R}(t=0)$ while the bottom one the same data but scaled according to (14). In the rescaling we use $\hat{\xi}=1\cdot\tau_{Q}^{\nu/(1+z\nu)}$ and $\hat{t}=1\cdot\hat{\xi}^{z}$ with the numerical coefficients set equal to $1$ for definiteness. For the single particle correlation function $2\Delta_{sp}=1+\eta$, where $\eta=0.038176(44)$ Elstner and Monien (1999); Capogrosso-Sansone et al. (2008); Krutitsky (2016). The collapse of the plots with different $\tau_{Q}$ demonstrates that we reached quench times long enough for the KZM scaling hypothesis to hold as their $\hat{t}$ is small enough to fall within the critical regime near the transition.

Although the correlation function is not quite exponential, an exponential profile seems to be a reasonably good first approximation that allows to characterize the range of correlations by a single number. In order to ignore numerical noise in the correlator’s long range tail we define the correlation length as $\xi(t)=\ln{C^{sp}_{1}(t)/C^{sp}_{2}(t)}$. The length is plotted in the top panel of Fig. 4 for several different quench times. Furthermore, motivated by a more general KZM scaling hypothesis,

[TABLE]

that should hold in the KZM regime after $-\hat{t}$, in the bottom panel of Fig. 4 we show scaled correlation length, $\xi(t)/\hat{\xi}$, in function of scaled time, $t/\hat{t}$. According to the hypothesis, for long enough $\tau_{Q}$ the scaled plots should collapse in the KZM regime and, indeed, this is what we can see for the slowest quenches. The collapse allows a linear fit to the collapsed sections of the plots after $-\hat{t}$. Our estimate of the slope is $v=1.8(5)J$. According to the causality version of KZM, the slope is upper bounded by twice the sound velocity at the critical point and, indeed, it is lower than the Lieb-Robinson velocity $6(2)J$ predicted and measured in Refs. Kaneko and Danshita, 2022; Takasu et al., 2020, respectively. However, it is strangely low as compared to the upper bound, at odds with many other examples Sadhukhan et al. (2020). We will come back to this issue below.

In the meantime, we observe that the collapse in the bottom panel of Fig. 4 is consistent with the general KZM scaling hypothesis (15). This conclusion is further corroborated by a direct test — without any assumption of an exponential or any other specific profile — made in Fig. 5, where scaled correlation functions for different $\tau_{Q}$, but for the same scaled time $t/\hat{t}=-1$, are plotted together. Their collapse appears even better than the later one at $t/\hat{t}=0$ in Fig. 3. These earlier states are less entangled, their correlations are shorter, hence their representation by the tensor network is more accurate. Similar collapses can be obtained in the whole range $t/\hat{t}\in[-1,0]$ completing demosntration of the KZM scaling hypothesis for the single particle correlation function.

The collapsed correlation functions in Figs. 3 and 5, equal to the scaling functions in (15), provide a more controlled way to estimate the propagation speed Sadhukhan et al. (2020). For a small threshold value $h>0$ and the two values of the scaled time, $t/\hat{t}=-1,0$, equation

[TABLE]

can be solved with respect to scaled distance $R/\hat{\xi}$. Given that for $z=1$ we have $\hat{t}=\hat{\xi}$, the increase of the scaled distance between $t/\hat{t}=-1$ and $t/\hat{t}=0$ is the propagation speed, $v(h)$. Gradually decreasing $h$ allows to probe the speed at which farther correlations are spreading and thus make contact with the Lieb-Robinson bound on the asymptote of the correlation function. Figure 6 shows graphic solution of (16), including its error bars, that results in a series of estimates: $v(0.3)=0.19(4)$, $v(0.2)=0.23(5)$, $v(0.1)=0.30(7)$, $v(0.05)=0.34(11)$ in our units where $U=1$. The same speed estimates for an arbitrary $U$ are listed in Table 1. The speed appears to increase as threshold $h$ is lowered but, at the same time, its error bars increase due to the growing relative significance of numerical uncertainties farther in the correlator tail. Within the error bars the speed is approaching the estimate $6(2)J$ Kaneko and Danshita (2022); Takasu et al. (2020) that is its upper speed limit according to the casual picture of the Kibble-Zurek mechanism.

In addition to the single particle correlation function we can also consider excitation energy per site:

[TABLE]

Here $E_{\rm GS}(t)$ is the ground state energy of the instantaneous Hamiltonian $H(t)$ and $N$ is the number of lattice sites. In the KZ regime after $-\hat{t}$ the excitation energy should satisfy a scaling hypothesis

[TABLE]

where $F_{Q}$ is a non-universal scaling function. On the one hand, with $z+d=3$ the dependence of $Q$ on $\tau_{Q}$ is very steep allowing for a clear-cut test but, on the other hand, with increasing $\tau_{Q}$ the excitation energy quickly becomes a small difference of two large numbers that is prone to numerical errors. Nevertheless, in the top panel of Fig.7 we plot the excitation energy in function of time for several values of the quench time and in the bottom panel we show the same plots but after the rescaling. The scaled plots demonstrate a rather convincing collapse in the KZM regime after $t/\hat{t}=-1$.

VI Thermalization

In order to follow thermalization in the non-integrable model the KZM ramp can be stopped either on the superfluid side of the transition or even right at the critical point where the thermalization should be the most expedient, unhampered by any gap in the energy spectrum. The following unitary evolution with the critical Hamiltonian conserves the KZM excitation energy density $Q\propto\hat{\xi}^{-(z+d)}$ while the state evolves into a thermal one with temperature $T$. The critical dispersion, $\omega\propto k^{z}$, means thermal excitations up to $k_{T}\propto T^{1/z}$ and thermal excitation energy $U_{T}\propto T^{(z+d)/z}$. Equating $Q$ with $U_{T}$ we obtain a “KZ temperature”

[TABLE]

and a thermal correlation range $\xi_{T}\propto k_{T}^{-1}\propto\hat{\xi}$. Despite this proportionality the thermal correlator is not the same as the KZ one immediately after stopping the ramp. Interestingly, similar thermalization at and near the critical point but after a sudden quench was considered in Ref. Haldar et al., 2021.

VII Conclusion

The state of the art quantum simulators of the Bose-Hubbard model at commensurate filling allow one to follow spreading of correlations after a sudden quench for times long enough to estimate their propagation velocities. Our 2D tensor network simulations demonstrate that the experimental times would also be long enough to test the quantum Kibble-Zurek mechanism by verifying the KZM scaling hypothesis for the single particle correlation function. The experiment could push this test beyond the limited range of quench times achievable by the classical simulation where the KZM scaling hypothesis should become even more convincing. It could also follow thermalization of the KZ excitations after the ramp is stopped which is a notoriously difficult task for the classical simulation due to rapid growth of entanglement. These are the challenges worthy a genuine quantum simulation.

Acknowledgements.

We are indebted to Ryui Kaneko for comments on the speed limit for correlations. This research was supported in part by the National Science Centre (NCN), Poland under project 2019/35/B/ST3/01028 (J.M.) and project 2021/03/Y/ST2/00184 within the QuantERA II Programme that has received funding from the European Union Horizon 2020 research and innovation programme under Grant Agreement No 101017733 (J.D.). The research was also supported by a grant from the Priority Research Area DigiWorld under the Strategic Programme Excellence Initiative at Jagiellonian University (J.D.).