Quantum phase transition dynamics in the two-dimensional transverse-field Ising model.

The quantum Kibble-Zurek mechanism (QKZM) predicts universal dynamical behavior near the quantum phase transitions (QPTs). It is now well understood for the one-dimensional quantum matter. Higher-dimensional systems, however, remain a challenge, complicated by the fundamentally different character of the associated QPTs and their underlying conformal field theories. In this work, we take the first steps toward theoretical exploration of the QKZM in two dimensions for interacting quantum matter. We study the dynamical crossing of the QPT in the paradigmatic Ising model by a joint effort of modern state-of-the-art numerical methods, including artificial neural networks and tensor networks. As a central result, we quantify universal QKZM behavior close to the QPT. We also note that, upon traversing further into the ferromagnetic regime, deviations from the QKZM prediction appear. We explain the observed behavior by proposing an extended QKZM taking into account spectral information as well as phase ordering. Our work provides a testing platform for higher-dimensional quantum simulators.

INTRODUCTION

The near-critical region of continuous symmetry breaking phase transition is characterized by critical slowing down reflected in an asymptotically divergent relaxation time scale and correlation length. When a many-body system is quenched—driven across the critical point at a fixed rate—critical slowing down will prevent its order parameter from keeping up with what would have been the instantaneous equilibrium. In particular, the correlation length—hence, the size of the fluctuating domains within which the order parameter is approximately uniform—will lag behind the state implied by the externally imposed conditions. Thus, instead of an infinite range order predicted in equilibrium for the posttransition broken symmetry phase, one is left with a mosaic of fluctuating domains. Independent choices of disparate broken symmetry states within each domain result in excitations and lead to the formation of stable topological defects.

Inevitable appearance of the topological defects in the wake of the cosmological phase transitions was pointed out by Kibble (–), who tied their emergence and stability to the homotopy group, and suggested that their posttransition density will be set by thermal activation (which would imply that their density is independent of the rate of quench). Kibble also noted that, in the wake of the Big Bang, the Hubble radius at the epoch of the transition will imply a lower bound on defect density and emphasized that, even at that density, they can have marked consequences for the subsequent evolution of the Universe.

The role of critical slowing down in the creation of topological defects in laboratory phase transition was pointed out by one of us (–). The key difference with the cosmological setting is that now the relativistic causality (reflected in the Hubble radius or the light cone) does not play any useful role. Rather, it is the combination of the critical slowing down and the limits imposed on the growth of the correlation length that decide the size of domains that can independently select how to break symmetry. Hence, one can infer properties of the posttransition broken symmetry state arising from this nonequilibrium process—including the dependence of the density of topological defects and other excitations deposited by the quench on the speed of the transition—from the universal equilibrium scalings and the quench rate (–). This leads to the Kibble-Zurek mechanism (KZM) that is still being tested both numerically and in laboratory experiments (–).

The KZM not only applies to thermal phase transitions but also has been extended to the case of quantum critical points of local Hamiltonians (–). By now, the corresponding quantum KZM (QKZM) has been extensively explored and understood for one-dimensional (1D) models, where, very recently, pioneering experiments in systems of Rydberg atoms have confirmed its predictions in interacting quantum matter (). Other experiments have already started to explore QKZM (–), considered effectively mean-field type systems (, –), or spin models immersed in the environment (, ).

Going beyond 1D for interacting quantum matter appears central not only because of the real-world relevance of 2D systems but also because of fundamental differences compared to 1D. Specifically, conformal field theories describing the universal properties at 2D quantum phase transitions are fundamentally different from their 1D counterparts in that they are generally of a strongly interacting nature. This leads to a change of the character of excitations and the dynamics in the vicinity of the critical point as well as subsequent evolution. While some results on nonequilibrium real-time evolution are available [see (, –)], the theoretical and numerical treatment of interacting quantum 2D systems poses severe challenges. The question to what extent the QKZM also applies to interacting 2D quantum many-body systems has so far remained largely open.

Very recent developments in the Rydberg atom quantum simulator platforms (–) or superconducting qubits () have nevertheless opened the way to access quantum dynamics in 2D at large scales and at a high level of control with the potential to target the outstanding challenge of QKZM in 2D (). A particularly notable step forward is a very recent experiment in a Rydberg atom array (). In the vicinity of the critical point, it is consistent with the QKZM scaling of the correlation length obtained through a fitting procedure to the postquench correlation function extracted from a 16 × 16 square lattice.

In this work, we provide a large-scale numerical analysis of the dynamics in the 2D transverse-field Ising model in systems of various sizes with open and periodic boundary conditions. Going beyond the previous analysis of experimental data, we access the full scaling form of the correlation function, which has the advantage that no prior assumptions on the correlation function are required besides the scaling hypothesis. Thereby, we take the first steps toward a numerical and theoretical study of the dynamics of quantum phase transitions and the QKZM for short-range interacting quantum matter in 2D. In a combined effort, we use a set of state-of-the-art numerical methods. These include time evolution via a time-dependent variational principle (TDVP) for matrix product states (MPS) () or neural quantum states (NQS) (, ) on finite lattices, and 2D infinite projected entangled pair states (iPEPS) operating directly in the thermodynamic limit (, ).

Profiting from the individual strengths of each of the numerical approaches, we find strong evidence for universal dynamical behavior consistent with the scaling properties predicted by the QKZM. We observe that the scaling regime is accessible already at moderate sweep rates and system sizes, in accord with the recent experiments on Rydberg atom quantum simulator platforms. While we can identify the predicted QKZM scaling in the vicinity of the quantum phase transitions with high fidelity, we also observe deviations upon sweeping deeper into the ferromagnetic phase. To account for it, we introduce an extended quantum KZM (xQKZM), which recognizes the effect of the modifications of the system’s spectrum during the sweep and their influence on the excitation energies.

RESULTS

Quantum Kibble-Zurek mechanism

Quantum phase transitions occur in ground states of quantum many-body systems. They mark the singular points where quantum phases of matter transform into each other. In the vicinity of these respective quantum critical points, the macroscopic physical properties become universal as a consequence of a divergent correlation length ξ. QKZM extends this universality—it applies not only to static but also to dynamical properties.

Of paradigmatic importance is the QKZM prediction of universal defect production upon quench through a quantum critical point. Such dynamical crossing can be parametrized by the distance from a quantum phase transition through a dimensionless Hamiltonian parameter ϵ. Close to the critical point, the correlation length ξ in the ground state diverges according to ξ ∝ ∣ϵ∣−ν and the energy gap closes according to ΔE ∝ ξ−, where ν and z are the universal correlation length and dynamical critical exponents, respectively.

To study QKZM, a quantum Ising system is initially prepared in a paramagnetic ground state of its Hamiltonian. It is subsequently smoothly ramped across a quantum critical point to the symmetry-broken state by varying that Hamiltonian. Close to the critical point, the ramp can be linearizedwith t denoting the time. Consequently, t = 0 corresponds to the quantum critical point in the fully adiabatic limit, and the quench time τ sets the speed of the ramp. As long as the evolution is adiabatic, the respective adiabatic correlation length, ξ ∝ ∣ϵ∣−ν, increases at the ratewhich diverges at the critical point. However, perturbations and excitations of the order parameter in a quantum many-body system have a limited maximal speed of propagation c (e.g., the speed of the relevant sound). Therefore, at some point, the actual speed at which correlation length can increase will not be able to keep up. This results in the so-called sonic horizon. It determines the size of the domains that can choose broken symmetry in unison, as shown schematically in Fig. 1.

Fig. 1.

Schematic depiction of the dynamics across a phase transition in a two-dimensional spin-1/2 model.

In the initial paramagnetic state (bottom), spins align with the direction of the transverse magnetic field. A measurement of the spin configuration in that state along the ordering direction would then typically yield a random pattern of spins pointing up (blue cones) or down (red cones). After a slow ramp across a quantum critical point, the system develops a quantum superposition of ferromagnetic domains, which, upon measuring spin configurations along the ordering direction, will yield typically a collapse onto a mosaic of such domains (top). On the front face, we include the growth of the ferromagnetic correlation range as a function of time t starting from t = −τ as the ramp progresses across the critical regime with the critical point located at t = 0. The healing length that determines the size of domains in the Kibble-Zurek (KZ) mechanism is set at the characteristic time , where the growth rate of the instantaneous ground state correlation length ξGS exceeds the maximal speed of the relevant sound, c, in the system.

Consequently, as the critical point is approached, there exists a time where the rate dξ/dt exceeds the relevant sound’s speed c. For z = 1, we have that c = const, whereas for , the velocity becomes scale dependent. For general z, we have that structures of size ξ are subject to the relevant velocity c ∝ ξ−( ∝ ϵν(. Comparing dξ/dt with c results in a characteristic time scaleat which the sonic horizon is determined. The scaling obtained in this way is the same as in the adiabatic-impulse approximation (–, ). The corresponding healing length is set at to be

At large length and time scales, and specify the quench-induced evolution of the system near the critical point.

Setting and methods

Motivated by the recent experiments in Rydberg atom quantum simulators (, ) and by its paradigmatic theoretical relevance, we consider in the following the 2D quantum transverse-field Ising Hamiltonian on a square lattice

Here, , α = x, y, z, denotes the Pauli matrices on lattice site l with the linear extent of the system, L, implying overall L2 lattice sites. The model exhibits a quantum phase transition at g/J = 3.04438 ().

In Fig. 2, we show results for the energy gap between the ground state and the first excited state in the zero momentum and even parity sector as a function of the transverse-field strength g for periodic boundary conditions obtained using NQS and the approach to obtain excited states introduced in (). In the vicinity of the quantum phase transition, as expected, one can observe substantial finite-size effects. Nonetheless, the collapse of the data after finite-size rescaling with the known critical exponents z = 1 and ν = 0.629971, which we show in the inset, reveals consistency with the expected universal behavior.

Fig. 2.

Energy gap as a function of the transverse field g in a periodic lattice.

We collect results for different linear system sizes L and fixed J = J = 1. The data were obtained using the NQS approach for excited states; see the main text and Methods. The inset shows the collapse of the data after finite-size rescaling with the known critical exponents z = 1 and ν ≈ 0.63. The collapse on the paramagnetic side was used to extrapolate the gap to L = 20 in the main panel. The black dashed lines on both sides of the transition—in the regimes near estimated , where the data are also converged in L—are consistent with exponent zν = 0.56 that is shifted by a small nonuniversal 7% correction from the exact 0.63.

In the analysis of the time-evolved system, however, we will later see that the best collapse of the data is achieved assuming zν = 0.56 to be the value of the product of the two critical exponents. The dashed lines in Fig. 2 show power laws of this form fitted to the data on both sides of the transition in regimes where finite-size effects are small. These fits show that the data on these finite intervals are consistent with zν ≈ 0.56. We attribute this to subleading corrections, which, for the ramp times τ that we can numerically achieve, still yield a noticeable contribution, and whose influence can be effectively captured by slightly modified critical exponents.

Of central importance is a quantitative estimate of prefactor in the time scale , whose general scaling form has already been presented in Eq. 3. For what follows, we define as the time at which the rate approximately equals the gap ΔE ∝ (g) (setting the prefactor in that equality to one). With our fitted value for the prefactor of the gap opening on the paramegnetic side, we obtain

The main uncertainty in the fit originates from varying the fitting range. The domain of the fitted curve on the paramagnetic side in Fig. 2 corresponds to the gap at for 0.8 ≤ τ ≤ 6.4. Notice that the value of the prefactor in Eq. 6 depends on the choice of an arbitrary 𝒪(1) prefactor when equating the rate with the gap.

To study the QKZM in our model, we initialize the system in the ground state ∣ψ(t = t)〉 = ∣→ → → … →〉 of the Hamiltonian (5) for J/g = 0, where all spins align along the transverse field. We then numerically solve the Schrödinger equation i∂∣ψ(t)〉 = H(t)∣ψ(t)〉. We fix the unit of time by setting J = 1 (and ℏ = 1). Throughout this work, we study different sweep protocols to ensure that our observations are independent of the protocol details. On the one hand, we consider a linear quenchwith ϵ(t) following Eq. 1, starting at t = −τ with J(t) = 0, crossing the critical point at t = 0, and ending at t = τ. As this sweep exhibits a nonanalytic temporal behavior at the starting t and therefore might potentially generate further excitations masking the targeted QKZ features, we complement our analysis also by a smooth ramp, where the dimensionless distance followsbetween and . It has a vanishing first time derivative at t limiting the generation of additional excitations at the start of the protocol. Both ramps exhibit the same slope in the vicinity of the quantum critical point around t ≈ 0, which, according to the general QKZM argument, is expected to result in identical universal scaling (we show in the Supplementary Materials that those additional excitations are negligible compared to QKZM excitations).

For the simulation of the quantum many-body dynamics, we use three different numerical techniques, allowing us, on the one hand, to perform mutual cross-checks and, on the other hand, to cover complementary regimes of applicability. For a summary of the individual strengths and limitations of the methods, see Table 1. First, we consider iPEPS as a tensor network method that operates directly in the thermodynamic limit of a 2D system. Notably, our iPEPS simulations have been done with the neighborhood tensor update code (), which provides a more stable upgrade of the full update code (). The second method is based on MPS, for which a 1D ordering of the lattice sites is chosen to represent the wave function. Thereby, we can simulate finite systems with open boundary conditions. The MPS wave function is evolved using the TDVP algorithm (). Last, we use neural network quantum states, a recently proposed class of variational wave functions (), to simulate finite systems with periodic boundary conditions. The time evolution was computed using convolutional neural networks (CNNs) and regularization techniques introduced in (). Furthermore, we also use NQS wave functions to variationally obtain the first excited states in addition to the ground state to extract information about the energy gap that we show in Fig. 2.

Table 1.

The used numerical methods have complementary strengths.

The table summarizes their rough ranges of applicability. Here, L is the linear lattice size, is the value of in Eq. 6 with approximated prefactor, and estimates the crossover to adiabatic transition at ; see Fig. 4.

Method	Boundary condition	L	τQ	tf
iPEPS	n.a.	∞	≲3.2	≲2tˆ_actual
NQS	Periodic	≲20	≲τQadiab	≲t^_actual
MPS	Open	≲14	Any	τQ

Fig. 4.

Kibble-Zurek dynamical scaling in 2D quantum Ising model; finite-size scaling.

We show the data collapse corroborating the dynamical scaling hypothesis, combining the results for open (OBC) and periodic boundary conditions (PBC). In (A), we show that the correlation function scaled according to Eq. 11, where finite-size effects are not appreciable in the presented range of . For OBC, we calculate the correlation function with respect to the central spin. In (B), we show the scaled excitation energy density following Eq. 14. For OBC, we compare the total energy (per spin) with the contribution to the excitation energy coming from the central spins. The latter closely follows the PBC results—extending them toward the adiabatic limit, marked by the change of power-law slope. We can estimate the quench rate marking the transition to adiabatic limit , based on crossover happening around . In this plot, we use the exponent with τ ≥ 0.28, as in Fig. 3, that provides the best collapse also in this case.

Universal behavior in the thermodynamic limit

We first consider the ramped quantum Ising model in the putative QKZ regime for times on an infinite lattice in the thermodynamic limit by means of iPEPS simulations. We probe the system’s properties mainly via two quantities. The first one is the ferromagnetic correlation functionwhere R is the distance between the spins m n. Concretely, we compute C(t, R) for spins aligned along one of the axes. Second, we consider the excitation energy densitywhere 〈H(t)〉 = 〈ψ(t)∣H(t)∣ψ(t)〉 denotes the time-dependent expectation value of the Hamiltonian H(t) at time t with ∣ψ(t)〉 being the numerically exact solution of Schrödinger’s equation. Moreover, E0(t) is the ground state energy of H(t) at parameter values g(t) and J(t).

In Fig. 3A, we show the correlation function C(t, R) at t = 0 in units rescaled by the correlation length for various ramp times τ. It is a central result of our work that we observe a data collapse for C(t, R) upon using the known scaling dimension Δ with 2Δ = 1 + η, where η = 0.036298(2) ().

Fig. 3.

Kibble-Zurek dynamical scaling in 2D quantum Ising model; infinite lattice.

In (A), we collect the scaled ferromagnetic correlation function at the critical point, , as a function of the scaled distance, , and in (B), a scaled excitation energy per site, , as a function of the scaled time, (the critical value of the field is reached at t = 0). The main panels show the collapse of the data for slower quenches with τ ≥ 0.28, in agreement with the dynamical scaling hypothesis. We obtain the best collapse for , where for rescaling we fix the prefactor in Eqs. 3 and 4 to one. The exponent that we obtain for the available (limited) τ’s is less than 10% below the expected one of 0.386. Insets show a full range of quench times τ = 0.1 · 2 = 0.1, …,3.2, with integer m in (A) and sparser data with integer m/5 in (B).

Our data directly align with the QKZM prediction implying the following scaling form (–)with F being a nonuniversal scaling function. Overall, this prediction is expected to be asymptotically exact in a coarse-grained sense in the long-wavelength and low-frequency limit, corresponding to long ramp times τ. In Fig. 3A, we find that the scaling regime can already be reached for rather small τ’s, which appears as a promising result from the experimental perspective. Let us emphasize, however, that we obtain the best data collapse for the slower quenches in the range 0.3 < τ < 3.2 for . The value 0.36 for the exponent aligns directly with the exponent identified for using the energy gap; see Eq. 6. As already discussed before, we see a discrepancy of about 7% as compared to the asymptotically expected zν/(1 + zν) = 0.3865, which suggests that subleading corrections beyond the asymptotic universal behavior still yield a weak but noticeable contribution (see the Supplementary Materials for a comparison of the qualities of collapse for the two values of the exponent).

In Fig. 3B, we quantify the number of defects as measured by the excitation energy density Q in the putative scaling regime. In the asymptotic limit of τ → ∞, it is expected from the QKZM in d spatial dimensions that the excitation energy density Q also follows a universal behavior according to (–)

Here, d is a number of space dimensions. In Fig. 3B, we observe a data collapse of the properly scaled Q in a time window in the vicinity of the critical point. As for the ferromagnetic correlation function, the best collapse for the slower quenches is obtained for , also consistent with Eq. 6.

Universal behavior in a finite system

After having explored the case of thermodynamically large systems, in the next step, we now address the experimentally relevant regime of large but finite system sizes with linear extent L. To obtain a comprehensive picture, including the behavior for both open and periodic boundary conditions, we base our analysis on results obtained using MPS and NQS wave functions, where, for MPS (open boundary conditions), we have a smooth ramp in Eq. 8, and for NQS (periodic boundary conditions), we have a linear ramp in Eq. 1. The key consequence of considering a finite system is that the energy gap does not close in the vicinity of the critical point, as illustrated by our numerical results in Fig. 2. This implies that, for sufficiently large τ, an additional adiabatic regime emerges, where the system asymptotically follows the ground state. This finite-size effect can provide a further test of a generalized KZ scaling. On a general level, the respective crossover from QKZ to adiabatic scaling occurs when . While the system follows the QKZ paradigm for , the adiabatic regime is recovered for .

In Fig. 4A, we focus first on the QKZ regime by showing results for the ferromagnetic correlation function C(t, R) for various τ resulting in different .

Here, we again observe a convincing data collapse irrespective of boundary conditions and quench details (smooth or linear), in line with the considerations outlined before, that is expected to yield the QKZ regime. We also note that we identify a collapse of similar quality for times within the predicted scaling regime.

Concerning the excitation energy density, we find power-law behavior for consistent withwhich follows directly from the general QKZM prediction that the excitation energy density is supposed to exhibit the following scaling form (–)

In particular, we emphasize that, for t = 0, our numerical finding of algebraic ξ−3 dependence extends over more than one decade. For larger τ and consequently upon approaching , we also observe the expected deviations toward the adiabatic regime, where . Let us note that the finite-size gap differs for open and periodic boundary conditions, implying that the crossover scale between QKZ and adiabatic is also slightly shifted with respect to each other.

Extended quantum Kibble-Zurek mechanism

After exploring the QKZ scaling in the vicinity of the quantum critical point, we now take the next step by continuing the parameter ramp deep into the ferromagnetic phase down to g = 0. Unlike short-range 1D systems, the 2D quantum Ising model supports a symmetry-broken phase at nonzero temperatures that can lead to other kinds of dynamics such as coarsening or phase-ordering kinetics. These simulations, which correspond to long evolution times, turned out to be numerically the most challenging, and we could not fully converge the iPEPS and NQS simulations for this purpose. Therefore, we rely mostly on the MPS technique in the following.

In Fig. 5A, we display our numerical results for the excitation energy per site Q as a function of τ in rescaled units for finite system sizes obtained using MPS. While we observe a kind of data collapse, as one might expect from the general picture of the QKZM, there does not appear a clear power-law behavior as the collapsed data exhibit a slight bending in the used double-logarithmic plot (apart from the expected crossover to the exponential scaling in the adiabatic limit at large τ).

Fig. 5.

Excitation energy measured deep in the ferromagnetic phase.

In (A), we plot the excitation energy per site at the final g = 0 for a system with OBC. The data for different system sizes are plotted as a function of , and the energy is scaled consistently with the relation that follows from Eq. 15. This results in a relatively good—but imperfect—collapse for available system sizes. We show both the total energy (per spin) and the contribution coming from the central spins. In (B), we test the xQKZM conjecture of adiabatic excitation energy rescaling in Eq. 15, plotting the ratio of the excitation energy to the energy gap at different values of ϵ, comparing it with the final ϵ = 1. The flattening of the curves in the central part of the panel supports the conjecture. Here, we use L = 14 with OBC. Last, in (C), we quantify the deviation from the pure power-law dependence on τ that follows from the simple conjecture in Eq. 15. Snapshots of excitation energies at other values of g show its buildup with an increasing energy gap (decreasing g). Here, we calculate energy contribution for central spins in a system with OBC and L = 14, which, in the KZ regime, we also compare with the available iPEPS results (crosses). We fit Eq. 16 with parameters (a, b) equal to (1.4, 0.16) for g = 0 (ϵ = 1), (1.2, 0.2) for ϵ = 0.5, and (1.0, 0.3) for ϵ = 0.3. The fit is restricted to the regime of validity of Eq. 16 in between the KZ regime and the adiabatic regime. We use calculated at twice the estimate in Eq. 6 that is also used to single out part of dynamics within the universal KZ regime (red filled area).

To understand and quantify the deviations from the expected power-law dependence, we now formulate a simple conjecture to predict the excitation energy at the end of the ramp at g = 0 far beyond the universal regime, which terminates around . We call this conjecture the xQKZM, generalizing a concept established in 1D () to our interacting 2D setting. After , the evolution becomes adiabatic; i.e., the gap becomes large enough to prevent any transfer of occupation between the instantaneous ground state and excited states. The xQKZM is based on two assumptions, whose validity and limitations for the considered parameter regimes will be discussed later on: (i) after , the redistribution of occupations among the instantaneous excited eigenstates can be neglected, so that energy changes can only emerge from the parametric dependence of the energy eigenvalues corresponding to the eigenstates during the sweep and (ii) the details of the parametric dependence of the relevant energy eigenvalues can also be neglected so that the eigenvalues exhibit roughly a global rescaling by a scale set by the gap ΔE(ϵ). A schematic depiction displaying the parametric dependence of the gap and the occupations is shown in Fig. 6.

Fig. 6.

Schematic depiction of the xQKZM.

The lines in the plot show the energies of the ground state and of the first excited state (combined data from Fig. 2); the shaded area represents further excited states. The circles indicate the occupation of the ground state and of the low-lying excited states. The assumption of the xQKZM is that, for a given rate τ, excited state occupation is generated during the time interval . Subsequently, the occupations are assumed to be time-independent, but excitation energy is increased because of adiabatic gap rescaling.

This yields the following estimate for the excitation energy Q at the end of the ramp

In the first step, we assume that the excitation energy at , corresponding to the time , is proportional to in accordance with the scaling hypothesis in Eq. 12 and the data collapse in Fig. 3B. In the second step, we assume a further simplification in that the gap scales as over the full range of considered ; this assumption holds for small enough (or, equivalently, slow enough quenches) and large enough system sizes such that the slowest quenches are still not adiabatic. In this regime, the simple conjecture predicts that, in particular, for g = 0, the excitation energy should scale as . This power law is indicated by dashed lines in Fig. 5 (A and C). While it captures the leading trend, our data show a deviation beyond the leading power-law behavior that we discuss later. An equation similar to Eq. 15 appeared in () but for small ϵ where ΔE(ϵ) ∝ ϵ.

In Fig. 5B, we test the xQKZM prediction including the ratio ΔE(ϵ = 1)/ΔE(ϵ) numerically by comparing the relative excitation energy change Q(ϵ = 1)/Q(ϵ) to the relative change ΔE(ϵ = 1)/ΔE(ϵ) of the gap between some intermediate values ϵ of the ramp and the end ϵ = 1. In accordance with Eq. 15, we can see flattening of the curves in the central part of Fig. 5B that lies between the KZ regime on the left (where ) and the adiabatic regime on the right.

Let us now discuss the regime of validity of the xQKZM conjecture. Concerning the first assumption of neglecting the redistribution of occupations among eigenstates, this clearly depends on the quench times. For sufficiently large τ, the general expectation would be that the system is supposed to undergo phase-ordering kinetics and coarsening dynamics, which would involve a second type of universal dynamical process and which originates precisely from redistribution of occupations. However, for the τ considered in the numerical computations of this work, we are operating in a different regime. While the τ values are still sufficiently large to observe numerical evidence for the QKZM, the overall time span of the dynamics in the ferromagnetic phase for times is limited so that such redistribution, for now, can be approximately neglected. Concerning the second assumption of a roughly uniform shift of the relevant energy eigenvalues deep in the ferromagnetic phase, the respective validity depends crucially on whether the dominant occupation of the instantaneous eigenstates originates from states of the order of the instantaneous gap or, more specifically, from eigenstates not too far up in the excitation spectrum. Although the QKZM describes the creation of occupations in excited states, it is still reasonable to assume that these excitations are not predominanltly located in the nonuniversal high-energy regime, where universality would anyway be out of reach. Furthermore, by increasing τ, the likelihood of generating high-energy excitations can be systematically decreased, so that assumption (ii) is more likely to be approximately valid.

In Fig. 5C, we test the conjecture in Eq. 15 against our numerical data. We find that the bare xQKZM already accounts for the main contributions to the excitation energies. Consistent with Fig. 5B, however, we also observe that a quantitative comparison requires further corrections. Empirically, we find that these are consistent with a logarithmic dependence on τwith a and b being constants. In Fig. 5C, we include the xQKZM in combination with these logarithmic corrections to the numerical MPS data and observe a close correspondence. The main influence of the logarithmic correction is to impose a bending of the excitation energy toward smaller τ. Let us note that, for the xQKZM data, we focus just on the bulk behavior to avoid boundary contributions, which are notable for the considered system sizes but are irrelevant in the thermodynamic limit. For that purpose, we only show the MPS data for the energy in the center of the 2D lattice. Furthermore, assuming that the bulk gap is not affected by boundary conditions, we use our numerically obtained gap from Fig. 2 and the exact value ΔE(g = 0) = 16J.

Beyond xQKZM

A central assumption in the xQKZM is that redistribution of occupations among eigenstates can be neglected. This naturally neglects all scattering processes leading to thermalization or phase-ordering kinetics, which might be especially relevant in the 2D context studied here. For that purpose, we study a slight variant of the dynamical protocol, which allows us to obtain some understanding of the influence of thermalization and phase-ordering kinetics onto the excitation energy density Q. Specifically, we interrupt the ramp when the transverse-field strength g reaches a value g/g = 1/2 for a waiting time t, where all Hamiltonian parameters are held constant, before continuing down toward g = 0. In this way, we provide further time for the system to relax and to redistribute occupations. The effect of the additional evolution on the outcome at g = 0 is shown in Fig. 7, where we include both the final excitation energy density ϵ and the final magnetization fluctuations 〈M2〉 = ∑(t, R) in the inset. We can see that the final excitation energy density tends to decrease with respect to t = 0, implying a kind of cooling effect due to the intermediate waiting interval. Furthermore, the magnetization fluctuations increase with t consistent with coarsening, i.e., the tendency of the system to develop ferromagnetic order at sufficiently low transverse fields.

Fig. 7.

Influence of free evolution in the ferromagnetic phase on the final state at g = 0.

Here, we consider a protocol as in Fig. 5, but with the ramp that stops at ϵ = 0.5 deep in the ferromagnetic phase for a waiting time t. The evolution then continues to g = 0 where the observables are measured. In the main panel, we show the total excitation energy per spin, and in the inset, we show correlations in the system measured by the variance of ferromagnetic magnetization. The energy decreases and quickly saturates with increasing t. In contrast, no such saturation is seen in the variance, where longer t provides more time for magnetic ordering.

We attribute this observed path dependence to thermalization dynamics and phase-ordering kinetics, which is not present in the 1D version of the model that is effectively noninteracting. It, however, becomes immediately relevant for the 2D case. One can understand the observed energy decrease by first deriving the equation of motion for the internal energy E(t) = 〈H(t)〉, which we will show here for the linear quench

While the energy itself is the sum of the spin-spin interaction and the transverse-field term, the energy change during the dynamics is governed by their difference. Now, it is crucial to realize that the thermalization dynamics and the resulting coarsening in our 2D Ising model are exactly characterized by a redistribution of energy between these individual contributions. It is a central consequence of the ramp starting on the paramagnetic side that the transverse magnetization is enhanced compared to the instantaneous equilibrium state. The accompanying thermalization dynamics is characterized by a transfer of energy from the transverse field to the interaction term. According to Eq. 17, this, in turn, also implies a negative change in the energy leading to a cooling effect. In other words, thermalization and coarsening dynamics, which naturally favor the interaction energy compared to the transverse-field contribution, directly affect the final energy at g = 0. These considerations also provide a direct interpretation of the results that we observe in Fig. 7. While the intermediate interruption of the ramp protocol itself does not change the energy directly, it gives the system the additional time t to thermalize toward the instantaneous equilibrium state, leading to an enhanced interaction and reduced field energy. However, this implies a larger right-hand side in the magnitude of Eq. 17 when the protocol is resumed after t, implying a stronger reduction of the total energy as compared to the case without stopping. Overall, Eq. 17 highlights, in combination with thermalization and coarsening properties of genuinely interacting 2D models, that there can be a noticeable energy change occurring during the parameter ramp in the ordered phase.

DISCUSSION

In this work, we have studied the QKZM in a 2D transverse-field Ising model using the combined effort of state-of-the-art numerical methods. As a main result, we have found universal defect production in the vicinity of the quantum critical point. For parameter sweeps deep into the ferromagnetic phase, we have introduced an xQKZM, which accounts for additional spectral information for the prediction of the final excitation energy densities.

The exponent that yields the best scaling collapse consistently deviates by about 7% from the precise values of the critical exponents that have been determined in previous numerical studies. We attribute this deviation to the fact that all simulations were limited by maximal feasible ramping rates τ or system sizes L. These limitations constrain our numerical experiments to probe a regime where corrections to the asymptotic universal scaling laws are relevant. In this regime, the gap opening is still approximately described by a power law, but the best-fitting exponent differs from the known value in the asymptotic limit. Hence, we expect better agreement of the dynamically observed critical exponents with results from studies in equilibrium for larger system sizes and slower ramping rates. Our results, among others, help to estimate the necessary parameters. For instance, avoiding the finite-size effects for τ ≃ 20, assuming in Fig. 4, would require the system of linear size L = 30.

Motivated in part by the remarkable progress in Rydberg atom and superconducting qubit quantum simulators, the theoretical analysis provided in this work exhibits a natural implementation in an experimental context. While evidence for QKZM has already been observed in systems of Rydberg atoms for a 1D quantum spin chain (), very recent developments enable the realization of the dynamics of transverse-field Ising models in 2D geometries involving hundreds of spin degrees of freedom (, ). Notice that although these systems typically realize antiferromagnetic spin interactions, the resulting dynamics is equivalent to that of an Ising ferromagnet upon transforming on every other lattice site. It is straightforward to implement the fully polarized initial condition (, , , ) and to temporally tune the couplings to realize the parameter sweep across the underlying quantum critical point (, , ). Furthermore, measurements in these systems yield spin configurations along a tunable axis, e.g., along σ or σ, which can give access to both the spin-spin correlation functions we considered (see Eq. 9) and the total Ising energy (see Eq. 5). Although quantum-optical systems such as Rydberg atoms exhibit remarkable isolation from the environment, decoherence nevertheless limits the experimentally accessible time scales. In this context, it is important to emphasize that the recent experiments have demonstrated already long coherence times, such as to observe the QKZM in one dimension () and close-to-adiabatic preparation of symmetry-broken low-energy states (, ). Overall, this makes our theoretical results directly experimentally accessible in Rydberg atom systems with the potential to address central questions that have remained open after the pioneering experiment on 2D QKZM in (). This concerns, for instance, the QKZ scaling deep in the ferromagnetic regime for larger system sizes and longer ramp times as well as immediate evidence for scaling behavior illustrated by a data collapse of the full correlation function.

METHODS

In this work, we use a combination of three state-of-the-art numerical algorithms. They provide mutual cross-checks where their ranges of applicability overlap and their full range is broader than any individual one.

The first, in order of appearance, is the 2D iPEPS tensor network on an infinite lattice. To simulate time evolution, we used the neighborhood tensor update (), which is a more efficient version of the code in (). Here, we use the second-order Suzuki-Trotter decomposition with time step dt = 0.01. Accuracy is limited in a controlled way by a bond dimension D of the iPEPS ansatz. All presented results appear converged for D = 8. Evaluation of expectation values requires approximating the infinite tensor environment whose accuracy is limited by an environmental bond dimension χ. The results for D = 8 were converged with χ ≤ 32. The iPEPS simulations fail after , where the convergence in D becomes insufficient.

The second method is based on representing the system’s wave function as a 1D MPS, where the 1D chain spans a 2D lattice. We obtain the best convergence (with respect to the required MPS bond dimensions D) using a diagonal steps–like covering. The time evolution is simulated using the TDVP algorithm of () (combining the one-site scheme with local application of two-site updates for enlarging bond dimensions) and the fourth-order time-dependent Suzuki-Trotter decomposition. We typically find the time step dt = 1/8 to be sufficient (with a few smaller time steps at the beginning to avoid instability of the TDVP when applied to the initial product state). We simulate the system sizes up to L = 14, with the bond dimension up to D = 512 to converge the presented results in various limits.

The third numerical approach uses NQS as variational ansatz for the wave function (). This approach was recently shown to be suited for the accurate simulation of quench dynamics in the 2D quantum Ising model (). The precision of NQS simulations is determined by the size of the neural network, which allows for systematic convergence checks. The time evolution was simulated using CNNs and regularization techniques introduced in (). For system sizes L ≤ 10, we used a single-layer network with fully connected filters. For larger systems, we used two-layer networks with square filters that have a diameter of L/2. We checked different network sizes for convergence and found that six channels, or four followed by three channels, were sufficient in the single- and two-layer case, respectively.

With the NQS, we can simulate finite systems with periodic boundary conditions in 2D, explicitly enforcing the invariance of the wave function under all lattice symmetries. This allows us to estimate boundary effects by comparing to MPS simulations, and the large system sizes reached convincingly reveal the power-law scaling of energy at the critical point. However, we found that the accuracy of the NQS simulations breaks down when continuing the ramp too far into the ferromagnetic phase; in particular, g = 0 is currently out of reach.

To compute the energy gap, we implemented the algorithm for excited state search introduced in (). The first step is to perform the ground-state search using the stochastic reconfiguration algorithm as already established in (). In the second step, the stochastic reconfiguration is modified such that the trial wave function is explicitly projected onto the subspace orthogonal to the ground state. We used the rescaled energy variance σ = 〈H2 − 〈H〉2〉/L4 to assert the accuracy of our result and reached for both the ground state and the first excited state in all cases σ < 10−4 and often σ ≈ 10−5 within 250 optimization steps of size δτ = 0.01.