Universal short time quantum critical dynamics of finite size systems

TL;DR

This paper studies the universal short-time quantum critical dynamics in finite systems using imaginary time relaxation and quantum Monte Carlo, revealing distinct critical exponents and locating quantum critical points in 1D and 2D transverse field Ising models.

Contribution

The paper introduces a method to determine quantum critical points and exponents from short-time imaginary time dynamics in finite quantum systems, validated by simulations on 1D and 2D Ising models.

Findings

Quantum critical points accurately located in 1D and 2D models.

Critical initial slip exponents differ from classical counterparts.

Universal scaling behaviors observed in short-time quantum dynamics.

Abstract

We investigate the short time quantum critical dynamics in the imaginary time relaxation processes of finite size systems. Universal scaling behaviors exist in the imaginary time evolution and in particular, the system undergoes a critical initial slip stage characterized by an exponent $\theta$, in which an initial power-law increase emerges in the imaginary time correlation function when the initial state has zero order parameter and vanishing correlation length. Under different initial conditions, the quantum critical point and critical exponents can be determined from the universal scaling behaviors. We apply the method to the one- and two-dimensional transverse field Ising models using quantum Monte Carlo simulations. In the one-dimensional case, we locate the quantum critical point at $(h/J)_{c}=1.00003(8)$ \thirdrevise{in the thermodynamic limit}, and estimate the critical initial slip exponent $\theta=0.3734(2)$, static exponent $\beta/\nu=0.1251(2)$ \thirdrevise{by analyzing data on chains of length $L=32\sim 256$ and $L=48\sim 256$, respectively}. For the two-dimensional square-lattice system, the critical coupling ratio is given by $3.04451(7)$ \thirdrevise{in the thermodynamic limit} while the critical exponents are \thirdrevise{$\theta=0.209(4)$ and $\beta/\nu=0.518(1)$ estimated by data on systems of size $L=24\sim 64$ and $L=32\sim 64$, correspondingly.} Remarkably, the critical initial slip exponents obtained in both models are notably distinct from their classical counterparts, owing to the essential differences between classical and quantum dynamics. The short time critical dynamics and the imaginary time relaxation QMC approach can be readily adapted to various models.

Figures

Equations

$$\frac{\partial }{\partial \tau }|\psi (\tau )⟩=-H|\psi (\tau )⟩,$$

$$|\psi (\tau )⟩=Z_0^{-1}{e}^{-\tau H}|\psi ({\tau }_0)⟩,$$

$$Z_0=||e^{-\tau H}|\psi ({\tau }_0)⟩||,$$

$$|\psi (\tau )⟩$$

$$P(\tau ,g,M_0,L)=b^{\varphi }P({\tau }^{\prime },g^{\prime },M_0^{\prime },L^{\prime }),$$

$$M^k(\tau ,g,M_0,L)=b^{-k\beta /\nu }{M}^k(b^{-z}\tau ,b^{1/\nu }g,b^{x_0}{M}_0,b^{-1}L).$$

$$M^k(\tau ,M_0,L)={\tau }^{-k\beta /\nu z}{f}_M({\tau }^{x_0/z}{M}_0,{\tau }^{-1/z}L),$$

$$M^2(\tau ,L)=L^{-2\beta /\nu }{f}_M(\tau /L^z),$$

$$C(\tau )\equiv \underset{M_0\to 0}{\lim }\frac{M(\tau )}{M_0},$$

$$C(\tau ,L)=L^{\theta z}{f}_C(\tau /L^z),$$

$$Q(\tau )\equiv ⟨\text{sgn}(M(\tau ))⟩$$

$$Q(\tau ,g,L)=f_Q(\tau /L^z,L^{1/\nu }g),$$

$$A(L)=a{L}^{\sigma }(1+b{L}^{-\omega })$$

$$Z=⟨\Psi (\tau )⟩\Psi (\tau )=⟨\psi ({\tau }_0)⟩e^{-\tau H}{e}^{-\tau H}\psi ({\tau }_0),$$

$$H=-\sum _{b=1}^{N_{\text{op}}}{H}_b,$$

$$Z=\sum _n^{\infty }\sum _{S_n}⟨\psi ({\tau }_0)⟩\frac{{\beta }^n}{n!}{S}_n\psi ({\tau }_0),$$

$$Z=\sum _n^{\infty }\sum _{S_n}\frac{{\beta }^n}{n!}⟨\psi ({\tau }_0)⟩\sum _n^{k}w(n,k)S_n^k\psi ({\tau }_0),$$

$$\sum _{n_1=0}^{n}w(n,n_1)=1,\forall n.$$

$$H=-J\sum _{⟨i,j⟩}^{L^d}{\sigma }_i^z{\sigma }_j^z-h\sum _i^{L^d}{\sigma }_i^x,$$

$$\theta$$

$$\beta /\nu$$

$$\theta$$

$$\beta /\nu$$

$$|\psi ({\tau }_0)⟩=\sum _z\sqrt{f(M_0)}|z⟩,$$

$$f(M_0)=\prod _j\frac{1}{2}(1+M_0{\sigma }_j^z),$$

$$f(M_0)=\frac{1}{2^N}(1+M_0\sum _j{\sigma }_j^z).$$

$$Z=\frac{1}{2^N}\sum _{z_1,z_2}⟨z_1⟩e^{-\tau H}\sqrt{f(M_0)f(M_0)}{e}^{-\tau H}{z}_2,$$

$$Z=\frac{1}{2^N}\sum _{z_1,z_2}⟨z_1⟩e^{-\tau H}(1+M_0\overline{M})e^{-\tau H}{z}_2,$$

$$M(\tau )$$

$$\underset{M_0\to 0}{\lim }\frac{M(\tau )}{M_0}$$

$$A(\Delta t)=\frac{⟨O_{i+\Delta t}{O}_i⟩-{⟨O_i⟩}^2}{⟨O_i^2⟩-{⟨O_i⟩}^2},$$

Full text

I Introduction

The understanding of non-equilibrium dynamics in interacting many-body quantum systems is a crucial issue in modern physics with increasing focus stimulated by the experimental developments in the field of cold atoms Bloch et al. (2008). Different approaches can be used to take a system out of equilibrium, such as applying a driving field or pumping energy into the system, among which one simple but interesting candidate is to impose a sudden quench on an equilibrium system Dziarmaga (2010); Polkovnikov et al. (2011). After the sudden quench, the system may go through different types of non-equilibrium processes, for instance, it can gradually relax back to equilibrium or enter a quasi-steady prethermal state and end up in a thermalized state Berges et al. (2004), depending on the setup of the system and the dynamics that governs the evolution. Experimental interests in the quantum quench dynamics are heating up Sadler et al. (2006); Ritter et al. (2007); Trotzky et al. (2012); Navon et al. (2015); Nicklas et al. (2015) and many theoretical efforts have also been dedicated, some of which focus on the postquench long time universal scaling behaviors Lamacraft (2007); Rossini et al. (2009); Campos Venuti and Zanardi (2010); De Grandi et al. (2010); Deng et al. (2011); Foini et al. (2012); Iyer and Andrei (2012); Karrasch et al. (2012); Karl et al. (2017) while others pay their attention to the transient dynamics in the short time regime after the quench Dalla Torre et al. (2013); Yin et al. (2014); Gagel et al. (2014, 2015); Calabrese et al. (2012a, b); Maraga et al. (2015); Chiocchetta et al. (2015, 2016, 2017); Liu and Täuber (2016).

In the late 1980’s, Janssen et al Janssen et al. (1989) and Huse Huse (1989) discovered respectively the universal short time critical dynamics (STCD) in classical phase transitions. Universal scaling behaviors are found in the relaxation process after performing a sudden quench to the critical point from high temperature with a small order parameter $M_{0}$ and vanishing correlation length Janssen et al. (1989); Li et al. (1995, 1996). The emergence of the initial power-law increase in the order parameter, characterized as the critical initial slip with an independent exponent $\theta$, is an important feature of the STCD. In the short time regime, after a transient microscopic period, the order parameter increases as $M(t)\propto M_{0}t^{\theta}$. The characteristic time scale of this short time regime depends on $M_{0}$ and scales as $t_{\text{cr}}\sim M_{0}^{-z/x_{0}}$, in which $x_{0}=\theta z+\beta/\nu$ is the scaling dimension of $M_{0}$ Janssen et al. (1989); Diehl and Ritschel (1993). After $t_{\text{cr}}$, the initial condition becomes irrelevant and the behavior of $M(t)$ crosses over to the well-known long time power-law decay $M(t)\sim t^{-\beta/\nu z}$ Janssen et al. (1989). In the past few decades, the STCD has become a powerful tool in studying critical properties Li et al. (1994, 1995, 1996); Zheng (1998) and successfully applied to various models Ying et al. (2001); Albano et al. (2011); Tomé and de Oliveira (1998); Zheng (1996); Zheng et al. (2003); da Silva et al. (2002); Arashiro and Drugowich de Felício (2003); Baumann and Gambassi (2007); da Silva et al. (2013, 2014); Schüke and Zheng (1995); Luo and Zheng (1995); Grassberger (1995).

The classical STCD has been studied extensively (see Ref. Albano et al. (2011) for a review) while the same issue in quantum systems is attracting increasing attention in recent years. Experimentally, temporal scaling crossover is observed in the time evolution of an isolated one-dimensional (1d) two-component Bose gas, when quenched into the vicinity of the critical point Nicklas et al. (2015). The critical initial slip and crossover to long time non-critical regime are found in the correlation length of the order parameter Nicklas et al. (2015). Theoretical interests focus on several aspects, including open and isolated quantum systems as well as real and imaginary time evolutions Calabrese et al. (2012a, b); Gagel et al. (2014, 2015); Maraga et al. (2015); Chiocchetta et al. (2015, 2016, 2017); Liu and Täuber (2016); Yin et al. (2014); Zhang et al. (2014); Weinberg and Sandvik (2017). The works in Refs. Gagel et al. (2014, 2015) present investigations on an $N$-component $\varphi^{4}$ model coupled to an external bath, where the critical initial slip and long time power-law decay are found in the real time relaxation of quantum open systems. A renormalization group analysis by Chiocchetta et al Chiocchetta et al. (2017) shows the existence of the critical initial slip and dynamical crossovers between different scaling regimes in prethermal states. Depending on dimension and the energy scale injected by the quench, crossover can exist between two critical initial slip stages with different exponents $\theta$’s governed by a quantum and a classical prethermal fixed point correspondingly Chiocchetta et al. (2017). Due to the unitary nature of real time evolution, only the prethermal stage has similar evolution behaviors as the classical case Maraga et al. (2015); Chiocchetta et al. (2015, 2016, 2017). Differently, in imaginary time, high energy excited states decay so fast that the evolution is governed by low energy levels close to the ground state and exhibits universal power-laws Yin et al. (2014); Zhang et al. (2014). A scaling theory for the quantum STCD in imaginary time was proposed in Ref. Yin et al. (2014), and realized in 1d Ising systems. The critical initial slip also appears in the imaginary time evolutions of quantum systems, i.e., the time dependence of the order parameter in Ref. Janssen et al. (1989) also holds in imaginary time $M(\tau)\propto M_{0}\tau^{\theta}$ for a small finite $M_{0}$ Yin et al. (2014). The critical initial slip terminates at time scale $\tau_{\text{cr}}\sim M_{0}^{-z/x_{0}}$, which is called the initial-value time Diehl and Ritschel (1993), and crosses over to the universal long time power-law decay $M(\tau)\sim\tau^{-\beta/\nu z}$. The value of the critical initial slip exponent is found to be $\theta=0.373$ Yin et al. (2014) in the 1d transverse field Ising model (TFIM), in sharp difference from its classical counterpart, the two-dimensional (2d) classical Ising model $\theta=0.191(1)$ Okano et al. (1997); Albano et al. (2011).

Even with the progress made by density-matrix renormalization group and related matrix-product state methods Schollwöck (2005, 2011); Daley et al. (2004), real time dynamics is difficult to deal with computationally for most of quantum systems, especially for higher dimensions. As demonstrated in Refs. De Grandi et al. (2011); De Grandi et al. (2013); Liu et al. (2013); Witczak-Krempa et al. (2014); Karrasch et al. (2012), real time and imaginary time dynamics bear so many similarities that some scaling properties of real time dynamics can be predicted by the imaginary time dynamics. Recent developments in quantum Monte Carlo (QMC) Liu et al. (2013); De Grandi et al. (2011); De Grandi et al. (2013) light up the possibility to simulate the imaginary time relaxation process in a wide range of “sign-problem” free quantum models, offering a bright way to study the quantum STCD.

In this work, we study the quantum STCD of finite size systems in imaginary time relaxation and apply it to both the 1d and 2d TFIMs using a ground state projection QMC approach. The quantum critical point and critical exponents ($\theta$ and $\beta/\nu$) are extracted from the universal scaling behaviors. The key result of our work is the capture of the critical initial slip of the 2d square-lattice TFIM. To our knowledge, the critical initial slip exponent $\theta$ of the 2d TFIM has not been computed yet. In the meantime, we obtain the $\theta$ of the 1d TFIM, soundly confirm the result found in Ref. Yin et al. (2014). In both cases, the critical initial slip exponents are distinct from their classical counterparts. Our results show the capability of the quantum STCD method in determining critical properties and shed light on future applications to other models.

The rest of the paper is organized in the following way. In Sec. II we first briefly review the imaginary time dynamics. In Sec. III we generalize the scaling theory for the quantum STCD to finite size systems in imaginary time. The key idea of the projector QMC method is outlined in Sec. IV while Sec. V introduces the properties of the models studied. Numerical results of the 1d and 2d TFIMs are presented in Sec. VI and a summary is given in Sec. VII.

II Imaginary time dynamics

In this section, we discuss the imaginary time evolution of a quantum state $\ket{\psi(\tau)}$. By performing the standard Wick rotation of $t\rightarrow-i\tau$, the Schrödinger equation describes the imaginary time evolution of the wavefunction as Zinn-Justin (1996); Altland and Simons (2006)

[TABLE]

with the Planck constant being $1$ and normalization condition $\braket{\psi(\tau)}{\psi(\tau)}=1$. A formal solution of Eq. (1) is given by

[TABLE]

in which $e^{-\tau H}$ is the imaginary time evolution operator and the normalization factor $Z_{0}$ is

[TABLE]

where $||...||$ denotes a modulo operation.

In the energy representation, Eq. (2) reads

[TABLE]

where the coefficients are given by the overlap of the initial state and the $i$th eigenstate of the Hamiltonian $c_{i}=\braket{i}{\psi(\tau_{0})}$ and $\Delta=E_{1}-E_{0}$ is the gap between the first excited state and the ground state. It is indicated in Eq. (4) that the high energy levels decay too fast to be non-negligible and $\ket{\psi(\tau)}$ is governed by the low energy levels during the imaginary time evolution Yin et al. (2014).

As pointed out in Ref. Zhang et al. (2014), a dissipative equation of the probability to find a given eigenstate can be derived from the Schrödinger equation and it is argued that this dissipative equation in fact exhibits similar evolution properties with the classical master equation, especially the critical initial slip in the short time stage, though, the dynamics are essentially different Zhang et al. (2014).

III Quantum short time critical dynamics in finite size systems

In this section, we discuss the quantum STCD in finite size systems during the imaginary time evolution under different quench protocols.

When a system is initially off-critical with vanishing correlation length, a sudden quench to the critical point triggers the relaxation Janssen et al. (1989); Yin et al. (2014). In imaginary time, universal scaling behaviors are found after a short period of non-universal microscopic time scale $\tau_{\text{mic}}$ Yin et al. (2014). Taking the system size into account, the scaling form of an observable $P$ is given by Yin et al. (2014)

[TABLE]

where $\tau^{\prime}=b^{-z}\tau$, $g^{\prime}=b^{1/\nu}g$, $M_{0}^{\prime}=b^{x_{0}}M_{0}$, and $L^{\prime}=b^{-1}L$ for an arbitrary scaling factor $b$. The arguments $\tau$, $g$, $M_{0}$ and $L$ refer to the imaginary time, distance to the critical point, initial value of the order parameter and system size, respectively. $x_{0}$ is the scaling dimension of $M_{0}$, satisfying $x_{0}=\theta z+\beta/\nu$ as in the classical case Janssen et al. (1989); Diehl and Ritschel (1993); Yin et al. (2014), and $\phi$ is related to the critical exponent of the quantity studied. We consider two initial conditions $M_{0}=0$ and $M_{0}=M_{\text{sat}}$ ($M_{\text{sat}}$ being the saturate value of the order parameter) since they are both fixed points of the scaling transformation.

With $\phi=-\beta/\nu$, the $k$th moment of the order parameter follows Janssen et al. (1989); Li et al. (1995, 1996); Yin et al. (2014)

[TABLE]

At the critical point $g=0$, setting the scaling factor $b=\tau^{1/z}$, one reaches

[TABLE]

where $f_{M}$ is a scaling function related to $M^{k}$ (similar notations apply in the context). We focus on $k=1$ (corresponding to the order parameter) and $k=2$. In infinite systems, when $M_{0}$ is small but finite, the critical initial slip $M(\tau)\propto M_{0}\tau^{\theta}$ appears by expanding $\tau^{x_{0}/z}M_{0}$ in the short time regime $\tau\ll M_{0}^{-z/x_{0}}$ Yin et al. (2014). The presence of finite system size in Eq. (7) implies that the relaxation involves another characteristic time scale $L^{z}$, which is known as the finite-size relaxation time, that controls the scale of the short time regime Diehl and Ritschel (1993); Li et al. (1994).

When $M_{0}=0$ and $g=0$, $M^{2}(\tau,L)$ obeys

[TABLE]

if one choose $b=L$. Only the finite-size relaxation time $L^{z}$ is involved in the scaling form of $M^{2}(\tau,L)$ since the initial-value time $M_{0}^{-z/x_{0}}$ diverges as $M_{0}\rightarrow 0$.

In analogy to the classical case Huse (1989); Tomé and de Oliveira (1998), we consider

[TABLE]

which is called the imaginary time correlation function since it measures the correlation of the order parameter between the initial state and the state at $\tau$ as indicated in Appendix A.

Combining Eqs. (6) and (9), the scaling form is given by

[TABLE]

with the scaling factor $b=L$. In infinite systems, as indicated in Eq. (9), the behavior of $C(\tau)$ is characterized by the critical initial slip exponent as well, $C(\tau)\propto\tau^{\theta}$ Tomé and de Oliveira (1998), while the order parameter is incapable of showing the critical initial slip when $M_{0}=0$, since symmetry restricts $M(\tau)=0$ Janssen et al. (1989); Diehl and Ritschel (1993); Li et al. (1994); Yin et al. (2014). In the presence of finite system size, $C(\tau,L)$ shows power-law increase with an exponent $\theta z$ for fixed $\tau/L^{z}$.

When $M_{0}=M_{\text{sat}}$ and $g$ varies, let us consider a quantity $Q$ related to the average sign of the order parameter, defined as de Oliveira (1992); Silvério Soares et al. (1997)

[TABLE]

in which sgn is the sign function and $\braket{...}$ denotes non-equilibrium average. $Q(\tau)$ records the system’s memory of the initial condition during the relaxation, which decays from the initial value $Q_{0}=1$ to the equilibrium value. As argued in Refs. de Oliveira (1992); Silvério Soares et al. (1997), in Ising systems, $\phi$ is [math], leading to the scaling form

[TABLE]

with $b=L$. It is clear that for fixed $\tau/L^{z}$, different system sizes undergo curve crossing of $Q(\tau,g,L)$ at $g=0$, in analogy to the Binder cumulants Binder (1981a); Li et al. (1995, 1996); Landau and Binder (2009).

The scaling behaviors of $M^{2}(\tau,L)$, $C(\tau,L)$ and $Q(\tau,g,L)$ are useful in detecting the quantum critical point and critical exponents.

(i) With $M_{0}=M_{\text{sat}}$, the scaling form of $Q$ indicates that for fixed $\tau/L^{z}$, a series of size-dependent critical points can be extracted by performing curve crossing analysis for different system sizes Landau and Binder (2009). The quantum critical point is therefore obtained by extrapolating to $L\rightarrow\infty$ Binder (1981a); Landau and Binder (2009).

(ii) With $M_{0}=0$, when $\tau/L^{z}$ is fixed, $M^{2}(L)$ and $C(L)$ have power-law behaviors, giving the critical exponents $-2\beta/\nu$ and $\theta z$. The requirements of small $M_{0}$ and extrapolation to $M_{0}\rightarrow 0$ to compute $\theta$ in previous investigations Grassberger (1995); Okano et al. (1997); Zheng (1998); Ying et al. (2001); Albano et al. (2011) is therefore avoided. Note that Eq. (8) also holds when $M_{0}=M_{\text{sat}}$ and $g=0$, meaning that $\beta/\nu$ can be measured in the long time regime ($\tau_{\text{cr}}\sim 1$) under this initial condition as well Yin et al. (2014).

(iii) In order to get reliable estimates of the critical exponents, the leading scaling correction should be taken into account. We use the following ansatz Pelissetto and Vicari (2002); Luo et al. (2001); Albano et al. (2011)

[TABLE]

to extract the critical exponents and the leading correction exponent $\omega$. In Eq. (13), $A(L)$ corresponds to $M^{2}(L)$ or $C(L)$ for fixed $\tau/L^{z}$ and $\sigma$ equals to $-2\beta/\nu$ and $\theta z$ accordingly. The choice of $\tau/L^{z}$ should not affect the asymptotic behavior of $A(L)$ but only bring different finite size corrections.

The quantum STCD provides access to the critical properties while the system is relaxing towards its ground state, and it overcomes the critical slowing down problem (in the sense of physical time) Li et al. (1995, 1996); Yin et al. (2014), requiring much less computational efforts comparing with the traditional finite size scaling studies.

IV Quantum Monte Carlo Method

The method we employ is the projector QMC algorithm based on the stochastic series expansion (SSE), which allows us to obtain the relaxation properties of a system after evolving an imaginary time $\tau$ De Grandi et al. (2011); Farhi et al. (2012); Liu et al. (2013); De Grandi et al. (2013). The idea of the QMC algorithm is to substitute the imaginary time evolution operator $e^{-\tau H}$ with its Taylor expansion series and perform importance samplings in the configuration space, which contains information needed to compute the expectation values of physical quantities.

Following Eqs. (1) and (2), the partition function $Z$ is defined as

[TABLE]

in which $\ket{\Psi(\tau)}=e^{-\tau H}\ket{\psi(\tau_{0})}$. Divide the Hamiltonian into a sum of bond operators

[TABLE]

and series expand the exponential factor $e^{-\tau H}$ with $\beta=2\tau$, Eq. (14) becomes Farhi et al. (2012); Sandvik (2010)

[TABLE]

where $S_{n}$ denotes the operator sequence of bond operators.

In order to calculate precise expectation values at imaginary time $\tau$, a binomial factor $w(n,k)=\binom{n}{k}/2^{n}$ should be inserted to Eq. (16) Farhi et al. (2012); Sandvik (2010)

[TABLE]

in which $S_{n}^{k}$ is an operator sequence of length $n$ with an imaginary “cut” at position $k$ running through it and $w(n,k)$ satisfies

[TABLE]

In fact, $w(n,k)$ comes out immediately if one expand the two exponential factors independently, concatenate and sample the operator sequences together. $w(n,k)$ describes contributions of different propagations along the imaginary time direction and plays an important role when it aims at computing precise results of the states after small length of evolutions Farhi et al. (2012). As the operator sequence grows longer, the influence of $w(n,k)$ becomes smaller and the partition function reduces to Eq. (16), suggesting that measurements can be taken in the middle of the operator sequence.

The expansion order $n$ can be truncated to some maximum order $n_{\text{max}}$ without causing detectable errors Sandvik (2010). The expectation value of an observable $O$ is given by $\braket{O(\tau)}=\braket{\Psi(\tau)}{O}{\Psi(\tau)}/Z$ as in the standard SSE algorithm. In the presence of $w(n,k)$, measurements should be taken at different positions of the operator sequence, bringing different contributions to $\braket{O(\tau)}$ Farhi et al. (2012).

V Model

The models we study here are the transverse field Ising models, defined by the Hamiltonian

[TABLE]

in which $\braket{i,j}$ are nearest-neighbor sites and $\sigma^{x}$ and $\sigma^{z}$ are the Pauli matrice. The couplings $J$ are chosen to be ferromagnetic interactions ($J>0$) and $h$ is a external field in the transverse direction. The order parameter is given by the magnetization $M=\braket{\sigma^{z}}$ and periodic boundary condition is used.

The 1d TFIM can be solved rigorously by mapping into the 2d classical Ising model Suzuki (1976). The critical point of the ordered/disordered quantum phase transition is $(h/J)_{c}=1$. The static exponents of interest here are $\beta=1/8$ and $\nu=1$ and the dynamic exponent is $z=1$ Sachdev (1999); Sondhi et al. (1997); Vojta (2003). Application of the STCD to the 2d classical Ising model has attracted many numerical attention Grassberger (1995); Okano et al. (1997); Zheng (1998); Ying et al. (2001); Albano et al. (2011). The critical initial slip arises in the short time regime with an exponent $\theta=0.191(1)$ Okano et al. (1997); Albano et al. (2011). For the 1d TFIM, Ref. Yin et al. (2014) extracts the critical initial slip exponent from the imaginary time relaxation, giving $\theta=0.373$, which is almost twice as large as its classical counterpart.

Owing to the fact that the 2d TFIM has no rigorous solution, many numerical efforts have been devoted to study its properties, including locating the critical point with high precision Hamer (2000); Liu et al. (2013); Blöte and Deng (2002); du Croo de Jongh and van Leeuwen (1998) and extracting critical exponents Pfeuty and Elliott (1971); du Croo de Jongh and van Leeuwen (1998); Albuquerque et al. (2010). For the 2d square-lattice TFIM, an exact diagonalization study on a $6\times 6$ lattice Hamer (2000) gives the critical point $(h/J)_{c}=3.04497(18)$ and a recent work using quasi-adiabatic QMC claims a higher precision estimate $(h/J)_{c}=3.04458(7)$ Liu et al. (2013). The dynamic exponent has an exact result $z=1$ Sachdev (1999); Sondhi et al. (1997); Vojta (2003), while the static exponents are shared with the 3d classical Ising model since they belong to the same universality class. A recent study of the 3d Ising universality class estimates $\beta/\nu=0.51814(5)$, which is not directly reported in Ref. Hasenbusch (2010) but can be obtained using the scaling relation $\beta/\nu=(1+\eta)/2$ with $\eta=0.03627(10)$ Hasenbusch (2010). In Ref. Jaster et al. (1999), an STCD study of the 3d classical Ising model gives $\beta/\nu=0.517(2)$ along with the critical initial slip exponent $\theta=0.108(2)$. To our knowledge, the critical initial slip of the 2d TFIM has not been studied yet.

VI Numerical results

In this section, we present results of the 1d and 2d TFIMs. Two different initial conditions $M_{0}$ both with vanishing correlation length are employed here. We shall first use $M_{0}=M_{\text{sat}}=1$, which is a state with all spins in the up direction $\ket{\Uparrow}$. Quenches from this ordered state to the vicinity of the critical point are imposed so that the critical coupling $h_{c}$ (for simplicity, $J=1$) can be detected through the finite size scaling analysis of $Q(\tau,g,L)$. The other initial condition is a state with all spins aligned along the transverse direction $\ket{\Rightarrow}$ with $M_{0}=0$. The critical initial slip can be observed during the relaxation and the exponents $\theta$, $\beta/\nu$ are extracted from the scaling behaviors of $C(\tau,L)$ and $M^{2}(\tau,L)$. The exact value of $z=1$ Sachdev (1999); Sondhi et al. (1997); Vojta (2003) for both models is used in the scaling relations.

In the QMC simulations, both local and cluster updates are carried out to perform efficient samplings Sandvik (2003). The computational time scales as $L^{d}\tau$ Sandvik (2003) and therefore proportional to $L^{d+z}$ (as the ratio $\tau/L^{z}$ is fixed), but still, large amount of efforts required to achieve the ground state, corresponding to much larger $\tau/L^{z}$, is saved by the quantum STCD. The convergence and autocorrelations of the QMC algorithm are discussed in Appendix B.

VI.1 1d transverse field Ising model

The 1d TFIM provides a rigorous testing ground for the QMC algorithm and the scaling for the quantum STCD in finite size systems.

When the system is quenched from $\ket{\Uparrow}$ to the vicinity of the critical point, the scaling behavior of $Q(\tau,g,L)$ is described by Eq. (12). We investigate $Q$’s dependence on the tunable parameter $h$ for different system sizes, as shown in the upper panel of Fig. 1. The argument $g$ in Eq. (12) is replaced by $h-h_{c}$ in order to carry out curve crossing analysis. We fixed the evolution time at $\tau=L/4$ for $L$ from $8$ to $128$ and extract the crossing points of system size pairs of $L$ and $2L$ Landau and Binder (2009) by fitting $Q(h)$ using polynomial forms (up to cubic terms). Panel (b) of Fig. 1 shows the size dependence of the crossing points. By fitting to the form of $h_{c}(L)=h_{c}+aL^{-\omega}$ Binder (1981b) and extrapolating to $L\rightarrow\infty$, we obtain the critical point $h_{c}=1.00003(8)$ with $\omega=2.33(1)$, in good agreement with the exact result $h_{c}=1$ Sachdev (1999); Vojta (2003); Sondhi et al. (1997).

The fittings are performed by searching for parameters to have the goodness of fit $\chi^{2}$ per degree of freedom $N_{\mathrm{DOF}}$ close to $1$. To estimate the statistical errors of the fitting parameters, we generate different sets of data by adding Gaussian noise to the data, with the magnitude of the noise equal to the error bars of the data Sandvik (2010). Repeating the fits for a large number of times, the standard deviation of the distributions give the error bars of the fitting parameters Sandvik (2010). The same way of error estimating is used in the following results.

Next we prepare the system in the state $\ket{\Rightarrow}$ and suddenly quench it to the critical point $h_{c}=1$ at $\tau=0$. After a short microscopic period, the critical initial slip starts to emerge. We compute the imaginary time correlation function $C(\tau,L)$ according to Appendix A in QMC simulations. As illustrated in Fig. 2, for $\tau=L/8$, $C(L)$ behaves as a power-law form with the critical initial slip exponent $\theta$. We perform power-law fitting using the form of $C(L)=aL^{\theta}$ to extract $\theta$. In Fig. 2(b), we plot $\theta$ as a function of $1/L_{\mathrm{min}}$ ($L_{\mathrm{min}}$ being the minimal system size included into the fit) to show the dependence on the fitting range. Each point is obtained by performing power-law fitting with/without correction in the range from $L_{\mathrm{min}}$ to $L_{\mathrm{min}+p}$. When there is no correction included, $p=4$ is used so as to satisfy the minimum requirement of $N_{\mathrm{DOF}}=1$ while in the presence of correction, $p$ is set to $6$. As can be seen in Fig. 2(b), $\theta$ only fluctuates slightly even for small system sizes. Therefore, we use the same form to fit the data within the range from $L=32$ to $256$, yielding the critical initial slip exponent $\theta=0.3734(2)$, in excellent agreement with Ref. Yin et al. (2014). No finite size correction is included since $C(L)$ is well described by the power-law form except two smallest system sizes. We also read off the the result $\theta=0.3733(7)$ given by the fitting range that consists of the largest system sizes $L=128$ to $256$ to provide knowledge about the magnitudes of statistical errors involved in Fig. 2(b).

The value of $\theta$ of the 1d TFIM is in great difference from the result of the 2d classical Ising model $\theta=0.191(1)$ Okano et al. (1997); Albano et al. (2011).

We also include $C(L)$ of $\tau=L/16$ in Fig. 2 to show the influence of different time-space ratios. Finite size correction is needed to produce good fits. As shown in Fig. 2(b), the dependence on fitting range is more obvious than the case $\tau/L=1/8$. Even though $\theta$ of $L_{\mathrm{min}}=40$ and $56$ are already close, to be on the safe side, we only use the range starting from the largest $L_{\mathrm{min}}$ (i.e. $L=56\sim 128$) to represent the result of $\tau/L=1/16$, which reads $\theta=0.3733(6)$, agreeing with the result of $\tau/L=1/8$ within errorbar. The choice of $\tau/L$ does not affect the asymptotic behavior but only cause different finite size corrections. Depending on the quantity studied, optimal ratio(s) that alleviates finite size correction may exist. For $C(L)$, $\tau/L=1/8$ is close to the optimal value so that there is little finite size effect remained.

The static critical exponent $\beta/\nu$ can be measured under the same initial condition, as shown in Fig. 3. For $\tau=L/8$, we compare the result in the presence and absence of correction, both with $N_{\mathrm{DOF}}=1$, as illustrated in Fig. 3(b). The exponent $\beta/\nu$ only fluctuates slightly when finite size correction is included while in the absence of correction, $\beta/\nu$ changes with varying fitting range. Using Eq. (13), the largest system sizes starting from $L_{\mathrm{min}}=96$ gives $\beta/\nu=0.1249(8)$. Since $\beta/\nu$ agrees within errorbar when $L_{\mathrm{min}}\geq 48$, we expand the range to $L=48\sim 256$ and obtain $\beta/\nu=0.1251(2)$, in good agreement with the exact result $\beta/\nu=1/8$ Sachdev (1999); Vojta (2003); Sondhi et al. (1997).

Likewise, we plot $M^{2}(L)$ for $\tau/L=1/16$ in Fig. 3(a) and the dependence of $\beta/\nu$ on fitting range in panel (b). Even in the presence of correction, $\beta/\nu$ only becomes close when $L_{\mathrm{min}}\geq 48$. To be safe, we use the result of the largest $L_{\mathrm{min}}=56$ to represent the final result $\beta/\nu=0.1253(5)$, consistent with the result of $\tau/L=1/8$ and the exact value as well. Figure 3(b) also implies that to get more reliable estimate, larger systems should be involved.

We therefore arrive at the final estimates of $\theta$ and $\beta/\nu$ for the 1d TFIM

[TABLE]

based on analysis on data of $\tau/L=1/8$ within the fitting range $L=32\sim 256$ and $L=48\sim 256$, respectively. Note that the errors here for the exponents $\theta$, $\beta/\nu$(as well as $\omega$, $a$ and $b$) are only statistical variations among the data points in a given range.

Our results agree nicely with the rigorous solutions Sondhi et al. (1997); Sachdev (1999); Vojta (2003) and the numerical results Yin et al. (2014), demonstrating the validity of the quantum STCD and the QMC algorithm in detecting critical properties.

VI.2 2d transverse field Ising model

In the following, we apply the quantum STCD to the 2d TFIM.

As in the 1d case, we shall first determine the critical point by carrying out curve crossing analysis of $Q(h,L)$ with $\tau=L/4$ and $M_{0}=1$, as illustrated in Fig. 4. Crossing points of different system size pairs of $L$ and $2L$ are extracted using polynomial fits (up to cubic terms) of the $Q(h,L)$ Landau and Binder (2009). Using the form of $h_{c}(L)=h_{c}+aL^{-\omega}$ Binder (1981a), we estimate the critical coupling $h_{c}=3.04451(7)$, in good agreement with the result in Ref. Liu et al. (2013), in which the calculations are carried out on systems of similar sizes to those in this study.

When the system is suddenly quenched from the state $\ket{\Rightarrow}$ to the critical point $h_{c}=3.04451$, as expected, the initial increase of $C(\tau,L)$ arises after a transient microscopic period. In Fig. 5, we plot $C(L)$ for $\tau=L/4$ and compare the exponent $\theta$ obtained using different fitting forms. When there is no correction, the power-law form only captures the behavior of $C(L)$ for large systems and small system sizes should be discarded in order to get good fits. With the leading correction, the results are improved and become converged (within errorbar) for large $L_{\mathrm{min}}$. For comparison, we marks the result $\theta=0.210(7)$ of $L_{\mathrm{min}}=44$ in the presence of correction using a dashed line in Fig. 5(b). Since $\theta$ converges when $L_{\mathrm{min}}\geq 24$, we use Eq. (13) to fit data in the range $L=24\sim 64$, yielding $\theta=0.209(4)$, along with the leading correction exponent $\omega=2.25(6)$. The value of $\theta$ is again distinct from its classical counterpart, the 3d classical Ising model, $\theta=0.108(2)$ Jaster et al. (1999).

In both 1d and 2d, the quantum models have a critical initial slip exponent approximately twice as large as their classical counterparts. In fact, there is no reason to expect the critical initial slip exponent are in the $d$ to $d+z$ classical/quantum correspondence. Though both dissipative, the dynamics governing the two cases are essentially different. The quantum systems are described by the Schrödinger’s equation while the classical systems follow the master equation, with different dynamic exponents $z$ Yin et al. (2014); Zhang et al. (2014). Specifically, the dynamic exponent of the 1d and 2d TFIMs both equal to $1$ Sachdev (1999); Sondhi et al. (1997); Vojta (2003) while in the 2d and 3d classical Ising models, the Metropolis dynamics gives $z$ equals to $2.1667(5)$ Nightingale and Blöte (2000) and $2.042(6)$ Jaster et al. (1999), respectively. The critical initial slip exponent $\theta$ is dependent on the dynamical equation Yin et al. (2014); Zhang et al. (2014) and different dynamic exponents $z$ can lead to different values of $\theta$. However, whether or not the similar ratios of the $\theta$’s between the quantum models and their classical counterparts are coincidental may need further investigations.

Next, we measure the static critical exponent $\beta/\nu$ in the same way, as seen in Fig. 6. The initial condition and evolution time are the same with the measurement of $C(L)$. Panle (b) of Fig. 6 compares the exponent $\beta/\nu$ obtained using different fitting forms. In the absence of correction, $\beta/\nu$ strongly depends on the fitting range. Including the leading correction produces steady results when $L_{\mathrm{min}}\geq 32$, among which, we read off the exponent $\beta/\nu=0.518(2)$ given by the largest $L_{\mathrm{min}}=44$. Expanding the fitting range to $L=32\sim 64$, we obtained $\beta/\nu=0.518(1)$, agreeing well with with the 3d classical Ising model result $\beta/\nu=0.51814(5)$ Hasenbusch (2010).

We reach our final estimates of $\theta$ and $\beta/\nu$ for the 2d TFIM

[TABLE]

using the results in the range $L=24\sim 64$ and $L=32\sim 64$, correspondingly. Again, the errors only represent statistical variations in the fitting range.

As can be seen in Figs. 5 and 6, with the leading finite size correction, small system sizes down to $L=24$ and $32$ (for $C(L)$ and $M^{2}(L)$, respectively) can be included in the fitting function but subleading finite size corrections are still needed in order to take all sizes into consideration.

VII Summary

In summary, we studied the quantum STCD in imaginary relaxation of finite size systems. By imposing different quench protocols on equilibrium systems, critical properties can be detected from the universal scaling behaviors and the critical slowing down is avoided Diehl and Ritschel (1993); Li et al. (1994, 1995, 1996); Yin et al. (2014). The method is applied to the 1d and 2d TFIMs using QMC simulations.

We located the quantum critical points of the 1d and 2d TFIMs with high precision, which are $h_{c}=1.00003(8)$ and $3.04451(7)$, respectively. The result of the 1d TFIM agrees excellently with the exact result Sachdev (1999); Sondhi et al. (1997); Vojta (2003) and the 2d case is in good agreement with the critical coupling ratio given in Ref. Liu et al. (2013). The emergence of the critical initial slip are observed in both models. In 1d, finite size scaling analysis on data between $L=32\sim 256$ gives the critical initial slip exponent $\theta=0.3734(2)$, in excellent agreement with Ref. Yin et al. (2014), but approximately twice as large as its classical counterpart $\theta=0.191(1)$ Okano et al. (1997); Albano et al. (2011). In 2d, though predicted from the quantum STCD, no previous numerical observation is reported. Remarkably, our calculations are able to capture the critical initial slip exponent of the 2d TFIM $\theta=0.209(4)$ within the range $L=24\sim 64$ , which is again different from the value $0.108(2)$ of the 3d classical Ising model Jaster et al. (1999). Different critical initial exponents are actually expected in the dissipative relaxation processes of quantum and classical systems. Even though one can find equations that bear similar evolution properties in both cases, the dynamics are essentially different. In classical systems, the master equation governs the evolution after a sudden quench while in quantum systems, the Schrödinger equation is responsible for the imaginary time evolution Yin et al. (2014); Zhang et al. (2014). In addition, we obtained the static exponent $\beta/\nu=0.1251(2)$ of the 1d TFIM using data within $L=48$ and $256$, agree well with the exact result $1/8$ Sachdev (1999); Sondhi et al. (1997); Vojta (2003). For the 2d TFIM, we estimate $\beta/\nu$ to be $0.518(1)$ in the range from $L=32$ to $64$, in good agreement with a recent high-precision numerical estimate $0.51814(5)$ Hasenbusch (2010). It is also helpful to notice that, in the 2d case, the non-monotonic behavior of $\beta/\nu$ as a function of $L_{\mathrm{min}}$ indicates that there are potential corrections to the scaling forms introduced by short time effect in addition to finite size effect. In general, when $\tau/L^{z}\ll 1$, finite size effect plays a minor role while for $\tau/L^{z}\gg 1$, the scaling forms are reduced to the static case. In our calculations, the time-space ratio $\tau/L^{z}$ is fixed at values that both effects can come into play while we have only considered finite size corrections at present. Even though the short time effect can be taken into account by including leading, subleading finite size correction terms to the scaling forms (since $\tau/L^{z}$ is fixed) and we have shown that consistent asymptotic behaviors can already be obtained with only the leading finite size correction included, different ratios should indeed be considered in order to analyze the corrections caused by finite size and short time effect to the scaling forms in further studies. Our results manifest the capability of the quantum STCD to determine critical properties and indicate broader applications of the method to many other “sign-problem” free models.

acknowledgement

We would like to thank A. Sandvik, S. Capponi and J. Marino for their valuable discussions and critical reading this manuscript. We also thank an anonymous referee for his/her sugguestions on data analysis. YRS thanks P. Weinberg, L. Wang and N. Xu for their useful discussions. This project is supported by NKRDPC-2017YFA0206203, NSFC-11574404, NSFC-11275279, NSFG-2015A030313176, Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase), National Supercomputer Center in Guangzhou, and Leading Talent Program of Guangdong Special Projects.

Appendix A Measurement of the imaginary time correlation function

Here we discuss the measurement of the imaginary time correlation function $C(\tau)$ in QMC simulations.

In the standard $\sigma^{z}$ basis, let us consider a state with a given magnetization $M_{0}$ at $\tau_{0}$

[TABLE]

which a superposition of all $2^{N}$ ($N=L^{d}$ being the total particle number) basis states with different amplitudes. The coefficients related to $M_{0}$ are given by

[TABLE]

if one consider constructing a state site by site Tomé and de Oliveira (1998). For small $M_{0}$, up to linear terms in $M_{0}$

[TABLE]

The partition function in Eq. (14) becomes

[TABLE]

Plugging in Eq. (22), one arrives at

[TABLE]

where $\overline{M}$ equals to $[\sum_{i}{\sigma_{i}^{z}}+\sum_{j}{\sigma_{j}^{z}}]/2$. The indices $i$ and $j$ correspond to sum over the basis state $\ket{z_{1}}$ and $\ket{z_{2}}$, respectively. Therefore, the magnetization at time $\tau$ is given by

[TABLE]

in which $M=\sum_{i}{\sigma^{z}_{i}(\tau)}/N$. The first term on r.h.s vanishes due to up-down symmetry. Divide $M(\tau)$ by $M_{0}$ before taking the limit $M_{0}\rightarrow 0$, one gets

[TABLE]

where $\overline{M}$ is averaged over the two initial states $\ket{z_{1}}$ and $\ket{z_{2}}$ on the boundaries of sampling space in QMC simulations.

As indicated by the second line of Eq. (26), $C(\tau)$ measures the correlation between the initial state and the state at imaginary time $\tau$. $C(\tau)$ actually corresponds to $\braket{\sigma_{k}^{z}(0)\sigma_{-k}^{z}(\tau)}$ with momentum $k=0$ Huse (1989), where $\sigma_{k}^{z}(\tau)$ is the Fourier transform of the spins at $\tau$ and $\sigma_{k}^{z}(0)$ is averaged over $\ket{z_{1}}$ and $\ket{z_{2}}$ at $\tau=0$. For small $k$, the scaling form of $C(\tau)$ can be generalized by including the scale transformation of the momentum $k$, which reads $k^{\prime}=bk$. However, for large $k$, the fluctuation modes are far away from the low energy regime, hence, the behavior is not described by the scaling form. Therefore, we restrict our calculations to the $k=0$ sector and use a disordered initial state with $M_{0}=0$ and vanishing correlation length, as required by the quantum STCD for the scaling forms to apply. The critical initial slip exponent can be conveniently extracted without suffering from the requirements of small but finite magnetization and extrapolation to $M_{0}\rightarrow 0$.

Appendix B Convergence and autocorrelations

In this section, we present the convergence and the autocorrelations for the 1d TFIM to show the performance of the projector QMC algorithm.

To distinguish from the real time $t$ in the main text, we use $t_{\text{MC}}$ to refer to the Monte Carlo (MC) steps in the following. One MC step consists of a full sweep of single operator updates followed by constructions and updates of the clusters, which are flipped with probability $1/2$ Sandvik (2003). The details of sampling are described in Refs. Sandvik (2003, 2010).

As a convergence test, we compute the imaginary time correlation as a function of $t_{\text{MC}}$ for $\tau=4$ and $L=128$ with an initially disordered state (which is $\ket{\Rightarrow}$), shown in Fig. B.1. $C(t_{\text{MC}})$ converges to the equilibrium result rapidly after $t_{\text{MC}}\geq 40$. Each data point in Fig. B.1 is averaged over $10^{3}$ to $10^{4}$ independent simulations. The growth of $C(t_{\text{MC}})$ at short $t_{\text{MC}}$ corresponds to the increases of the cut-off expansion order $n_{\text{max}}$ (in Sec. IV) and the number of flipped clusters at the beginning of the simulations. As $t_{\text{MC}}$ gets longer, $C(t_{\text{MC}})$ approaches to its equilibrium value.

To accelerate equilibration, we use a $\tau$-doubling trick, which is very similar to the $\beta$-doubling in standard SSE algorithm Sandvik (2002). For a given evolution time $\tau_{\text{max}}$, the $\tau$-doubling process starts from sampling a much shorter time $\tau=\tau_{\text{min}}$ with $t_{\text{MC}}$ (typically a few hundreds) and doubles the evolution time to $\tau=2\tau_{\text{min}}$ by sampling the doubled operator sequence $S_{2n}$ Sandvik (2002). Repeating the doubling until $\tau=\tau_{\text{max}}$, the initial configuration is almost equilirated since it has already rather long history at shorter times Sandvik (2002). For the data presented in the main text, to be safe, at the last doubling, we use $t_{\text{MC}}=10^{4}\sim 10^{5}$ to perform the equilibrations and measurements for $50\sim 200$ bins ( one bin consists of $10^{4}\sim 10^{5}$ MC steps, and the bin number depends on the system sizes and the quantity studied). The statistical errors are estimated as one standard deviation of the average of the bin averages, which are reliable only when $t_{\text{MC}}$ is much longer than the autocorrelation time so that the bin averages can be regarded as statistically independent.

Next we discuss autocorrelation function, which, for an observable $O$, is defined as Landau and Binder (2009)

[TABLE]

where the averages are over the MC time index $i$ and $\Delta t$ is the separation of the two configurations. $A(0)$ is normalized to $1$ and when $\Delta t$ is large, $A(\Delta t)$ approaches [math] exponentially

[TABLE]

where $\Theta$ is called the autocorrelation time. Figure B.2(a) shows the autocorrelation functions of $C(\tau,L)$ at $\tau=16$ for the 1d TFIM with system sizes $L=16,32,64$ and $128$ under the initial condition of $M_{0}=0$ and a very short correlation length. The autocorrelation times are extracted by fitting the long time behaviors of $A(\Delta t)$ using the form in Eq. (28), as shown in the inset of Fig. B.2(a). In fact, we expect similar behaviors of $A(\Delta t)$ and $\Theta$ for the 2d TFIM, as shown in Fig. B.2(b), in which $A(\Delta t)$ is computed for $C(\tau,L)$ with $\tau=4$ and $L=12,16,20,24$. Since the initial state has a vanishing correlation length, after a short length of imaginary time evolution (not sufficient to reach the ground state), the correlation length is still finite in the short time regime. When the system size grows larger than the correlation length, the autocorrelation time becomes insensitive to the system size, indicating that the critical slowing down problem is overcome. The autocorrelation times are only up to a few MC steps so that in the simulations, the number of steps we used $t_{\text{MC}}=10^{4}$ to $10^{5}$ is sufficient to generate statistically independent configurations.