# «A quantum Monte Carlo algorithm realizing an intrinsic relaxation Tota Nakamura and Yoshiyuki Ito Department of Applied Physics, Tohoku University, ...»

typeset using JPSJ.sty ver.1.0b

arXiv:cond-mat/0301560v2 [cond-mat.stat-mech] 14 May 2003

A quantum Monte Carlo algorithm realizing an intrinsic relaxation

Tota Nakamura and Yoshiyuki Ito

Department of Applied Physics, Tohoku University, Aoba, Sendai, Miyagi 980-8579

(Received )

We propose a new quantum Monte Carlo algorithm which realizes a relaxation intrinsic to

the original quantum system. The Monte Carlo dynamics satisﬁes the dynamic scaling relation

τ ∼ ξ z and is independent of the Trotter number. Finiteness of the Trotter number just appears as the ﬁnite-size eﬀect. An inﬁnite Trotter number version of the algorithm is also formulated, which enables us to observe a true relaxation of the original system. The strategy of the algorithm is a compromise between the conventional worldline local ﬂip and the modern cluster loop ﬂip. It is a local ﬂip in the real-space direction and is a cluster ﬂip in the Trotter direction. The new algorithm is tested by the transverse-ﬁeld Ising model in two dimensions.

An accurate phase diagram is obtained.

KEYWORDS: quantum Monte Carlo method, nonequilibrium relaxation method, quantum dynamics, transverse- ﬁeld Ising model The Monte Carlo method is a powerful tool to investigate the condensed matter physics. The method directly beneﬁts from a progress of microprocessors. It is also suitable for the parallel computing. Developments of various Monte Carlo algorithms consist of an important part of the statistical physics in these decades

1 are simultaneously ﬂipped. This loop corresponds to the correlated spins in this (d+1)-dimensional lattice. Therefore, the algorithm drastically reduces the correlation time.

The loop algorithms, or generally cluster algorithms, are a trend of recent developments in Monte Carlo methods. It accelerates the Monte Carlo simulations suﬀering from the slow dynamics.

Another trend is the nonequilibrium relaxation (NER) method.7–9) The NER method positively makes use of the critical slowing down to detect the phase transition. One observes a relaxation function of a physical quantity at each temperature. If it exhibits an algebraically slow relaxation, the temperature is judged to be a critical point. The method works very well in the slow-dynamic systems.10–12) A Monte Carlo investigation on the relaxation process is a good approach to the critical phenomena.

An updating algorithm in the NER simulation must realize an intrinsic relaxation of the system.

For example, an algebraic relaxation should be observed at the critical point. A relaxation in the paramagnetic phase should be exponential with a ﬁnite correlation time corresponding to a ﬁnite correlation length. The requirement is guaranteed by updating degrees of freedom with the same length at each time step. One realization is a single spin ﬂip algorithm in the classical spin systems.

In regard to the QMC method, however, an algorithm which realizes a relaxation intrinsic to the original quantum system has not been invented yet. A relaxation function by the loop algorithm has nothing to do with the physics. The NER analysis is possible by using the worldline local ﬂip.13) A ratio between the inverse temperature J/T and the Trotter number m is ﬁxed ﬁnite. The simulation exhibits a relaxation of the (d + 1)-dimensional classical system. It suﬀers from the freezing as the Trotter number increases.

Aim of this Letter is to propose a new QMC algorithm realizing the relaxation intrinsic to the original quantum system. A shortcoming of the worldline local ﬂip algorithm is a mixture of a physical relaxation and an unphysical relaxation causing the freezing. The new algorithm extracts only a physical relaxation part. The unphysical one is eliminated by using an idea of the loop algorithm. The inﬁnite Trotter number version is also formulated.

A basic idea of the present algorithm is to make an update local in the real-space direction and global in the Trotter direction. The former one is to ensure the Monte Carlo dynamics reﬂecting the physics of the original quantum system. We typically choose one real-space interaction bond or one real-space spin as a local unit of the updating. One Monte Carlo step corresponds to a time unit in which a change of this real-space length occurs. The latter idea is realized by a cluster ﬂip extending only in the Trotter direction. Each updating unit is connected along the Trotter direction in a same manner as the Swendsen-Wang algorithm14) in one dimension. A cluster of the connected updating units is ﬂipped simultaneously. A length of the cluster corresponds to the correlation in the Trotter direction. Therefore, the freezing due to the Suzuki-Trotter decomposition is solely eliminated. The present update algorithm can be considered as a single quantum-spin ﬂip.

Fig. 1. (a)The Suzuki-Trotter decomposition of the Heisenberg chain. The worldlines are depicted by bold lines.

Pseudo-spins are depicted by arrows. Those updated in this trial are depicted by bold arrows. (b) Deﬁnitions of a pseudo-spin in a blank plaquette due to the worldline conﬁguration. No pseudo-spin is assigned to the worldline vacancy and the double occupancy. (c) The Suzuki-Trotter decomposition of the transverse-ﬁeld Ising chain. Spins updated in this trial is depicted by solid circles. Rectangles indicate clusters.

We consider an S = 1/2 Heisenberg chain for a simple explanation of the algorithm.

The system is decomposed into a checkerboard plane as shown in Fig. 1(a). An updating unit is four spins consisting of a blank plaquette in this ﬁgure. We deﬁne a pseudo-spin in a blank plaquette as depicted in Fig. 1(b). If a worldline runs vertically in the left (right) side of the plaquette, a pseudo-spin pointing to the right (left) is assigned. If there is no vertical worldline or there are two worldlines in both sides, no pseudo-spin is assigned. A pseudo-spin represents a possible worldline move. It changes the direction, if the spin ﬂip is accepted. The worldline local ﬂip is regarded as a single spin ﬂip of this pseudo-spin.

An actual procedure of the updating is as follows. We choose one real-space bond to try an update (the double-line bond connecting two solid spins in Fig. 1(a)). Pseudo-spins in regard with these two real-space spins are updated. (Bold arrows in the ﬁgure.) A Boltzmann weight W

One consists of weights of plaquettes on the updating bond ( ) as denoted by Wp. These plaquettes can be considered as eﬀective bonds connecting pseudo-spins on the updating bond. They are used to deﬁne a cluster of pseudo-spins. The other part consists of weights of plaquettes on the neighboring bonds ( ) as denoted by A. The weights determine a probability to accept the cluster ﬂip.

If two pseudo-spins neighboring along the Trotter direction point to the same direction, they are connected to form a cluster with a probability 1 − p of

** p≡ / = tanh(J/mT ) ∼ J/mT. (3)**

Pseudo-spins in a cluster are ﬂipped simultaneously with a probability calculated from the weights of plaquettes on the neighboring bonds A. A product of the weights of plaquettes adjacent to each pseudo-spin in the cluster is calculated before (Ainitial ) and after (Aﬁnal ) the ﬂip. We accept the ﬂip with a probability Aﬁnal. (4) Ainitial + Aﬁnal We try this ﬂip for each cluster independently. A new worldline conﬁguration is obtained by new pseudo-spin conﬁgurations.

The ergodicity of this update algorithm is same as the worldline local ﬂip. The worldline global ﬂip is necessary to change the magnetization. The present algorithm is an adoption of the SwendsenWang algorithm under the external ﬁeld. Here, a molecular ﬁeld from the neighboring real-space spins is regarded as the external ﬁeld. Therefore, the detailed balance is also satisﬁed as is guaranteed in the Swendsen-Wang algorithm. An important notice is that the acceptance of a ﬂip does not depend on the Trotter number nor the temperature. An average size of the cluster is about 1/p ∼ mT /J. Then, the contribution from the neighboring interaction bonds A is an order of 1 J exp ∼ O(1). The freezing due to the decomposition is solely eliminated by this relation.

× mT p A cluster analysis is easy because it is performed only in one dimension. It enables us to write a fast program without diﬃculties for various systems. An actual computational time of one Monte Carlo step is almost same as that of the conventional local ﬂip one.

It is possible to take the inﬁnite Trotter number limit beforehand. In this scheme we introduce a ‘breakup’, which is a domain wall of the pseudo-spins, i.e., where a worldline hops.3–5) A location of a breakup is stored in memory to deﬁne a state. A probability of a breakup to exist is p ∼ J/mT. An average number of the breakups is J/T in the inﬁnite m limit. Therefore, among the pseudo-spins on a selected updating bond we put breakups by Poisson random numbers with the expectation number J/T. The pseudo-spins between the neighboring breakups are ﬂipped by using the probability of 4 Eq. (4). Here, we do not discriminate between the newly-assigned breakups and the already-existing breakups.

An updating procedure for the transverse-ﬁeld Ising model is almost same as for the Heisenberg model. Only a diﬀerence is that there is no local spin conservation in this model. An updating unit, which is a plaquette in the Heisenberg model, shrinks to a single spin. Therefore, the single-spin ﬂip is possible in the (d + 1)-dimensional lattice. The Hamiltonian is

The Suzuki-Trotter decomposition is shown in Fig. 1(c). A Boltzmann weight is assigned to each interaction bond. First, we choose one real-space spin to try an update. Spins along a line in the Trotter direction will be updated. (Solid circles in Fig. 1 (c).) A cluster is deﬁned by connecting two spins neighboring in the Trotter direction with a probability 1 − p of

** p = tanh(Γ/mT ) ∼ Γ/mT, (6)**

if they take same spin value. The clusters are shown by rectangles in Fig. 1(c). They are ﬂipped with a probability of Eq. (4). The weights Ainitial and Aﬁnal are calculated by molecular ﬁelds from the spins neighboring in the real-space direction to the cluster. The inﬁnite Trotter number version in this model is also trivial. We put breakups with the expectation number Γ/T along a line in the Trotter direction of the updated spin. A cluster between two breakups is ﬂipped by Eq. (4).

As a demonstration of the new algorithm we consider the transverse-ﬁeld Ising model in two dimensions. The classical phase transition occurs at a temperature which varies with Γ/J. The quantum phase transition occurs at T = 0 for a ﬁnite Γ/J value.

Figure 2 shows NER plots of the magnetization M at T /J = 0.1 and Γ/J = 3.05. This point is in the vicinity of the transition point. The relaxation functions depend on the Trotter number in the conventional algorithm as shown in Fig. 2(a). As the Trotter number increases, it becomes later for the algebraic relaxation to begin. On the other hand, the relaxation functions by the new algorithm essentially do not depend on the Trotter numbers as shown in Fig. 2(b). The relaxation function of the inﬁnite Trotter number is depicted by a line. The ﬁgure clearly shows that the Trotter number dependence only appears as the ﬁnite-size eﬀect. Before it appears, the relaxation coincides with that of the inﬁnite Trotter number. The relaxation of the original quantum system is considered to be realized. The equilibrium data of small Trotter numbers are consistent between two algorithms. The logarithmic slope of the relaxation −β/zν also coincides with each other.

Therefore, we can conclude that two algorithms are observing the same critical phenomenon. A computational time for one Monte Carlo step is almost same between two algorithms. The total necessary computation time to achieve the same resolution is reduced to 1/10.

5 0 10

Fig. 2. NER plots of the magnetization of 2D transverse-ﬁeld Ising model near the critical point. Γ/J = 3.05, T /J = 0.1. The real-space system size is 400 × 400. (a) By the conventional local ﬂip algorithm. (b) By the new ﬂip algorithm proposed in this Letter. Data of the inﬁnite Trotter number version (m = ∞) is depicted by a line.

Figure 3 shows a phase diagram obtained by the NER analysis of the magnetization using the inﬁnite Trotter number version of the present algorithm. The smallest value of Γ at which the relaxation exhibits the exponential decay is the lower bound of the paramagnetic phase. The largest value of Γ at which the relaxation converges to a ﬁnite value is the upper bound of the ferromagnetic phase. The phase boundary line is in between. The lowest temperature we have carried out the simulation is T /J = 0.01. Within the system size (149 × 150) and the time range (1000 steps) of the simulation at this temperature there is no temperature eﬀect and the data can be considered as those of the ground state. Our estimate of the quantum transition point is Γc = 3.044(2). In the classical limit (Γ/J = 0) the phase transition of the 2D classical Ising 6 model occurs at T /J = 2.269. We have also veriﬁed this temperature and the critical exponents by using the same program. The phase diagram is obtained very accurately compared with previous investigations.15) The shape of the phase diagram qualitatively agrees with the experimental result of LiHoF4,16) which is considered to realize the transverse-ﬁeld Ising model in three dimensions.

Fig. 3. A phase diagram of 2D transverse-ﬁeld Ising model obtained by the NER using the inﬁnite Trotter number algorithm. The real-space system size is 199 × 200 (T /J ≥ 0.1) and 149 × 150 (T /J = 0.01). The time range is 1000 steps.

Figure 4 shows an NER plot on the phase boundary line when changing the temperature. The logarithmic slopes −β/zν converge to the same value in the long time limit at high temperatures.