Testing the event-chain algorithm in asymptotically free models

We apply the event-chain algorithm proposed by Bernard, Krauth and Wilson in 2009 to toy models of lattice QCD. We give a formal prove of stability of the algorithm. We study its performance at the example of the massive Gaussian model on the square and the simple cubic lattice, the $O(3)$-invariant non-linear $\sigma$-model and the $SU(3) \times SU(3)$ principle chiral model on the square lattice. In all these cases we find that critical slowing down is essentially eliminated.


I. INTRODUCTION
Monte Carlo simulation of statistical field theory systems is a standard tool in their study. While the development of algorithms has a long tradition, it is still an active field and source of substantial progress. Here we test a new class of recently proposed algorithms, the so-called event-chain Monte Carlo [1][2][3][4], in systems with asymptotic freedom. In particular, we study the two-dimensional OðNÞ-invariant nonlinear σ-model and the two-dimensional SUðNÞ × SUðNÞ principal chiral model. In our simulations, we take N ¼ 3 in both cases. As a preliminary step, we study free field theory in two and three dimensions.
The cost of a simulation is determined essentially by two factors: the numerical cost of one update of the whole system and the autocorrelation time τ in units of these updates. For definitions of the autocorrelation time see Appendix B. It typically scales with a power of the correlation length ξ, The dynamical critical exponent z therefore characterizes the algorithm's performance close to criticality. Algorithms, such as the local Metropolis algorithm, which essentially perform a random walk in configuration space have z ≈ 2, whereas for some models also update strategies whose autocorrelation times increase much slower with the correlation length have been devised. Examples are the cluster [5,6] and multigrid [7] algorithms. They are characterized by a coherent update of a large fraction of the field variables in one step. Such advanced algorithms, however, rely on special features of the theory to be simulated and are difficult to generalize: Algorithm performance depends strongly on the model under investigation. For a new proposal, it is therefore pivotal to test it in many cases to understand better its dynamics. An interesting new entry in the toolbox is the so-called event-chain algorithm pioneered in Refs. [1][2][3]. It is basically a walker on the lattice which updates the local field variable in Monte Carlo time until a so-called event occurs, with probabilities given by the theory. Then it moves to a neighboring site, which is determined by the event. The details are given below. It is remarkable that for some theories this leads to a small value of z and not to the random walk one might naively expect.

II. EVENT DRIVEN ALGORITHM
In this section, the algorithm is discussed for a model with a single real variable at each site of the lattice. In order to apply the algorithm to more complicated models, the simple model is embedded, as we explain below. Note that the idea of embedding is also used in the case of cluster [6,8] and multigrid [9] algorithms. To define the algorithm, one first discusses infinitesimal updates with a step size ϵ > 0. These updates are integrated analytically, leading to the version of the algorithm that is implemented at the end.
We consider a general scalar field theory with the action where s x ðϕ x Þ and s hx;yi ðϕ x ; ϕ y Þ are functions of the real variables ϕ x . Here we consider analytic functions. However the example of hard spheres considered in Ref. [1] shows that this requirement can be relaxed. The collection of all sites of the lattice is denoted by Λ and hx; yi is a pair of interacting sites. In our numerical study, we consider nearest neighbor interactions only. We require that holds for the interaction terms. The configuration space ϕ is now extended by two "lifting" variables: the position of the walkerx ∈ Λ at which field updates are performed and the direction σ ∈ f−1; 1g of these changes. Both are drawn from a uniform distribution, in particular, uncorrelated to the fields ϕ. Let us denote an enlarged configuration by The action pertaining to a single site, keeping the fields at all other sites fixed, can be written as where i ¼ 0 indicates the single site term and i ¼ 1; 2; …; n labels the interaction partners y of x.
The enlarged configuration X is always changed. Either ϕx is replaced by ϕx þ σϵ, σ is replaced by −σ, orx is replaced by one of its interacting partners y i . The decisions are taken for each term of the single site action (5) separately. With the standard Metropolis probability, applied to s i;x the proposal ϕx þ σϵ is accepted. Expanding in ϵ we get The update ϕx þ σϵ is rejected, if it is rejected by at least one s i;x . If the proposal is not accepted the following is done: In the case of i ¼ 0, σ is replaced by −σ, while for i > 0 the walker is set to the site y i . At finite ϵ there is a probability Oðϵ 2 Þ that the proposal is rejected by more than one s i;x . Hence for a sequence of n updates with t ¼ nϵ finite, the probability for conflicting decisions is OðϵÞ and hence the algorithm becomes well defined in the limit ϵ → 0. A proof of the correctness of the algorithm is given in Appendix A. The algorithm that is implemented in the program is obtained by integrating infinitesimal steps over a finite fictitious time t. The probability that there is no rejection by s i;x in n ¼ t=ϵ subsequent steps is Taking the logarithm and the limit ϵ → 0 one arrives at Performing the integral we get where k labels the intervals with σs 0 i;x being positive throughout within ½min½ϕx; ϕx þ σt; max½ϕx; ϕx þ σt. ϕ 0;k and ϕ 1;k are the lower/upper and upper/lower bounds of these intervals, depending on the value of σ. The field ϕx is updated until the update is rejected by one of the s i;x .
The times t ðeventÞ i when these events occur are determined in the following way: One draws a uniform random number r i ∈ ð0; 1 for each i to fix P i ðtÞ. We arrive at Eq. (11) of [2], which has to be solved for t . The field is updated ϕx → ϕx þ σt ðeventÞ . For i min ¼ 0 we replace σ by −σ andx remains unchanged. Or else, for i min > 0 the walker assumes the new position y i min and σ keeps its sign.
The algorithm evolves in a fictitious time t MC . For each event the fictitious time increases as t MC → t MC þ t ðeventÞ . A sequence of updates is started at t MC ¼ 0 by drawingx and σ from a uniform distribution. The sequence of updates is stopped at some fixed t f . To this end, events are generated until t MC overshoots t f for the first time. This last event is not taken into account in the update of the field and t MC is not increased by this last t ðeventÞ . Instead ϕx → ϕx þ σðt f − t MC Þ is performed. Furthermore, measurements should be performed in equal intervals of the fictitious time. Mostly, we performed a measurement after a complete update sequence of the length t f or after a fixed number of these sequences. Note that measuring at events leads to a bias that, at least for small system sizes, can be easily seen in the averages of estimators.
In order to arrive at a formal proof of ergodicity one has to introduce some randomness in the evolution time t f . An alternative would be to add Metropolis updates that ensure ergodicity. In our explorative study here, similar to Refs. [1][2][3], we simply ignore this question.

III. THE MODELS
We consider square or simple cubic lattices. The sites of the lattice are denoted by x ¼ ðx 0 ; x 1 ; …; x d−1 Þ with x i ∈ f0; 1; 2; …L i − 1g, where d is the dimension of the system, setting the lattice spacing to unity, a ¼ 1. The directions of the lattice are denoted by μ ∈ f0;1;…;d−1g. There are 2d nearest neighbors y i ¼ x AEμ, whereμ is a unit vector in μ-direction. In our simulations we take

A. Free field theory
First we study the two-and three-dimensional scalar field theory. It is defined by the action In the following we only discuss the free case λ ¼ 0.
Here the functions (9) are easy to evaluate. For the single site term of the action we get For i > 0 we get There is always a solution for a positive t.

B. Nonlinear σ-model
The action of the nonlinear OðNÞ-invariant σ-model is where ⃗ s x is a unit vector with N real components. We study the model on the square lattice for N ¼ 2 and 3. For N ¼ 2 we get the so-called XY-model. It undergoes a Kosterlitz-Thouless phase transition at β KT ¼ 1.1199ð1Þ [10]. For temperatures above the transition, there is a finite correlation length. At lower temperatures there is no ordering and the two-point correlation function is decaying with a power law.
For references see e.g., [11]. For the event-chain algorithm it is useful to write the XY-model in terms of angles α x , The pair terms of the action are In the update, α˜x is incremented. The functions (9) are given in Ref. [2]. To simplify the notation, let us assume σ ¼ 1 in the following. Let us define where the integer m is chosen such that −π < δ i ≤ π. 2n where the integer n ≥ 0 is chosen such that −π < δ i þ t − 2nπ ≤ π. Note that 2n is the contribution from n complete cycles, where cosð0Þ − cosðπÞ ¼ 2. For an illustration see Fig. 4 of Ref. [2].
For N ¼ 3 the model has a finite correlation length at any finite value of β. The model is asymptotically free. The divergence of the correlation length as β → ∞ is governed by the so-called β-function. For a discussion of the physics of this model see e.g., [12].
In the literature there are a number of Monte Carlo studies of this lattice model. It can be efficiently simulated by using the cluster algorithm [13] and the multigrid algorithm [9]. Also the microcanonical overrelaxation algorithm [14] has been applied successfully.
For simulations with the worm algorithm see Ref. [15]. The simulation for N > 2 with the event-chain Monte Carlo algorithm is performed, by an embedding of a generalized XY model. To this end, for one sequence of updates, one picks out a pair ðl; kÞ of components from the N components of the field. In the update, only rotations in this plane are performed. The algorithm is run as for the XY model, with the exception that in Eq. (19) the prefactor β is replaced by C. SUðNÞ × SUðNÞ principal chiral model The principal chiral model on the square lattice is defined by the action where U x ∈ SU (N). The action restricted to a single site is The model is asymptotically free. For a discussion of the physics of this model see e.g., [16]. No efficient implementation of the cluster algorithm has been devised for this model so far. However the multigrid algorithm has been implemented successfully [16,17]. To this end, in Ref. [16], for one sequence of updates, a one-parameter subgroup of SUðNÞ is considered. Here we apply the same idea. Following [16], a general one-parameter update can be composed of a random SUðNÞ rotation and a rotation by a variable angle along one, fixed element λ of the algebra. With leading to an embedded action for site x and component i, where

IV. NUMERICAL RESULTS
Below we discuss our numerical results for the scalar free field theory, the Oð3Þ-invariant nonlinear σ-model and the SUð3Þ × SUð3Þ-invariant principal chiral model. We abstain from a discussion of our simulations of the twodimensional XY-model, since our results are fully consistent with those of Ref. [18].

A. Scalar free field theory
We simulated the Gaussian model on the square and simple cubic lattice. In our simulations, we have taken L ∝ m −1 , evolving the fictitious Monte Carlo time for the period t f using the event-chain algorithm. Once this is reached, a new sitex and direction σ are selected. Note that t f is the only free parameter of the algorithm.
Let us define the quantities that we measured in our simulations. The first two quantities are related to the terms in the action, Furthermore we determined the magnetization and a proxy of the second moment correlation length, with the Fourier transformφ i ¼ P x expð−i2πx i =LÞϕ x . In the Gaussian model, these quantities can be easily computed exactly by Fourier transformation. We compared the outcome of our simulations with these exact results. The largest deviation had been 3.4 standard deviations for a single observable.
We performed a measurement after n m sequences of the length t f . For a given L we kept t meas ¼ n m t f constant. It is chosen such that the integrated autocorrelation times in units of t meas are of Oð10Þ. Studying local algorithms, such as the local Metropolis algorithm, one typically quotes the autocorrelation times in units of sweeps over the whole lattice. In our case, the artificial time that is needed on average for L d events could be used as an analog. As discussed below, we find that the number of events per unit of the artificial time is slightly larger than one for both the square as well as the simple cubic lattice. Similar observations are made in the case of the Oð3Þ-invariant nonlinear σ-model and the SUð3Þ × SUð3Þ-invariant principal chiral model.

Two dimensions
In the two-dimensional case we simulated the linear lattice sizes L ¼ 32, 64, 128, and 256. Correspondingly the masses are taken as m ¼ 0.3, 0.15, 0.075, and 0.0375, respectively. Each series of simulations consists of 10 6 measurements, separated by t meas ¼ 100, 400, 2000 and 8000 for L ¼ 32, 64, 128, and 256, respectively. We focus on the dependence of the autocorrelations on the length t f of the update sequence.
In the left part of Fig. 1 we give the integrated autocorrelation time of O 3 . We determined the integrated autocorrelation time in units of measurements. In the plot, we give τ int in units of L 2 . To this end, we multiply the numbers that we compute by t meas =L 2 . Plotting the autocorrelation time as a function of t f =L 2 , we see a reasonable collapse of the data for the four different lattice sizes. We conclude that t f should be chosen such that t f ∝ m −2 ∝ L 2 . The integrated autocorrelation time of O 3 decreases with increasing t f and seems to approach a plateau value for large t f . We cannot exclude that τ int increases for very large t f again, as it is the case for the principal chiral model, discussed below.
A perfect collapse of the data for the different lattices sizes would also imply that critical slowing down is eliminated completely. However we see a slight increase of τ int with increasing L. This is studied below for a fixed value of t f =L 2 using higher statistics. Now let us briefly comment on the behavior of τ int for the observables that we do not show in the plot. The behavior of τ int for O 4 is similar to that of O 3 . In the case of O 1 and O 2 the behavior is different. For O 1 we see an increase of τ int with increasing t f . This increase is not dramatic and τ int seems to approach a plateau. For O 2 we see a decrease of τ int with increasing t f . It is however less pronounced as for O 3 and O 4 .
Next we performed simulations for t f =L 2 ¼ 125=4096≈ 0.030517578… fixed. Measurements were performed after t meas ¼ 2t f . Here we performed 5 × 10 6 measurements for each lattice size. In addition to L ¼ 32, 64, 128, and 256 we have simulated L ¼ 512. The simulation for L ¼ 512 took 8 hours on one core of a Intel(R) Xeon(R) CPU E5-2660 v3 2.60 GHz. The results for the integrated autocorrelation times are summarized in Table I.
Looking at the numbers we see that for each observable, the integrated autocorrelation time is increasing with increasing lattice size. Fitting the numbers for the observable O 1 with the ansatz τ int;1 ¼ cL z , taking into account the linear lattice sizes L ≥ 32, we get z ¼ 0.112ð4Þ and χ 2 =d:o:f: ¼ 1. 15. Fits for the other three observables result in even smaller values for z. Fits with the ansatz τ int ¼ a þ c lnðLÞ seem to be even more adequate. For example, fitting all data for τ int;1 results in c ¼ 1.05ð2Þ and χ 2 =d:o:f: ¼ 0.68. We conclude that slowing down is essentially eliminated.
The CPU time that is used is proportional to the number of events. There is a small dependence on the mass and the lattice size of the average number of events within one update sequence divided by t f . Roughly this ratio equals 1.128 for all our masses and lattice sizes.
Note that our results are consistent with those of Ref. [18], where the massless case was studied.

Three dimensions
Next we simulated the three-dimensional Gaussian model. We simulated the linear lattice sizes L ¼ 16, 32, 64, and 128 for various values of t f . Similar to the twodimensional case, we use the masses m ¼ 0.6, 0.3, 0.15, and 0.075. In Fig. 1 on the right, we give the integrated autocorrelation time of the observable O 3 . Also in three dimensions it turns out that the curves for the different L fall approximately on top of each other for plotting τ int as a function of t f =L 2 . Since we have taken L ∝ m −1 , it is likely that actually m 2 t f should be kept constant.
As for two dimensions, we have measured the number of events per t MC . We get here ≈1.38. For t f =L 2 ¼ 1=4 ¼ 0.25 we performed simulations with 5 × 10 6 measurements. The results of these simulations are summarized in Table II. In contrast to two dimensions, the autocorrelation times are slightly decreasing with increasing lattice size (decreasing mass). Hence slowing down is eliminated.
Our results for the observables and the corresponding integrated autocorrelation times are summarized in Table III. The integrated autocorrelation times are given in units of three XY embeddings. In the last column we give the average number of events divided by the number of sites for one cycle of three XY embeddings. The number slowly increases with increasing lattice size. This has to be taken into account when judging the performance of the algorithm. The fact that the number is close to 1 in all cases means that in one cycle, each site is touched roughly once.
In the case of the energy density E, the integrated autocorrelation time even decreases with increasing correlation length ξ. On the other hand, for the susceptibility χ, the integrated autocorrelation time is increasing with increasing ξ. Fitting with the ansatz τ int;χ N ev =L 2 ¼ cξ z gives z ¼ 0.386ð5Þ and χ 2 =d:o:f: ¼ 7.74, taking all data into account and z ¼ 0.367ð7Þ and χ 2 =d:o:f: ¼ 4.93, taking all data with β ≥ 1.5 into account. On the other hand, fitting with the ansatz τ int;χ N ev =L 2 ¼ a þ c lnðξÞ gives c ¼ 0.74ð1Þ and χ 2 =d:o:f: ¼ 7.82, taking all data into account and c ¼ 0.80ð2Þ and χ 2 =d:o:f: ¼ 1.33, taking all data with β ≥ 1.5 into account. The logarithmic increase of the integrated autocorrelation time is favored by the fits. In any case, slowing down is essentially eliminated, since also the fits with a power law ansatz result in small values for z.
Note that in Ref. [3] slowing down with z ≈ 1 has been observed for the Oð3Þ-invariant model on the threedimensional simple cubic lattice. Furthermore one should notice that the single cluster algorithm [13] is more efficient than the event-chain Monte Carlo algorithm and on top of that provides variance reduced estimators for various observables.
C. The SUð3Þ × SUð3Þ-invariant principal chiral model For an update sequence of the length t f a fixed random rotation matrix R is used. Then for the next update sequence a new R is generated. We measure the following observables, constructed from fundamental and adjoint correlation functions [16] The energies E F and E A and the susceptibilities χ F and χ A as well as the second moment correlation length In Table IV we give the basic parameters of our runs and our result for the correlation length ξ F . We have chosen the linear lattice size L such that L=ξ F ≈ 11 for all values of β. We performed a large number of simulations, varying the length t f of the update sequence. We find that the results from various runs are compatible and agree also with results from the literature [16,17], confirming the correctness of our implementation of the algorithm.
In this particular study, we also measure in time intervals which are smaller than t f in order to reveal the scaling behavior with long update sequences. Our estimates of the integrated autocorrelation times are shown in Fig. 2. We give both the autocorrelation time and the length t f of the update sequence normalized by the average number of events divided by the number of sites L 2 . We observe a rather weak dependence on the length t f of the update sequence, with a shallow minimum. In particular for sequences with a length t f such that 0.1 up to 1 events per site occur, we see very little variation.
Next we study in more detail how the autocorrelation times scale with the correlation length ξ F . To this end, for each value of β, we determine t f at which the integrated autocorrelation time is minimal. In Fig. 3, the corresponding minimal integrated autocorrelation times for four different observables are plotted against the correlation length.
The behavior of τ int;χ A and τ int;χ F can be described reasonably well by a power law behavior with z ≈ 0.3 and 0.5, respectively. This behavior is indicated by the two solid lines in Fig. 3. In the case of the energies E A and E F , the integrated autocorrelation is even decreasing, going from β ¼ 1.65 to 1.985. We conclude that the dynamical critical exponent z is well below 1 and hence critical slowing down is essentially eliminated.

V. CONCLUSIONS AND OUTLOOK
We applied the event-chain Monte Carlo algorithm [1][2][3][4] to asymptotically free models in two dimensions. In the literature, these models are considered as toy models of lattice QCD that allow for testing of algorithmic and theoretical ideas in a simplified setting. In a preliminary study, we applied the event-chain Monte Carlo to free field theory on a square and a simple cubic lattice. We find that critical slowing down is eliminated. This is quite astonishing for an algorithm, where local changes of the field variables are performed. Next we studied the twodimensional XY-model. We do not discuss our results in detail. They corroborate the findings of [18]: At low temperatures, in the spin wave phase, slowing down is eliminated, which is consistent with our free field theory result. On the other hand, in the high temperature phase, the physics is governed by vortices. Here we find slowing down with z ≈ 2, even though the amplitudes of the integrated autocorrelation times are considerably reduced compared with the local Metropolis algorithm. One should note however that for this model the over-relaxation algorithm [14] gives z ≈ 1 in both phases and the cluster algorithm eliminates slowing down completely [13,20]. In order to apply the event-chain Monte algorithm to the Oð3Þ-invariant nonlinear σ-model and the SUð3Þ × SUð3Þ chiral model, we use an embedding of the XY-model. In the case of the Oð3Þ-invariant model this simply means that one picks out two of the three components of the field. During one sequence of updates only these two components are updated. For the next sequence a new pair of components is taken. In the case of the SUð3Þ × SUð3Þ principal chiral model the embedding is a bit more complicated. We follow Ref. [9] whose authors used such an embedding in the context of the multigrid algorithm. For details see Sec. III C.
We find for both the

APPENDIX A: PROOF OF STABILITY
We have to show that the infinitesimal update by ϵ preserves the target distribution where Z is the partition function of the lattice model and the factors 2 and V give the number of possible values of σ andx. Let us write the update probabilities defined in Sec. II as TðY ← XÞ. Then stability means where P X is a shorthand for the integral over the field ϕ and the sums overx and σ. As a first step, we list the configurations X that can end up in a given Y ¼ ðϕ; σ;xÞ, after an update step.
(i) X 0 given by ϕ the same as for Y, σ is replaced by −σ and the position of the walkerx remains the same. (ii) X i given by ϕ is the same as for Y, σ is the same as for Y andx ¼ y i , meaning that the walker is hopping to its interaction partner y i . (iii) X nþ1 given by ϕx − σϵ, ϕ x with x ≠x are the same as for Y, σ andx keep their value. Hence Eq. (A2) can be written as TðY ← X i ÞπðX i Þ: ðA3Þ Note that the probability density of X i for i ≤ n is the same as that of Y, since the field ϕ is the same. The probability density of X nþ1 is Now let us work out the transition probabilities, Note that for i ¼ 0 the sign of σ changes, while for 1 ≤ i ≤ n the hopping of the walker along with Eq. (3) produces a minus sign. Finally Now we can put things together, Contributions OðϵÞ exactly cancel each other.

APPENDIX B: AUTOCORRELATION TIMES
The autocorrelation function of an observable A is given by where the index of A gives the position in the Markov chain. The modulus of the autocorrelation function is bounded from above by an exponentially decaying function. Following Ref. [21] we define and the exponential autocorrelation time as which characterizes the Markov chain. Since τ exp is hard to compute, we determine integrated autocorrelation times. The integrated autocorrelation time of the observable A is given by The statistical error of the estimate of hAi is where n is the number of measurements and is the variance. Computing τ A , we use the estimate of ρ A ðtÞ obtained from our simulation. Therefore the summation in Eq. (B4) has to be truncated at some finite t max . Since ρ A ðtÞ is falling off exponentially at large distances, the relative statistical error becomes large at large distances. Therefore it is mandatory to truncate the summation at some point that is typically much smaller than the total length of the simulation. In the literature one can find various recommendations on how this upper bound should be chosen. Sokal [21] suggests to take t max ≈ Mτ int , where the value of M should be at least 6.