Quasi-adiabatic quantum Monte Carlo algorithm for quantum evolution in imaginary time

Quasi-adiabatic quantum Monte Carlo algorithm for quantum evolution in imaginary time

Cheng-Wei Liu    Anatoli Polkovnikov    Anders W. Sandvik Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA
July 3, 2019
Abstract

We propose a quantum Monte Carlo (QMC) algorithm for non-equilibrium dynamics in a system with a parameter varying as a function of imaginary time. The method is based on successive applications of an evolving Hamiltonian to an initial state and delivers results for a whole range of the tuning parameter in a single run, allowing for access to both static and dynamic properties of the system. This approach reduces to the standard Schrödinger dynamics in imaginary time for quasi-adiabatic evolutions, i.e., including the leading non-adiabatic correction to the adiabatic limit. We here demonstrate this quasi-adiabatic QMC (QAQMC) method for linear ramps of the transverse-field Ising model across its quantum-critical point in one and two dimensions. The critical behavior can be described by generalized dynamic scaling. For the two-dimensional square-lattice system we use the method to obtain a high-precision estimate of the quantum-critical point , where is the transverse magnetic field and the nearest-neighbor Ising coupling. The QAQMC method can also be used to extract the Berry curvature and the metric tensor.

pacs:
05.30.âd, 03.67.Ac, 05.10.âa, 05.70.Ln

I Introduction

Quantum Monte Carlo (QMC) methods assaad07 (); sandvik10 () have become indispensable tools for ground-state and finite-temperature studies of many classes of interacting quantum systems, in particular those for which the infamous “sign problem” can be circumvented.kaul12 () In ground-state projector methods, an operator is applied to a “trial state” , such that approaches the ground state of the Hamiltonian when and an expectation value , with the norm , approaches its true ground-state value, . For the projector, one can use or a high power of the Hamiltonianbetanote (), . Here we will discuss a modification of the latter projector for studies of dynamical properties of systems out of equilibrium.

Real-time dynamics for interacting quantum systems is difficult to deal with computationally. Solving the Schrödinger equation directly, computations are restricted to very small system sizes by the limits of exact diagonalization. Despite progress with the Density-Matrix Renormalization Group (DMRG) PhysRevLett.69.2863 (); RevModPhys.77.259 () and related methods based on matrix-product states, this approach is in practice limited to one-dimensional systems and relatively short times. Efficiently studying long-time dynamics of generic interacting quantum systems in higher dimensions is still an elusive goal. However, recently, in Ref. [degrandi11, ], it was demonstrated that real-time and imaginary-time dynamics bear considerable similarities, and in the latter case, powerful and high-precision QMC calculations can be carried out on large system sizes for the class of systems where sign problems can be avoided.

Our work reported here is a further development of the method introduced in Ref. [degrandi11, ], where it was realized that a modification of the ground-state projector Monte Carlo approach with can be used to study non-equilibrium set-ups in quantum quenches (or ramps), where a parameter of the Hamiltonian depends on time according to an arbitrary protocol. By performing a standard Wick rotation of the time axis, a wave function is governed by the Shrödinger equation in imaginary time ( being real),

 ∂τ|ψ(τ)⟩=−H[λ(τ)]|ψ(τ)⟩. (1)

Here the Hamiltonian depends on the parameter through time, e.g.,

 H=H0+λ(τ)V, (2)

where and typically do not commute. The method is not limited to this form, however, and any evolution of can be considered. The Schrödinger equation has the formal solution

 |ψ(τ)⟩=U(τ)|ψ(τ0)⟩, (3)

where the imaginary-time evolution operator is given by

 U(τ)=Tτexp[−∫ττ0dτ′H[λ(τ′)]], (4)

where indicates time ordering. A time-evolved state and associated expectation values can be sampled using a generalized projector QMC algorithm. In Ref. [degrandi11, ] it was demonstrated that this non-equilibrium QMC (NEQMC) approach can be applied to study dynamic scaling at quantum phase transitions, and there are many other potential applications as well, e.g., when going beyond studies of finite-size gaps in “glassy” quantum dynamics and the quantum-adiabatic paradigm for quantum computing.

Here we introduce a different approach to QMC studies of quantum quenches, which gives results for a whole range of parameters in a single run (instead of just the final time), at a computational effort comparable to the previous approach. Instead of using the conventional time-evolution operator Eq. (4), we consider a generalization of the equilibrium QMC scheme based on projection with , acting on the initial ground state of with a product of evolving Hamiltonians:

 PM,1=[−H(λM)]....[−H(λ2)][−H(λ1)], (5)

where

 λt=λ0+tΔλ, (6)

and is the single-step change in the tuning parameter.gridnote () Here we will consider a case where the ground state of is known and easy to generate (stochastically or otherwise) and the ground states for other -values of interest are non-trivial. The stochastic sampling used to compute the evolution then takes place in a space representing path-integral-like terms contributing to the matrix element (the norm) . We will also later consider a modification of the method in which the ground state at the final point is known as well, in which case contributions to are sampled.

Staying with the doubly-evolved situation for now, we evaluate generalized expectation values after out of the operators in the product (5) have acted:

 ⟨A⟩t=⟨Ψ(λ0)|P1,MPM,t+1APt,1|Ψ(λ0)⟩⟨Ψ(λ0)|P1,MPM,1|Ψ(λ0)⟩. (7)

We will refer to this matrix element as an asymmetric expectation value, with the special case corresponding to a true quantum-mechanical expectation value taken with respect to an evolved wave function,

 |ψM⟩=PM,1|Ψ(λ0)⟩√⟨Ψ(λ0)|P1,MPM,1|Ψ(λ0)⟩, (8)

which approaches the ground state of the Hamiltonian for .

Away from the adiabatic limit, the evolved wave function Eq. (8) is, generally speaking, not the ground state of the equilibrium system. Nevertheless, as we demonstrate in detail in Sec. II, a quench velocity can be defined such that the symmetric expectation value in Eq. (7) approaches the expectation value after a conventional linear imaginary-time quantum quench with Eq. (4) done with the same velocity , if is low enough. In fact, the two quantities are the same to leading (linear) order in , not only in the strict adiabatic limit . We therefore name this scheme the quasi-adiabatic QMC (QAQMC) algorithm. Importantly, the leading corrections to the adiabatic evolution of the asymmetric expectation values for any contain important information about non-equal time correlation functions, very similar to the imaginary-time evolution.

The principal advantage of QAQMC over the NEQMC approach is that expectation values of diagonal operators in the basis used can be obtained simultaneously for the whole evolution path , by measuring in Eq. (7) at arbitrary points averagenote () (and one can also extend this to general off-diagonal operators, along the lines of Ref. [dorneich01, ], but we here limit studies to diagonal operators). The QAQMC scheme is also easier to implement in practice than the NEQMC method because there are no time integrals to sample.

As mentioned above, we will here have in mind a situation where the initial state is in some practical sense “simple,” but this is not necessary for the method to work—any state that can be simulated with standard equilibrium QMC methods can be used as the initial state for the dynamical evolution. The final evolved state can be very complex, e.g., for a system in the vicinity of a quantum-critical point or in a “quantum glass” (loosely speaking, a system with slow intrinsic dynamics due to spatial disorder and frustration effects). Here, as a demonstration of the correctness and utility of the QAQMC approach, we study generalized dynamic scaling in the neighborhood of the quantum phase transitions in the standard one-dimensional (1D) and 2D transverse-field Ising models (TFIMs).

As noted first in Ref. [degrandi11, ], the NEQMC method can be used to extract the components of the quantum metric tensor,provost_80 () the diagonal elements of which are the more familiar fidelity susceptibilities. Thanks to its ability to capture the leading non-adiabatic corrections to physical observables, the QAQMC approach can also be used for this purpose, and, as we will discuss briefly here and in more detail in Ref. [adi_long, ], one can also extract the Berry curvature through the imaginary antisymmetric components of the geometric tensor

The rest of the paper is organized in the following way. In Sec. II, we use adiabatic perturbation theory (APT) to demonstrate the ability of the QAQMC scheme to correctly capture the standard Schrödinger evolution in imaginary time, not only in the adiabatic limit but also including the leading corrections in the quench velocity. We show how these leading corrections correspond to the geometric tensor. In Sec. III, we discuss tests of the QAQMC scheme on 1D and 2D TFIMs, and also present a high-precision result for the critical field in the 2D model. In Sec. IV, we summarize our main conclusions and discuss future potential applications of the algorithm.

The key question we address in this section is whether the matrix element in Eq. (7) can give useful dynamical information for arbitrary “time” points in the sequence of operators. The expression only reduces to a conventional expectation value at the symmetric point , and even there it is not clear from the outset how computed for different relates to the velocity dependence of the expectation value based on the Schrödinger time-evolution operator in Eq. (4). Going away from the symmetric point brings in further issues to be addressed. For instance, there is no variational property of the asymmetric expectation value of the Hamiltonian for . Nevertheless, the approach to the adiabatic limit is well behaved and we can associate the leading deviations from adiabaticity with well defined dynamical correlation functions that appear as physical response in real time protocols. We show here, for the linear evolution Eq. (6), that one can identify a velocity such that a linear imaginary-time quench with in Eq. (6) gives the same results in the two approaches when , including the leading (linear) corrections in . For , the relevant susceptibilities in QAQMC defining non-adiabatic response are different than at but still well defined, contain useful information, and obey generic scaling properties.

In order to facilitate the discussion of the QAQMC method, we here first review the previous APT approach for Schrödinger imaginary-time dynamics degrandi11 (); adi_long () and then derive analogous expressions for the product-evolution. After this, we discuss some properties of the symmetric and asymmetric expectation values.

ii.1 Imaginary-time Schrödinger dynamics

The NEQMC method degrandi11 () uses a path-integral-like Monte Carlo sampling to solve the imaginary-time Shcrödinger equation Eq. (1) for a Hamiltonian with a time-dependent coupling. The formal solution at time is given by the evolution operator Eq. (4). In the strict adiabatic limit, the system will follow the instantaneous ground state, while in the slow limit one can anticipate deviations from adiabaticity, which will become more severe in gapless systems and, in particular, near phase transitions. Let us discuss the leading non-adiabatic correction to this imaginary-time evolution. The natural way to address this question is to use APT, similar to that developed in Refs. [ortiz_2008, ] and [PhysRevB.81.224301, ] in real time. We here follow closely the discussion of the generalization to imaginary time in Ref. [degrandi11, ].

We first write the wave function in the instantaneous eigenbasis of the time-dependent Hamiltonian :

 |ψ(τ)⟩=∑nan(τ)|n(λ(τ))⟩. (9)

We then substitute this expansion into Eq. (1),

 dandτ+∑mam(τ)⟨n|∂τ|m⟩=−En(λ)an(τ), (10)

where are the eigenenergies of the Hamiltonian corresponding to the states for this value of . Making the transformation

 an(τ)=αn(τ)exp[∫0τEn(τ′)dτ′], (11)

we can rewrite Eq. (1) as an integral equation;

 αn(τ)=αn(0)+∑m∫0τdτ′⟨n|∂τ′|m⟩αm(τ′) ×exp[−∫0τ′dτ′′(En(τ′′)−Em(τ′′))], (12)

where it should be noted that . In principle we should supply this equation with initial conditions at , but this is not necessary if is sufficiently large, since the sensitivity to the initial condition will then be exponentially suppressed. Instead, we can impose the asymptotic condition , which implies that in the distant past the system was in its ground state.

Eq. (12) is ideally suited for an analysis with the APT. In particular, if the rate of change is very small, , then to leading order in the system remains in its ground state; (except during the initial transient, which is not important because we are interested in large ). In the next higher order, the transition amplitudes to the states are given by;

 αn(0)≈−0∫−∞dτ⟨n|∂τ|0⟩exp[−∫0τdτ′Δn0(τ′)], (13)

where . The matrix element above for non-degenerate states can also be written as

 ⟨n|∂τ|0⟩=−⟨n|∂τH(τ)|0⟩/Δn0(τ). (14)

In what follows we will assume that we are dealing with a non-degenerate ground state.

To make further progress in analyzing the transition amplitudes Eq. (13), we consider the very slow asymptotic limit . To be specific, we assume that near the tuning parameter has the form (see also Ref. [PhysRevB.81.224301, ])

 λ(τ)≈λ(0)+vλ|τr|r!Θ(−τ). (15)

The parameter , which controls the adiabaticity, plays the role of the quench amplitude if , the velocity for , the acceleration for , etc. It is easy to check that in the asymptotic limit , Eq. (13) gives

 αn≈vλ⟨n|∂λ|0⟩(En−E0)r=−vλ⟨n|∂λH|0⟩(En−E0)r+1, (16)

where all matrix elements and energies are evaluated at . From this perturbative result we can in principle evaluate the leading non-adiabatic response of various observables and define the corresponding susceptibilities. For the purposes of comparing with the QAQMC approach, Eq. (16) suffices.

ii.2 Operator-product evolution

The quasi-adiabatic QMC method may appear very different from NEQMC but has a similar underlying idea. Instead of imaginary time propagation with Eq. (4), we apply a simple operator product to evolve the initial state. We first examine the state propagated with the first operators in the sequence in Eq. (5),

 |ψt⟩=[−H(λt)]…[−H(λ2)][−H(λ1)]|ψ0⟩, (17)

and after that we will consider symmetric expectation values of the standard form as well as the asymmetric expectation values in Eq. (7). We assume that the spectrum of is strictly positive, which is accomplished with a suitable constant offset to if needed.

ii.2.1 Linear protocols

The coupling can depend on the index in an arbitrary way. It is convenient to define

 τi=iT, (18)

where is the overall time scale, which can be set to unity. The leading non-adiabatic corrections will be determined by the system properties and by the behavior of near the point of measurement . The most generic is the linear dependence , where is related to the quench velocity (see below). In the end of this section we will briefly consider also more general nonlinear quench protocols.

Our strategy to analyze Eq. (17) in the adiabatic limit will be the same as in the preceding subsection. We first go to the instantaneous basis and rewrite

 |ψ(τi)⟩≡|ψi⟩=∑nan(τi)|n(λi)⟩≡∑nain|ni⟩. (19)

In the instantaneous basis, the discrete Schrödinger-like equation reads

 ai+1n=−∑maimEi+1n⟨ni+1|mi⟩, (20)

and it is instructive to compare this with Eq. (10). It is convenient to first make a transformation

 ain=t∏j=i+11(−Ejn)αin. (21)

This transformation does not affect the transition amplitude at the time of measurement : . Then the equation above becomes

 αi+1n=∑mαim[t∏j=i+1EjnEjm]⟨ni+1|mi⟩. (22)

Let us introduce a discrete derivative

 ⟨ni|←Δ≡⟨ni+1|−⟨ni|, (23)

and write the Schrödinger-like equation as

 αi+1n=αin+∑mαim[t∏j=i+1EjnEjm]⟨ni|←Δ|mi⟩. (24)

In the adiabatic limit, the solution of this equation is , i.e., the instantaneous ground state. To leading order of deviations from adiabaticity we find

 αi+1n=Cn+i∑k=0⎡⎣t∏j=k+1EjnEj0⎤⎦⟨nk|←Δ|0k⟩, (25)

where can be determined from the initial condition. In the limit of sufficiently large the initial state is not important so we should have for , so that . Therefore we find that the amplitude of the transition to the excited state is approximately

 αtn≈t−1∑k=0⎡⎣t∏j=k+1EjnEj0⎤⎦⟨nk|←Δ|0k⟩. (26)

Changing the summation index to we have

 αtn≈t∑p=1[t∏j=t+1−pEjnEj0]⟨nt−p|←Δ|0t−p⟩. (27)

It is clear that for large only terms contribute to the sum. In the extreme adiabatic limit one can thus move the matrix element outside of the summation and use the spectrum of the final Hamiltonian. In this case we find

 αtn ≈ EnE0⟨n|←Δ|0⟩1−En/E0 (28) = −EnΔλ⟨n|←∂λ|0⟩En−E0=EnΔλ⟨n|∂λ|0⟩En−E0,

where . By comparing Eqs. (16) and (28) we see that near the adiabatic limit QAQMC and NEQMC are very similar if . This can in principle always be ensured by having a sufficiently large energy offset, but even with a small offset we expect the ratio to be essentially constant for the range of contributing significantly when the spectrum becomes gapless close to a quantum-critical point. If the condition indeed is properly satisfied, then from Eqs. (16) and (28), we identify the quench velocity as

 vλ=E0Δλ. (29)

This is the main result of this section. We will confirm its validity explicitly in numerical studies with the QAQMC method in Sec. III. Since , where is the system size, we can also see that for a given total change in over the operators in the product.

Let us point out that Eq. (28) can be also rewritten as

 αtn≈−E0Δλ⟨n|←∂λ|0⟩En−E0−Δλ⟨n|←∂λ|0⟩. (30)

The first contribution here exactly matches that of Eq. (16) while the second term is an additional contribution corresponding to a sudden quench.

ii.2.2 Nonlinear protocols

We can extend the above result, Eq. (30), to arbitrary quench protocols. In particular, consider

 λt−p=λt+vλ(−E0)rpr−1(r−1)!, (31)

where (not necessarily an integer). For , we recover the linear protocol analyzed above. Then we can still rely on Eq. (27) but need to take into account that

 ⟨nt−p|←Δ|0t−p⟩ ≈ Δλt−p⟨nt|←∂λ|0t⟩ (32) = vλ(−E0)rpr−1(r−1)!⟨nt|←∂λ|0t⟩.

Thus, we find that

 αtn≈vλ(−E0)r(r−1)!Li1−r(En/E0)⟨nt|←∂λ|0t⟩, (33)

where is the Polylog function. In particular,

 Li0(x)=11−x, (34) Li−1(x)=x(1−x)2, (35) Li−2(x)=x(x+1)(1−x)3. (36)

Under the conditions discussed above (large offset or small energy gap) we again have and then we recover the continuum result using the fact that

 Li1−r(1−ϵ)≈(r−1)!ϵr.

Then, indeed,

 αtn=atn ≈ vλ(−E0)r(r−1)!(r−1)!(1−En/E0)r⟨nt|←∂λ|0t⟩ (37) = vλ⟨nt|←∂λ|0t⟩(En−E0)r,

which exactly matches Eq. (16).

ii.3 Expectation values

While asymptotically Eq. (7) gives the ground state of the observable in the adiabatic limit for all values of , the approach to this limit as is qualitatively different depending on whether is equal to or not. More precisely, if where as , we encounter two different asymptotic regimes for and .

ii.3.1 Symmetric expectation values; t=M

In this limit the expectation value of the observable in the leading order of the adiabatic perturbation theory reduces to

 ⟨A⟩t=M≈⟨ψ(vλ)|A|ψ(vλ)⟩, (38)

where is the imaginary time velocity identified earlier. For generic observables not commuting with the Hamiltonian, we find

 ⟨A⟩t=M≈⟨A⟩0+vλχ′Aλ, (39)

where

 χ′Aλ=∑n≠0⟨0|A|n⟩⟨n|∂λ|0⟩En−E0+c.c. (40)

is the susceptibility. All energies and matrix elements are evaluated at “time” .

For diagonal observables , like the energy or energy fluctuations, we have

 ⟨A⟩t=M≈⟨A⟩0+v2λ∑n≠0|⟨n|∂λ|0⟩|2(En−E0)2⟨n|A|n⟩. (41)

In particular, the correction to the energy is always positive as it should be for any choice of wave function deviating from the ground state. Let us emphasize that for diagonal observables the leading non-adiabatic response at the symmetric point in imaginary time coincides with that in real time, and, thus QAQMC or NEQMC can be used to analyze real time deviations from adiabaticity, as was pointed out in the case of NEQMC in Ref. [degrandi11, ].

ii.3.2 Asymmetric expectation value, t≠M

It turns out that the asymptotic approach to the adiabatic limit is quite different for non-symmetric points with . Without loss of generality we can focus on (since all expectation values are symmetric with respect to for the symmetric protocol we consider averagenote ()). Then the expectation value of is evaluated with respect to different eigenstates

 ⟨A⟩t=⟨ψL|A|ψR⟩⟨ψL|ψR⟩, (42)

where

 |ψR⟩ = H(λt)⋯H(λ2)H(λ1)|ψ0⟩, |ψL⟩ = H(λt+1)⋯H(λM−1)H(λM)PM,1|ψ0⟩. (43)

Note that the overlap is independent of by construction and is real.

It is easy to see that for diagonal observables we obtain a leading asymptotic as in Eq. (41) but with the opposite sign in the second term

 ⟨A⟩t≠M≈⟨A⟩0−v2λ∑n≠0|⟨n|∂λ|0⟩|2(En−E0)2⟨n|A|n⟩. (44)

In particular, the leading correction to the ground state energy is negative when deviates sufficiently from the symmetric point, i.e., . There is no contradiction here since the left and right states are different (i.e., we are not evaluating a true expectation value and there is no variational principle). Both Eqs. (41) and (44) recover the exact result in the adiabatic limit. Since the correction up to the sign exactly matches the real time result, we can still use the non-symmetric expectation value for diagonal observables to extract the real time non-adiabatic response. For , the sign of the correction should change, to connect smoothly to the variational expectation value. The crossover between positive and negative corrections to the energy takes place around a point that asymptotically converges to in the adiabatic limit (where the deviation from the ground-state energy at vanishes). We will illustrate this with numerical results in Sec. III.1 (see Fig. 1).

As in the symmetric case, using the APT discussed in the previous section the results derived here easily extend to other values of the exponent .

ii.3.3 The metric tensor and Berry curvature

If , then the susceptibility Eq. (40) reduces to the component of the metric tensor,degrandi11 (); adi_long () which, thus, can be readily extracted using the QAQMC algorithm. In particular, the diagonal components of the metric tensor define the more familiar fidelity susceptibility.

Next, we observe that for sufficiently different from , the approach to the ground state in the left function in Eq. (43) effectively corresponds to a change in sign of the velocity, and, thus, we find

 ⟨A⟩t≈⟨ψ(−v)|A|ψ(v)⟩⟨ψ(−v)|ψ(v)⟩, (45)

where the wave functions and are evaluated at the same value of the coupling determined by the value of . We can use the results of the previous section to find that for off-diagonal observables

 ⟨A⟩t≈⟨A⟩0−ivλχ′′Aλ, (46)
 χ′′Aλ=i∑n≠0⟨0|A|n⟩⟨n|∂λ|0⟩En−E0−c.c. (47)

Based on this result we conclude that the leading non-adiabatic correction is imaginary and coincides, up to the factor of imaginary , with the real-time non-adiabatic correction.adi_long () In particular, for the susceptibility is proportional to the Berry curvature.

The fact that we are getting the opposite sign (compared to the real time protocol) in the susceptibility for diagonal observables and the Berry curvature for off-diagonal observables away from the symmetric points in Eqs (44) and (46) is a consequence of general analytic properties of the asymmetric expectation values. As we discuss in Ref. [adi_long, ] the expectation value Eq.(45) is the analytic continuation of the real time expectation value to the imaginary velocity . This continuation is valid in all orders of expansion of the expectation value of in .

Iii Results

As a demonstration of the utility of QAQMC and the behaviors derived in the previous section we here study the TFIM, defined by the Hamiltonian

 H=−s∑⟨i,j⟩σziσzj−(1−s)∑iσxi, (48)

where are nearest-neighbor sites, and and are Pauli matrices. Here, plays the role of the tuning parameter, which in the simulations reported below will vary between (where the ground state is trivial) to a value exceeding the quantum-critical point; in a 1D chain and in the 2D square lattice.J.Phys.A.33.6683 ()

We work in the standard basis of eigenstates of all . The simulation algorithm samples strings of diagonal and off-diagonal terms in Eq. (48), in a way very similar to the stochastic series expansion (SSE) method, which has been discussed in detail in the case of the TFIM in Ref. [PhysRevE.68.056701, ]. The modifications for the QAQMC primarily concern the sampling of the initial state, here , which essentially amounts to a particular boundary condition replacing the periodic boundaries in finite-temperature simulations. An SSE-like scheme with such modified boundaries was also implemented for the NEQMC method in Ref. [degrandi11, ], and recently also in a study of combinatorial optimization problems in Ref. [farhi12, ]. We here follow the same scheme, using cluster updates in which clusters can be terminated at the boundaries. The implementation for the product with varying coupling is even simpler than SSE or NEQMC, with the fixed-length product replacing the series expansion of Eq. (4). The changes relative to Refs. [degrandi11, ] and [PhysRevE.68.056701, ] are straightforward and we therefore do not discuss the sampling scheme further here.

iii.1 Cross-over of the energy correction

As we discussed in Sec. (II), the asymmetric expectation value (7) of the Hamiltonian has a negative correction to the ground-state energy when is sufficiently away from the symmetric point . In Fig. 1 we illustrate this property and the convergence to the ground-state energy for all with increasing with simulation data for a small 1D TFIM system. We here plot the results versus the rescaled propagation power . The region of negative deviations move toward the symmetric point with increasing . Note that the deviations here are not strongly influenced by the critical point (which is within the parameter simulated but away from the symmetric point), although the rate of convergence should also be slow due to criticality. The rate of convergence to the ground state can be expected to be (and is here seen to be) most rapid for and .

iii.2 Quantum-critical dynamic scaling

The idea of dynamic scaling at a critical point dates back to Kibble and Zurek for quenches (also called ramps, since the parameter does not have to change suddenly, but linearly with arbitrary velocity as a function of time) of systems through classical phase transitions. J.Phys.A.9.1387 (); Nature.317.505 () Here, the focus was on the density of defects. The ideas were later generalized also to quantities more easily accessible in experiments, such as order parameters, and the scaling arguments were also extended to quantum systems.PhysRevB.72.161201 (); RevModPhys.83.863 (); PhysRevLett.95.105701 (); PhysRevLett.95.245701 () The basic notion is that the system has a relaxation time , and if some parameter (here a parameter of the Hamiltonian) is changed such that a critical point is approached, the system can stay adiabatic (or in equilibrium) only if the remaining time to reach the critical point is much larger than the relaxation time, . In general, one expects , where is the correlation length, the exponent governing its divergence with the distance to the critical point, and the dynamic exponent. For a system of finite size (length) , is maximally of order and, thus, for a linear quench the critical velocity separating slow and fast power-law quenches according to Eq. (15) should heuristically be given by , and for a power-law quench with exponent according to Eq. (15) this generalizes to PhysRevB.81.224301 ()

 vcrit∼L−(zr+1/ν). (49)

One then also expects a generalized finite-size scaling form for singular quantities ,

 A(L,ϵ)=Lκf(ϵL1/ν,vLzr+1/ν), (50)

where characterizes the leading size-dependence at the critical point of the quantity considered. For , Eq. (50) reduces to the standard equilibrium finite-size scaling hypothesis. This scaling was recently suggested and tested in different systems, both quantum deng_08 (); PhysRevB.81.224301 (); kolodrubetz_12 () and classical chandran_12 ().

The above expression Eq. (50) combined with the product-evolution Eq. (5) allows us to study a phase transition based on different combinations of scaling in the system size and the velocity in non-equilibrium setups. For example, if one wants to find the critical point for the phase transition and the exponent is known, one can carry out the evolution under the critical-velocity condition:

 vLz+1/ν=c, (51)

where is a constant. In this paper, we focus on linear quench protocols and set henceforth. As we discussed in Sec. II.2, the QAQMC method applied to a system of size (volume) based on evolution with operators in the sequence and change between each successive operator corresponds to a velocity , with the prefactor depending on the ground state energy (at the critical point). The exact prefactor will not be important for the calculations reported below, and for convenience in this section, we define

 v=sfNM, (52)

where is the final value of the parameter in Eq. (48) over the evolution (which is also the total change in , since we start with the eigenstate at ). The critical product-length is, thus, given by

 M=1cNLz+1/ν=1cLd+z+1/ν, (53)

where we have also for simplicity absorbed into .

Using an arbitrary of order in Eq. (51), the critical point can be obtained based on a scaling function with the single argument in Eq. (50). We will test this approach here, in Secs. III.4 and III.5, and later, in Sec. III.6, we will show that exact knowledge of the exponents in Eq. (51) is actually not needed. First, we discuss the quantities we consider in these studies.

iii.3 Quantities studied

We will focus our studies here on the squared -component magnetization (order parameter),

 (54)

We can also define a susceptibility-like quantity (which we will henceforth refer to as the susceptibility) measuring the magnetization fluctuations:

 χ=N(⟨m2z⟩−⟨|mz|⟩2). (55)

Here both terms have the same critical size-scaling as the equal-time correlation function;

 ⟨mz⟩2∼⟨|mz|⟩2∼L−(d+z−2+η), (56)

where is the spatial dimensionality. The prefactors for the two quantities are different, however, a divergent peak remains in Eq. (55) at the transition. Away from the critical point with increasing system size.

To clarify our use of , we point out that we could also just study the scaling of , but the peak produced when subtracting off the second term in Eq. (55) is helpful in the scaling analysis. According to Eq. (50) and using in Eq. (56), the full scaling behavior of the fluctuation around the critical point should follow the form

 χ∼L1−ηf((s−sc)L1/ν,vL1+1/ν), (57)

for any dimensionality .

We should point out here that the true thermodynamic susceptibility based on the Kubo formula PhysRevB.43.5950 () (imaginary-time integral) yields a stronger divergence . This quantity is, however, more difficult to study with the QAQMC algorithm, because, unlike in standard finite- QMC methods, the time integration cannot simply be carried out within the space of time-evolving Hamiltonians in Eq. (5) and Eq. (7). The standard Feynman-Suzuki correspondence between the -dimensional quantum and -dimensional classical systems is not realized in our scheme. The configuration space of time-evolving Hamiltonians builds in the relaxation time, , in a different way, not just in terms of equilibrium fluctuations in the time direction, but in terms of evolution as a function of a time-dependent parameter.

A useful quantity to consider for extracting the critical point is the Binder cumulant, PhysRevLett.47.693 (),

 U=32(1−13⟨m4z⟩⟨m2z⟩2). (58)

For a continuous phase transition, converges to a step function as . The standard way to analyze this quantity for finite is to graph it versus the argument for different and extract crossing points, which approach the critical point with increasing . Normally, this is done in the equilibrium, either by taking the limit of the temperature for each first, or by fixing if is known. Here, the latter condition is replaced by Eq. (51), but, as we will discuss further below, the condition can be relaxed and the exponents do not have to be known accurately a priori. Our approach can also be used to determine the exponents, either in a combined procedure of simultaneously determining the critical point and the exponents, or with a simpler analysis after first determining the critical point.

We have up until now only considered calculations of equal-time observables, but, in principle, it is also possible and interesting to study correlations in the evolution direction, which also can be used to define susceptibilities.

In the following we will illustrate various scaling procedures using results for the 1D and 2D TFIMs. The dynamic exponent is known for both cases, and in the 1D case all the exponents are rigorously known since they coincide with those of the classical 2D Ising model. For the 2D TFIM, the exponents are know rather accurately based on numerics for the 3D classical model.

iii.4 1D transverse-field Ising model.

The 1D TFIM provides a rigorous testing ground for the new algorithm and scaling procedures since it can be solved exactly.sachdev_qpt () The critical point corresponds to the ratio between the transverse field and the spin-spin coupling equaling , i.e., in the Hamiltonian Eq. (48). The critical exponents, known through the mapping to the 2D Ising model,Prog.Theor.Phys.56.1454 () are and .

The results presented here were obtained in simulations with the parameter spanning the range with , i.e., going from the trivial ground state of the field term to well above the critical point. Fig. 2 shows examples of results for the susceptibility and the Binder cumulant. The operator-sequence length , Eq. (5), was scaled with the system size in order to stay at the critical velocity according to Eq. (53). We emphasize again that a single run produces a full curve within the -range used. In order to focus on the behavior close to criticality, we have left out the results for small in Fig. 2. Since is very large (up to for the largest in the cases shown in the figure), we also do not compute expectation values for each in Eq. (7), but typically spacing measurements by operators.

Extracting Binder curve-crossings using system-size pairs and , with , and extrapolating the results to , we find , as illustrated in Fig.(3). Thus, the procedure produces results in full agreement with the known critical point.

The dynamical scaling of the susceptibility is illustrated in Fig. 4. Here, there are no adjustable parameters at all, since all exponents and the critical coupling are known (and we use the exact critical coupling , although the numerical result extracted below is very close to this value and produces an almost identical scaling collapse). While some deviations from a common scaling function are seen for the smaller systems and far away from the scaled critical point , the results for larger sizes and close to the peak rapidly approach a common scaling function. This behavior confirms in practice our discussion of the definition of the velocity and the ability of the QAQMC method to correctly take into account at least the first corrections to the adiabatic evolution.

iii.5 2D transverse-field Ising model

The 2D transverse-field Ising model provides a more serious test for our algorithm since it cannot be solved exactly. Among many previous numerical studies,J.Phys.C.4.2370 (); PhysRevB.57.8494 (); J.Phys.A.33.6683 () Ref. [J.Phys.A.33.6683, ] arguably has the highest precision so far for the value of the critical coupling ratio. Exact diagonalization was there carried out for up to lattice size. In terms of the critical field in units of the coupling , the critical point was determined to , where the error bar reflects estimated uncertainties in finite-size extrapolations. Results based on QMC simulations J.Phys.C.4.2370 (); PhysRevB.57.8494 () are in agreement with this value, but the statistical errors are larger than the above extrapolation uncertainty. One might worry that the system sizes are very small and the extrapolations may not reflect the true asymptotic size behavior. However, the data points do follow the functional forms expected based on the corresponding low-energy field theory, and there is therefore no a priory reason to question the results. It is still useful to try to reach similar or higher precision with other approaches, as we will do here with the QAQMC method combined with dynamic scaling.

In this case we simulate the linear quench in the window of , which contains the previous estimates for the critical value as discussed above. Although we could also carry out an independent scaling analysis to extract the critical exponents, we here choose to just use their values based on previous work on the classical 3D Ising model; , and .PhysRevB.59.11471 () Our goal here is to extract a high-precision estimate of the critical coupling, and, at the same time, to further test the ability of QAQMC to capture the correct critical scaling behavior. We again scale with according to Eq. (53), with the constant .

As in the 1D case, we extract Binder-cumulant crossing points based on linear system sizes and with . Fig. 5 shows the results versus along with a fit to a power-law correction PhysRevLett.47.693 () for . Extrapolating to infinite size gives , which corresponds to a critical field (in unit of ) . This is in reasonable good agreement with the value obtained in Ref. [J.Phys.A.33.6683, ] and quoted above, with our (statistical) error bar being somewhat smaller. To our knowledge, this is the most precise value for the critical coupling of this model obtained to date. We emphasize that we here relied on the non-equilibrium scaling ansatz to extract the equilibrium critical point. Allowing for deviations from adiabaticity in a controlled way and utilizing the advantage of the QAQMC algorithm allowed us to extract observables in the whole range of couplings in a single run. This requires considerably less computational resources than standard equilibrium simulations, which must be repeated for several different couplings in order to carry out the crossing-point analysis.

Fig. 6 shows the susceptibility scaled according to the behavior expected with Eq. (50) when the second argument is held constant. As in the 1D case, the data converge rapidly with increasing size toward a common scaling function in the neighborhood of the transition point, again confirming the correct quasi-adiabatic nature of the QAQMC method.

iii.6 Further tests

The results discussed in the preceding subsections were obtained with the KZ velocity condition Eq. (51), applied in the form of Eq. (53) tailored to the QAQMC approach, with specific values for the constant . In principle, the constant is arbitrary, but the non-universal details of the scaling behavior depend on it. This is in analogy with a dependence on the shape, e.g., an aspect ratio, of a system in equilibrium simulations at finite temperature, or to the way the inverse temperature is scaled as with arbitrary in studies of quantum phase transitions (as an alternative to taking the limit for each lattice size). The critical point and the critical exponents should not depend on the choices of such shape factors or limiting procedures.

To extract the critical coupling, in the preceding subsections, we fixed the exponents and at their (approximately) known values, and one may at first sight assume that it is necessary to use their correct values. It is certainly some times convenient to do so, in order to set the second argument of the scaling function Eq. (50) to a constant and, thus, obtain a simpler scaling function depending on a single argument. However, one can study critical properties based on the scaling approach discussed above as long as the velocity approaches zero as the system size increases. This observation can be important in cases where the critical exponents are not known and one would like to obtain an accurate estimate of the critical coupling before carrying out a scaling analysis to study exponents. We will test this in practice here. As we will discuss further below, one should use a different power in the scaling ansatz Eq. (50) if the velocity is brought to zero slower than the critical form.

In cases where we use the “wrong” values of the exponents, we formally replace by a free parameter ,

 v∼L−α/c, (59)

and the corresponding substitution in Eq. (53). To understand the scaling of the observables for arbitrary , we return to the general scaling form given by Eq. (50). In the case of the Binder cumulant and for linear quench protocol, this form reads

 U=f((s−sc)L1/ν,vLz+1/ν). (60)

As we pointed out above, when the velocity scales exactly as , the dependence on the second argument in the scaling function drops out and we can find the crossing point in a standard way as we did in Figs. 3 and 5. Suppose that we do not know the exponents and a priory and instead scale as in Eq. (59). Then there are three possible situations: (i) , (ii) , and (iii) , where we already have analyzed scenario (i). In scenario (ii), where velocity scales to zero faster than the critical KZ velocity, the second argument of the scaling function approaches zero as the system size increases and, thus, the scaling function effectively approaches the equilibrium velocity-independent form. We can then extract the crossing point as in the first scenario, and this gives the correct critical coupling in the limit of large system sizes. Finally, in case (iii) the velocity scales zero slower than the critical KZ value and the second argument in Eq. (60) diverges, which implies that the system enters a strongly non-equilibrium regime. This scenario effectively corresponds to taking the thermodynamic limit first and the adiabatic limit second. Then, if the system is initially on the disordered side of the transition, the Binder cumulant vanishes in the thermodynamic limit. At finite but large system sizes its approach to zero should be given by the standard Gaussian theory:

 U≈CLd. (61)

Combining this with the scaling ansatz Eq. (60) we find that for , the expected asymptotic of the Binder cumulant is

 U≈L−dv−d/(z+1/ν)~f((s−sc)L1/ν), (62)

where is some other velocity independent scaling function. Thus we can find the correct transition point by finding crossing points of . Similar considerations apply to the ordered side of the transition, where the Binder cumulant approaches one as the inverse volume.

The three cases are illustrated in the lower panel of Fig. 7, which shows Binder-cumulant crossings extracted from appropriately scaled data in cases (i), (ii), and (iii) above. Additionally, to illustrate the insensitivity to the choice of the constant in the scaled sequence length in Eq. (53), results based on two different constants are shown for case (i). In all cases, the extrapolated critical couplings agree with each other to within statistical errors. Note that, on the one hand, if the exponent gets very large, then the time of simulations, which scales as , rapidly increases with the system size and the algorithm becomes inefficient. On the other hand, if is very small, our results indicate that the size dependence is larger and it is more difficult to carry out the extrapolation to infinite size. The optimal value of should be as close as possible to the critical KZ power, but to be on the safe side when scaling according to the standard KZ critical form, case (i), one may choose a somewhat larger value, since the subcritical velocity in case (ii) has the same scaling form.

Next we illustrate how the same idea works in the case of the order parameter. Around the critical point (, ), the squared magnetization [see Eq. (54)] can be written as

 m2z = L−2β/νf((s−sc)L1/ν,vLz+1/ν). (63)

As in the previous discussion we scale and depending on the exponent there are two different asymptotics of the scaling function. For the second argument vanishes or approaches constant so we effectively get the equilibrium scaling

 m2z = L−2β/νf((s−sc)L1/ν) (64)

If, conversely, then on the disordered side of the transition scales as . This immediately determines the asymptotic of the scaling function in Eq. (63):

 m2z=L−dv(2β/ν)−dz+1/ν~f((s−sc)L1/ν). (65)

Equation (65) can be used in the same way as the Binder cumulant to extrapolate the critical point, using the standard form PhysRevLett.47.693 () for the rescaled . As shown in the top panel of Fig. (7), after rescaling the order parameter and extrapolating the crossing points between the appropriately rescaled curves to the thermodynamics limit, all curves, obtained from below or above the adiabatic limit Eq. (49), converge to the same value . This approach also suggests a way to determine the transition point in experiment, since one can sweep through the critical point at different velocities, the crossing point can then be extracted through the measurement of the order parameter. It is also worth mentioning that since one can extrapolate the critical point independently without knowing the actual exponent prior to the simulation, an optimization procedure can be carried out to determine the exponents posterior to the simulation.classical_quench.unpublished ()

For completeness we also briefly discuss the role of the final point of the evolution. Fig. 8 shows 2D results for the squared magnetization Eq. (54) and susceptibility Eq. (55) obtained for a range of final points above the critical value. Here the velocity was kept constant for all the cases. The values of the computed quantities at some fixed , e.g., at , show a weak dependence on for the lowest- runs. The deviations are caused by contributions of order and higher, which are non-universal as discussed in Sec. II.2. For very high velocities the dependence on can be much more dramatic than in Fig. 8, but this is not the regime in which the QAQMC should be applied to study universal physics.

Iv Summary and Discussion

We have presented a nonequilibrium QAQMC approach to study quantum dynamics, with a simple product of operators with evolving coupling replacing the standard Schrödinger time evolution. We showed that this approach captures the leading non-adiabatic corrections to the adiabatic limit, both by analytical calculations based on the APT and by explicit simulations of quantum-critical systems with the QAQMC algorithm. The simulation results obey the expected generalized dynamic scaling with known static and dynamic critical exponents. We also extended the scaling formalism beyond results obtained previously in Ref. [degrandi11, ]. We analyzed the leading non-adiabatic corrections within this method and showed that they can be used to extract various non-equal time correlation functions, in particular, the Berry curvature and the components of the metric tensor. A clear advantage of the QAQMC approach is that one can access the whole range of couplings in a single run. Being a simple modification of projector QMC, the QAQMC method is applicable to the same class of models as this conventional class of QMC schemes—essentially models for which “sign problems” can be avoided.

As an illustration of the utility of QAQMC, we applied the algorithm and the scaling procedures to the 1D and 2D TFIMs. The expected scaling behaviors are observed very clearly. In the 1D case we extracted a critical coupling in full agreement with the known value, and in 2D we obtained an estimate with unprecedented (to our knowledge) precision (small error bars); . Based on repeating the fitting procedures with different subsets of the data, we believe that any systematical errors due to subleading corrections neglected in the extrapolations should be much smaller than the statistical errors, and, thus, we consider the above result as unbiased.

The QAQMC approach bears some similarities to previous implementations of quantum annealing within QMC algorithms.Santoro06 (); Bapst12 () However, the previous works have mainly considered standard equilibrium QMC approaches in which some system parameter is changed as a function of the simulation time. This evolution is not directly related to true quantum dynamics (and, thus, is not really quantum annealing), but is dependent on the particular method used to update the configurations. In contrast, in our scheme, as in the NEQMC method introduced in Ref. [degrandi11, ], the evolution takes place within the individual configurations, and there is a direct relationship to true Schrödinger evolution in imaginary time.

In Green’s function (GF) QMC simulations the gradual change of a system parameter with the simulation time is rather closely related to the QAQMC scheme (since also there one applies a series of terms of the Hamiltonian to a state), with the difference being that QAQMC uses true importance sampling of configurations, with no need for guiding wave functions and no potential problems related to mixed estimators. Our asymmetric expectation values could be considered as a kind of mixed estimator as well, but we have completely characterized them within the APT. In addition, the previous uses of GFQMC with time-evolving Hamiltonians have, to our knowledge, never addressed the exact meaning of the velocity of the parameter evolution. The correct definition of the velocity is of paramount importance when applying quantum-critical scaling methods, as we have discussed here. We have here computed the velocity within APT for the QAQMC scheme. The same relationship with Schrödinger dynamics may possibly hold for GFQMC as well, but, we have not applied the APT to this case and it is therefore not yet clear whether GFQMC can capture correctly the same universal non-equilibrium susceptibilities as the QAQMC and NEQMC methods. We expect QAQMC to be superior to time-evolving GFQMC, because of its better control over measured symmetric and asymmetric expectation values and fully realized importance sampling.

Some variants of GFQMC use true importance sampling, e.g., the Reptation QMC (RQMC) method,PhysRevLett.82.4745 () which also avoids mixed estimators. The configuration space and sampling in the QAQMC method bears some similarities with RQMC, recent lattice versions of which also use SSE-inspired updating schemes.PhysRevE.82.046710 () However, to our knowledge, imaginary-time evolving Hamiltonians have not been considered in RQMC and in other related variants of GFQMC, nor has the role played by the velocity when crossing the quantum critical point been stressed. This has so far been our focus in applications of the QAQMC and NEQMC methods. In principle one could also implement the ideas of time-evolution similar to QAQMC within the RQMC approach.

We also stress that we have here not focused on optimization. Previous works on quantum annealing within QMC schemes have typically focused on their abilities to optimize difficult classical problems. While the QAQMC may potentially also offer some opportunities in this direction, our primary interest in the method is to use it to extract challenging dynamical information under various circumstances. A recent theoretical analysis of optimization within sign-problem free QMC approaches Hastings13 () is not directly applicable to the QAQMC and NEQMC approaches but generalizations should be possible.

The QAQMC and NEQMC methods provide correct realizations of quantum annealing in imaginary time. Besides their ability to study dynamic scaling, with exponents identical to those in real-time Schrödinger dynamics,degrandi11 () it will be interesting to explore what other aspects of real-time dynamics can be extracted with these methods. In particular, their applicability to quantum glasses, of interest in the context of quantum adiabatic computing farhi12 () as well as in condensed matter physics, deserves further studies.

The ability of the QAQMC to produce results for a whole evolution path in a single run can in principle also be carried over to the conventional Schrödinger imaginary-time evolution with in Eq. (4). By “slicing” the time evolution into successive evolutions over a time-segment ,

 U(τ)=K∏n=1Tτexp[−∫τnτn−1dτH[λ(τ)]], (66)

where , one can evaluate matrix elements analogous to Eq. (7) by inserting the operator of interest at any point within the product of time-slice operators in . In this case, the symmetric expectation value, evaluated at the mid-point, is identical to the NEQMC method,degrandi11 () and the asymmetric expectation values will exhibit properties similar to those discussed in Sec. II.3.2. We have not yet explored this approach, and it is not clear whether it would have any other advantage besides the exact reduction to Schrödinger dynamics of the symmetric expectation values. In practice the simulations will be more complex than the QAQMC approach because of the need to sample integrals, but not much more so than the NEQMC method. It should be relatively easy to adapt the RQMC method with an evolving Hamiltonian in this formulation of the time-evolution.

Finally, we point out that, in principle, one can also carry out a one-way evolutions with the QAQMC algorithm. Instead of starting with the eigenstate at both and and then projecting them to the eigenstate using two sequences of the form Eq. (5), one can make correspond to and let it evolve to corresponding to with only a single operator sequence of length . In the case of the TFIM Eq. (48), the obvious choice is then to evolve from to (the classical Ising model), so that both edge states are trivial. All our conclusions regarding the definition of the velocity and applicability of scaling form remain valid in this one-way QAQMC. Results demonstrating this in the case of the 1D TFIM are sown in Fig. 9. We anticipate that this approach may be better than the two-way evolution in some cases, but we have not yet compared the two approaches extensively.

Acknowledgements.
We acknowledge collaboration with Claudia De Grandi in a related work and would like to thank Mike Kolodrubetz for valuable comments. This work was supported by the NSF under Grant No. PHY-1211284.

References

• (1) F. F. Assaad and H. G. Evertz, Lect. Notes Phys. 739, 277 (2007).
• (2) A. W. Sandvik, AIP Conf. Proc. 1297, 135 (2010).
• (3) R. K. Kaul, R. G. Melko, and A. W. Sandvik, Annu. Rev. Cond. Mat. Phys. 4, 179 (2013).
• (4) Here gives the same rate of convergence for the two choices for a given system size . This can be shown, e.g., by a Taylor-expansion of the exponential, which for large is dominated by powers , where is the ground state energy.
• (5) S. R. White, Phys. Rev. Lett. 69, 2863 (1992)
• (6) U. Schollwöck, Rev. Mod. Phys. 77, 259 (2005).
• (7) C. De Grandi, A. Polkovnikov, and A. W. Sandvik, Phys. Rev. B 84, 224303 (2011).
• (8) In principle one can also consider a nonlinear grid of “time” points, but here we will consider the simplest case of a uniform grid.
• (9) Eq. (7) has been defined for , but we can also define the asymmetric expectation value for by placing the operator within the product . Clearly, we have the symmetry .
• (10) A. Dorneich and M. Troyer, Phys. Rev. E 64, 066701 (2001).
• (11) J. P. Provost and G. Vallee, Comm. Math. Phys., 76, 289 (1980).
• (12) C. De Grandi, A. Polkovnikov and A. Sandvik, arXiv:1301.2329.
• (13) G. Rigolin, G. Ortiz, and V. H. Ponce, Phys. Rev. A 78, 052508 (2008).
• (14) C. De Grandi, V. Gritsev, and A. Polkovnikov, Phys. Rev. B 81, 224301 (2010).
• (15) C. J. Hamer, J. Phys. A: Math. Gen. 33, 6683 (2000).
• (16) A. W. Sandvik, Phys. Rev. E 68, 056701 (2003).
• (17) E. Farhi, D. Gosset, I. Hen, A. W. Sandvik, P. Shor, A. P. Young, and F. Zamponi, Phys. Rev. A 86, 052334 (2012).
• (18) T. W. B. Kibble, J. Phys. A: Math. Gen. 9, 1387 (1976).
• (19) W. H. Zurek, Nature 317, 505 (1985).
• (20) A. Polkovnikov, Phys. Rev. B 72, 161201 (2005).
• (21) A. Polkovnikov, K. Sengupta, A. Silva and M. Vengalattore, Rev. Mod. Phys. 83, 863 (2011)
• (22) W. H. Zurek, U. Dorner, and P. Zoller, Phys. Rev. Lett. 95, 105701 (2005).
• (23) J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005).
• (24) S. Deng, G. Ortiz, L. Viola, Europhys. Lett. 84, 67008 (2008).
• (25) M. Kolodrubetz, B. K. Clark, and D. A. Huse, Phys. Rev. Lett., 109, 015701, (2012).
• (26) A. Chandran, A. Erez, S. S. Gubser, S. L. Sondhi, Phys. Rev. B 86, 064304 (2012)
• (27) A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991).
• (28) K. Binder, Phys. Rev. Lett. 47, 693 (1981).
• (29) S. Sachdev, Quantum Phase Transitions (Cambridge University Press, Cambridge, 2012) .
• (30) M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976).
• (31) P. Pfeuty and R. J. Elliott, J. Phys. C 4, 2370 (1971).
• (32) M. S. L. du Croo de Jongh and J. M. J. van Leeuwen, Phys. Rev. B 57, 8494 (1998).
• (33) M. Hasenbusch, K. Pinn, and S. Vinti, Phys. Rev. B 59, 11471 (1999).
• (34) C.-W. Liu, A. Polkovnikov, and A. W. Sandvik (unpublished).
• (35) G. E. Santoro and E. Tosatti, J. Phys. A: Math. Gen. 39 R393 (2006).
• (36) V. Bapst, L. Foini, F. Krzakala, G. Semerjian, and F. Zamponi, Phys. Rep. 523 127 (2013).
• (37) S. Baroni and S. Moroni, Phys. Rev. Lett. 82 4745 (1999),
• (38) G. Carleo, F. Becca, S. Moroni, and S. Baroni, Phys. Re. E 82 046710 (2010).
• (39) M. B. Hastings, M. H. Freedman, arXiv:1302.5733.