Antiferromagnetic Ising model in an imaginary magnetic field

# Antiferromagnetic Ising model in an imaginary magnetic field

Vicente Azcoiti Departamento de Física Teórica, Facultad de Ciencias, Universidad de Zaragoza
C/Pedro Cerbuna 12, E-50009, Zaragoza (Spain)
Giuseppe Di Carlo INFN, Laboratori Nazionali del Gran Sasso,
I-67100 Assergi, L’Aquila (Italy)
Eduardo Follana Departamento de Física Teórica, Facultad de Ciencias, Universidad de Zaragoza
C/Pedro Cerbuna 12, E-50009, Zaragoza (Spain)
Eduardo Royo-Amondarain Departamento de Física Teórica, Facultad de Ciencias, Universidad de Zaragoza
C/Pedro Cerbuna 12, E-50009, Zaragoza (Spain)
###### Abstract

We study the two-dimensional antiferromagnetic Ising model with a purely imaginary magnetic field, which can be thought of as a toy model for the usual physics. Our motivation is to have a benchmark calculation in a system which suffers from a strong sign problem, so that our results can be used to test Monte Carlo methods developed to tackle such problems.

We analyze here this model by means of analytical techniques, computing exactly the first eight cumulants of the expansion of the effective Hamiltonian in powers of the inverse temperature, and calculating physical observables for a large number of degrees of freedom with the help of standard multi-precision algorithms. We report accurate results for the free energy density, internal energy, standard and staggered magnetization, and the position and nature of the critical line, which confirm the mean-field qualitative picture, and which should be quantitatively reliable, at least in the high-temperature regime, including the entire critical line.

## I Introduction

One of the major challenges for high-energy and solid-state theorists is the numerical simulation of systems with a severe sign problem. If we denote the microscopic states of a given physical system by , and the thermodynamics of such system is described by a partition function of the form , we say that the system in question presents a sign problem if the “weights” are not real and positive: This implies that we cannot interpret as a proper probability distribution, and the standard, efficient Monte Carlo algorithms cannot be applied. Not all sign problems are equally severe. Let us restrict ourselves for simplicity to the case where the are real but not positive definite111The discussion for complex weights does not add any fundamental difficulty.. One can easily devise a reweighting algorithm that uses the absolute value as the weight of each state, and shifts the sign of into the observables. Now a standard Monte Carlo method is applicable, and in the limit of infinite statistics we should obtain the correct result. With finite statistics, however, a key quantity is the thermodynamic average of the sign of each contribution to the partition function, that is, . If this quantity goes to zero exponentially with the volume, , then we would need an exponential amount (in the volume of the system ) of statistics to get correct results, which is of course impossible in practice. In this case we say that the sign problem is severe.

QCD at finite baryon density, QCD with a topological term in the action, chains of quantum spins with antiferromagnetic interactions, the two-dimensional O(3) non linear sigma model with a topological term, and the Hubbard model are some of the most popular examples of relevant physical systems where a SSP is present. The existence of a SSP is the main reason for the little progress made on the theoretical understanding of these physical systems outside of phenomenological models.

In order to check novel Monte Carlo methods designed to tackle such problems, it is highly desirable to have a set of benchmark calculations as extensive as possible. For very few systems an analytic solution is known, for example, the one-dimensional antiferromagnetic Ising model with an imaginary magnetic field, the two-dimensional compact U(1) model with topological term, or the two-dimensional Ising model with an imaginary magnetic field . In a few other cases the sign problem can be avoided by reformulating the physical system with new degrees of freedom, taking advantage of the fact that a good choice of these degrees of freedom provides an equivalent physical system free from the sign problem, which can therefore be simulated by standard methods.222Unfortunately this idea works only in a few cases which, until now, are not the most interesting physical systems. Indeed none of the examples previously mentioned have been solved with this idea.

Our motivation for this paper is to provide a benchmark calculation for a system for which we do not have an analytic solution available, nor a reformulation that avoids the sign problem. We study the two-dimensional antiferromagnetic Ising model with a purely imaginary magnetic field, which can be thought of as a toy model for the usual physics. Indeed the Euclidean partition function for QCD with a nonvanishing term can be written in the form

 ZV(θ)=∑npV(n)eiθn (1)

where , the topological charge, is an integer, and is, up to a normalization, the probability of the topological sector at . This has the same structure as the partition function of the antiferromagnetic Ising model in an external purely imaginary magnetic field, as we will see in detail later on, and we expect that the SSP in both systems should also be similar.

This system was studied in Matveev and Shrock (2008) by locating the zeros of the partition function in the complex temperature-magnetic field plane, and they find, for purely imaginary magnetic field, a rich phase structure with two phases characterized by a vanishing (nonvanishing) staggered magnetization, separated by a phase transition line. We study this system by an exact cumulant expansion to eighth order, followed by the analytic computation of the partition function and other physical quantities for a large number of degrees of freedom with the help of a standard multiprecision algorithm. This amounts essentially to the computation of the effective Hamiltonian up to order , and therefore is expected to work well in the high-temperature regime, and we provide strong evidence that this is indeed the case. Our results are consistent with Matveev and Shrock (2008), and extend the results of Azcoiti et al. (2011), obtained through the application of algorithms developed in Azcoiti et al. (2002, 2003), and through a mean-field analysis. We are able to obtain a more precise quantitative determination of the transition line separating the paramagnetic and antiferromagnetic phases of the model.

For some systems with a SSP, we know a priori that the partition function will be positive, for example systems in thermal equilibrium with a (Hermitian) Hamiltonian description. Such is the case in a quantum field theory with a term. In the toy model we study here, although we do not have a rigorous proof in this case,333This would imply a nontrivial restriction on the position of the Lee-Yang zeros for the antiferromagnetic Ising model. To the best of our knowledge, very little is rigorously known about such zeros. we have evidence that, at least in the region where the approximation we use is valid, the partition function is indeed positive (it is trivially always real).

Such evidence is twofold. First, we can prove rigorously that up to the fifth cumulant, the partition function is indeed positive. Unfortunately we have not been able to extend this proof to higher cumulants, but in our multiprecision calculations with up to eight cumulants, we have never seen an instance where the partition function is negative or vanishes. This is highly nontrivial: If instead of a constant imaginary magnetic field we try, for example, to put a staggered imaginary field in our lattice (this is of course equivalent to the ferromagnetic model with a constant imaginary field), we immediately get a fluctuating sign for the partition function.

Second, there have been studies locating the Lee-Yang zeros of the two-dimensional antiferromagnetic Ising model up to lattices Kim (2004), and in lattices Matveev and Shrock (2008). Up to that size there is no sign of any zeros cutting the imaginary axis at any temperature.

Whereas this by no means amounts to a rigorous proof, we believe it provides a strong indication that, at least in the region of interest for this paper, this model should have a positive partition function.

This paper is organized as follows. Section II is devoted to formulate the model and to recall the main ingredients and results of the mean-field approximation developed in Azcoiti et al. (2011). In Sec. III we introduce the cumulant expansion, report the analytical results for the first eight cumulants in the two-dimensional model, and write the analytical expressions for the free energy and mean values of interesting physical quantities. The results for the staggered magnetization, susceptibility, and phase diagram of the model are reported in Sec. IV, where we also compare our results at and with the analytical solutions of Onsager (1944); Lee and Yang (1952); McCoy and Wu (1967). In Sec. V we report our conclusions. The technical details of the analytical computation of the cumulant expansion can be found in Appendix A, and several tables with numerical results can be found in Appendix B.

## Ii Two-dimensional Ising model

The Ising model Ising (1925); Onsager (1944); Lee and Yang (1952); McCoy and Wu (1967); Matveev and Shrock (1995); McCoy and Wu (1973); Matveev and Shrock (2008) has been studied for a long time now, and it has known analytical solutions in the one-dimensional case at any external magnetic field Ising (1925), and in two dimensions only for the case without magnetic field Onsager (1944) and for Lee and Yang (1952); McCoy and Wu (1967). The model with a pure imaginary magnetic field suffers from a SSP in any number of dimensions. In addition to that, the expected phase diagram for is non trivial Azcoiti et al. (2011), making the reconstruction of the dependence of the observables even more challenging. All this makes the model a good theoretical laboratory to test new methods designed to deal with the SSP. It is therefore worthwhile to carry out a detailed study of this model at purely imaginary magnetic field, particularly because little progress has been achieved on reconstructing the dependence of the observables, apart from the analysis of Azcoiti et al. (2011) and the recent study in de Forcrand and Rindlisbacher (2016).

The partition function of the model, following the conventions of Azcoiti et al. (2011), is:

 Z=∑{si}exp(F∑sisj+iθ12∑isi). (2)

The half magnetization

 M2≡12∑isi, (3)

is an integer taking any value between and , where is an even number denoting the total number of spins in the lattice. It is in this sense that we identify with a topological charge and regard the imaginary magnetic field term in the action as a term. It is important to mention that, from now on, we will consider only the antiferromagnetic case , since the model with imaginary field does not define a unitary theory for arbitrary values of the ferromagnetic coupling Lee and Yang (1952); Azcoiti and Galante (1999).

As we shall see in detail in the next section, by dividing the rectangular lattice into two sublattices, introducing the respective magnetizations and , making a cumulant expansion and keeping only the first cumulant, we arrive at the following approximation to the partition function (where denotes the dimensionality of the lattice):

 Z1c(F,θ)=∑{si}exp(iθM1+M22+4FdNM1M2). (4)

We recall now the mean-field analysis carried out in Azcoiti et al. (2011). The resulting partition function,

 ZMF(F,θ)=∑{si}exp(iθM1+M22−FdN(M1−M2)2), (5)

is different from Eq. (4). However, it can be seen to give the same qualitative results for the observables and the phase diagram. In this regard, we will consider the first-cumulant expansion as a mean-field approximation to , and the general expansion itself as an improvement of it, at least for small , where the expansion is expected to converge.

Applying standard saddle-point techniques to the mean-field partition function Azcoiti et al. (2011), one obtains the phase diagram shown in Fig. 1. A second order critical line,

 dFc=12cos2θc2, (6)

separates two different phases: a staggered one, with , for , and a paramagnetic one, with , for .

## Iii Cumulant expansion and observables

Our interest is focused on the antiferromagnetic model, where the staggered magnetization is a good order parameter. From now on we will work with a rectangular two-dimensional lattice, although the method is easily generalizable to any number of dimensions. We divide the lattice into two sublattices and in a chessboard fashion. In the two-dimensional lattice this means that if and index, respectively, the row and the column of a given spin, this spin will be in the first (second) sublattice if the sum is even (odd). For simplicity we will require that both lengths of the lattice be even. Denoting by the total number of points in the lattice, we define the magnetization densities and as

 mj≡MjN/2≡∑i∈ΩjsiN/2j=1,2, (7)

and the density of staggered magnetization is

 ms≡m1−m22. (8)

Let us denote by the number of microstates with magnetization densities and in sublattices and , respectively, that is,

 g(m1,m2)=∑{si}δ⎛⎝∑i∈Ω1si−M1⎞⎠δ⎛⎝∑i∈Ω2si−M2⎞⎠. (9)

A trivial computation gives:

 g(m1,m2)=(N/2N1+)(N/2N2+), (10)

with for . Defining now the expected value at fixed as:

 ⟨O⟩m1,m2 ≡ 1g(m1,m2) ×∑{si}δ(∑i∈Ω1si−M1)δ(∑i∈Ω2si−M2)O

we can rewrite the partition function (2) in the form:

 Z = ∑m1,m2g(m1,m2) (12) ×⟨exp(iθ2∑isi+F∑sisj)⟩m1,m2.

The term in Eq. (12) is just , and therefore constant under fixed and ; we can take it out of the expected value, arriving at

 Z = ∑m1,m2g(m1,m2)e14Niθ(m1+m2) (13) ×⟨exp(F∑sisj)⟩m1,m2.

We cannot evaluate exactly the expectation value in Eq. (13), as that would be equivalent to solving exactly the model for arbitrary values of the external field. Instead we perform a cumulant expansion and truncate at a given order. Let us recall the definition:

 ⟨etX⟩≡exp(∞∑n=1κntnn!), (14)

where the th cumulant is an th degree polynomial in the first noncentral moments of , given by the following recursion formula:

 κn=μ′n−n−1∑m=1(n−1m−1)κmμ′n−m,μ′n≡⟨Xn⟩. (15)

By expanding in cumulants in our partition function, taking and , we obtain

 Z = ∑m1,m2g(m1,m2) ×exp(14Niθ(m1+m2)+∞∑n=1κn(m1,m2)Fnn!),

where now the moments are given by

 μ′n=⟨(∑sisj)n⟩m1,m2. (17)

The computation of these quantities is somewhat involved, and we relegate the details to Appendix A. We calculate the cumulants using a numerical (but exact) method, up to . The results, at leading order in 444We can calculate the subleading terms also, but they become irrelevant as we approach the thermodynamic limit., for , are

 κ1 = 2Nm1m2, κ2 = 2N(m12−1)(m22−1), κ3 = 8Nm1m2(m12−1)(m22−1), κ4 = 4N(21m12m22−9(m12+m22)+5) ×(m12−1)(m22−1), κ5 = 32N(51m12m22−39m12−39m22+31) ×m1m2(m12−1)(m22−1), κ6 = 64N(675m14m24−690[m14m22+m12m24] +705m12m22+75[m14+m24−m12−m22]+8) ×(m12−1)(m22−1), κ7 = 128N(10935m41m42−13950[m41m22+m21m42] +3375[m41+m42]+17760m21m22−4290[m21+m22] +1051)m1m2(m21−1)(m22−1). κ8 = 32N(1685565m16m26−2604735[m16m24 (18) +m14m26]+994455[m16m22+m12m26] −55125[m16+m26]+4026645m14m24 −1541085[m14m22+m12m24]+85575[m14+m24] +595077m12m22−33663[m12+m22]+2125) ×(m21−1)(m22−1)

Now we can compute an approximation to the expected value of any observable of the form as follows:

 ⟨O⟩ = 1Z∑m1,m2O(m1,m2)g(m1,m2) ×exp{iθM1+M22+nmax∑n=1Fnn!κn(m1,m2)},

where depends implicitly on the number of cumulants included in the approximation, , and on the number of spins of the system . Taking the limit of both and to infinity, we should recover the exact result in the thermodynamic limit. Using this technique, we have computed several observables, such as the density of free energy , the density of internal energy , the specific heat and both the usual and the staggered magnetization and , respectively. The precise definitions of the computed observables are the following:

 ϕ ≡ −1NFlogZ, (20) e ≡ −12NdlogZdF,cv≡−F2ddFe, (21) ⟨m⟩ ≡ ⟨m1+m22⟩,⟨ms⟩≡⟨m1−m22⟩. (22)

It must be noted that at , where the model has an analytical solution, the free energy has a singularity at Lee and Yang (1952); McCoy and Wu (1967). In the next section we will talk about its nonsingular part, which is simply the result of subtracting the singular term from the full expression:

 ϕ≡ϕns−12Flog(1−e4F). (23)

As we have mentioned before, the complex-valued exponentials in Eq. (LABEL:sumExponentials) give rise to a severe sign problem. To deal with it we use a multiprecision algorithm, which allows us to keep as many digits as needed. In order to crosscheck our calculations we have used several multiprecision libraries (GMP, GNU MPFR, GNU MPC, gmpy2) to do the sum over and . The computational cost when computing the observables grows on one hand with due to the number of summands in (LABEL:sumExponentials). In addition to that, the number of digits needed grows linearly with , increasing the cost of each multiprecision operation.

## Iv Results

At and we know the analytical solution for the two-dimensional Ising model Onsager (1944); Lee and Yang (1952); McCoy and Wu (1967), and therefore we can compare the exact results with the approximations obtained from Eq. (LABEL:sumExponentials). We can see in Figs. 2 and 3 the density of free energy as a function of the coupling , for different approximations. Concretely we show the approximations obtained by keeping only the first, up to the fourth, and up to the eighth cumulant. For clarity we show only the results corresponding to the largest size that we have calculated, although we have carefully checked that the finite-size effects are tiny at that value of . We can see that the agreement with the exact result, especially for the fourth and eighth approximations, is excellent at small , where we can expect the cumulant expansion to be well behaved. At the approximations start to drift away from the analytic result, especially the eighth, possibly indicating the lack of convergence of the cumulant expansion at such larger couplings.

The above results are consistent with those of the density of internal energy, which we can see in Figs. 4-6. The same can be said about the specific heat for , in Fig. 7. The results of the specific heat for , in Fig. 8, show also a good agreement with the analytical solution, as long as we are far from the critical point. In the neighborhood of the critical point we can see that keeping a finite number of cumulants has a strong impact. However, the results seem to converge to the exact solution quickly when we increase the number of cumulants, and indeed the peak when including all eight cumulants is not far from the analytic result.

The agreement with the exact results both at and at suggests that the cumulant expansion can be trusted at all values of , as long as .

We expect a nonvanishing value of to signal the transition from the paramagnetic to the staggered phase. Because of translational symmetry, we cannot simply compute this observable, since for a finite system it is always zero [permuting and leaves Eq. (LABEL:sumExponentials) invariant]. However, we can compute , which also separates the weak and strong coupling phases.

In Fig. 9 we show results for at . One can see how, as we approach the thermodynamic limit, becomes a steeper function of . To obtain the critical line for a given cumulant approximation, we numerically calculate the quantity (which should diverge in the thermodynamic limit at the critical line), and find the maximum along lines of constant . This gives us, for each size and each value of , . We can see in Fig. 10 the behavior of such quantity as a function of and , for the specific value , in the eight cumulant approximation. The height of the peak does not scale as , at least at the volumes we have been able to calculate, therefore suggesting a continuous phase transition; however, our data are not extensive enough to calculate the critical exponents.

The phase diagram obtained in this way is shown in Fig. 11, for several truncation orders of the cumulant expansion. The transition lines that we obtain lie entirely below , where we have good evidence that the cumulant expansion works well. The change from the line corresponding to and is very large, but the results seem to stabilize quickly with the order of the expansion, and the lines corresponding to and are quite close together. Therefore we expect the phase diagram for to be a quite accurate approximation to the exact one. Further evidence of this is the agreement with the few maximal values for estimated in Matveev and Shrock (2008) from the computation of the zeros of the partition function of the model in the complex temperature-magnetic field plane. As can be seen in the plot, they lie above but quite close to our line.

As another crosscheck we show in Fig. 12 results for the specific heat at in the eight cumulant approximation, computed for several system sizes. The behavior is similar to the one in Fig. 10: a peak of increasing height in the vicinity of the critical point, and smooth behavior and small finite N effects elsewhere.

## V Conclusions

We have analyzed the two-dimensional antiferromagnetic Ising model with an imaginary magnetic field by analytical techniques. We have calculated the first eight cumulants of what is essentially the expansion of the effective Hamiltonian in powers of the inverse temperature, and computed physical quantities for a large number of degrees of freedom with the help of multiprecision algorithms. The motivation for such a calculation was to have an example of a physical system with SSP and nontrivial phase structure, the dynamics of which is well known, at least in the high-temperature region.

Our results confirm the qualitative picture described in Azcoiti et al. (2011), and predict the existence of two phases in this model, which can be characterized by the staggered magnetization as an order parameter. The finite-size scaling suggests that the two phases are separated by a continuous phase transition line. The position of the critical point at is in very good agreement with the exact result , and the free and internal energy densities at agree also well with the analytical prediction, at least in the high-temperature regime, thus giving reliability to our results in this region. Therefore this model could be a good laboratory to check proposals to simulate physical systems afflicted by a SSP.

###### Acknowledgements.
This work was funded by Ministerio de Economía y Competitividad/Fondo Europeo de Desarrollo Regional Grants No. FPA2012-35453 and No. FPA2015-65745-P and Diputación General de Aragón-Fondo Social Europeo Grant No. 2015-E24/2.

## Appendix A Computation of the cumulants κn

In order to use expressions (LABEL:ZenKappas) and (LABEL:sumExponentials), we need to compute the cumulants . The th cumulant can be calculated in terms of the first noncentral moments ,

 μ′n≡⟨(∑sisj)n⟩m1,m2, (24)

by means of the recursion relation (15). The summation over runs over each couple of neighboring spins, or in other words, over each link. Two neighboring spins always belong to different sublattices.

Before going further, let us comment on two intermediate results. First, we consider a lattice of spins, the magnetization of which is the sum , and ask about the expected value of the product of of these spins at fixed (or fixed , the number of positive spins), that is, . One can perform this calculation by means of the microcanonical formalism, arriving at

 ⟨s1s2⋯sn⟩m=1(NN+)n∑k=0(−1)k(nk)(N−nN+−n+k). (25)

In the above expression, can be read as the number of negative spins in the product . In this way, the first summand, , counts the number of states with zero negative spins in the product and multiplies it by the expected value of the product in this case, . The second one, , does the same for one negative spin in , and so on. Dividing the sum by the total number of configurations with magnetization , one obtains the previous expected value at fixed . Secondly, consider an observable in our two sublattice system, with a dependence on and such as we can write it as . In this case, from the definition (LABEL:defmicro) of the expected value at fixed and , we have

 ⟨O1(m1)O2(m2)⟩m1,m2=⟨O1(m1)⟩m1⟨O2(m2)⟩m2. (26)

This immediately applies to the spin product . We can always divide it into two products and , each one containing the spins of one of the sublattices, and then

 ⟨s1s2⋯sn⟩m1,m2=⟨sa⋯sb⟩m1⟨sα⋯sβ⟩m2. (27)

With the previous couple of results, we come back to Eq. (24), and apply the linearity of the expected value, arriving at

 μ′n=∑,,⋯,⟨sisjsksl⋯spsq⟩m1,m2, (28)

which is the sum of the expected values of the product of links, running over all permutations with repetitions of these links. Then, in every summand we have the product of spins, in some cases with some of them identical. Taking into account that , each summand can be reduced to the expected value of the product of different spins, and being the number of spins in each sublattice. Since by means of Eq. (25) we already have an expression that computes , the problem is reduced to count how many summands in Eq. (28) have spins. We call these numbers geometrical factors, and denote them by . Following this convention, we can write the th central moment as

 μ′n=∑{n1,n2}G(n1,n2)⟨sa⋯sbn1spins⟩m1⟨sα⋯sβn2spins⟩m2, (29)

where the sum runs over the couples of integers the sum of which is even and less than or equal to .

The computation of the geometrical factors can be done by hand for the first few cumulants. As an example, for the second noncentral moment we have to compute four cases: the two links being the same (sharing both spins), sharing only one spin belonging to the first or the second sublattice, and finally not sharing any spin at all. That is, in terms of the previous notation,

 {(n1,n2)}={(0,0),(2,0),(0,2),(2,2)}. (30)

The factors can be computed easily in this case, even for an hypercubic lattice of arbitrary dimension , arriving at the following expression for the second moment

 μ′2 = Nd⟨1⟩+Nd(d−1)(⟨s1s2⟩m1+⟨s1s2⟩m2) +Nd(Nd−2(d−1)−1)⟨s1s2⟩m1⟨s1s2⟩m2.

We can use this expression to calculate the second cumulant ,

 κ2=μ′2−μ′21N→∞−−−−→Nd(m21−1)(m22−1), (32)

where we have taken the thermodynamic limit, keeping only the terms of order , which is the leading order for all cumulants. Subleading orders can be preserved if needed, but they are not relevant for our paper. The difficulty of the previous computation escalates quickly with the order of the cumulant, and it is quite cumbersome for just . In order to get beyond this limitation, we have developed a program which computes the geometrical factors numerically for a finite bidimensional lattice. Since these factors are polynomials in of order (and with integer coefficients), we can run the program for lattices of different sizes, obtaining a set of (,) points, which we can use to recover the exact integer coefficients of each geometrical factor, by means of the Lagrange interpolation formula.

The basic idea of the program is very simple. We just construct a periodic rectangular lattice, with , being the order of the cumulant we want to compute. With this restriction we avoid products of links crossing the entire lattice, that would not appear in the thermodynamic limit for any finite cumulant. Once we have this, we start a loop running over all the permutations with repetitions of links, and perform the following steps,

• We have a product of links, or equivalently spins, .

• Recursively, we remove couples of equal spins from this product.

• We classify the remaining product by the number of spins in each sublattice, .

• We add one to the geometric factor and proceed to the next iteration.

When the algorithm finishes, we obtain all the values for a given . The computational cost is associated to the number of iterations of the main loop, which grows as , that is, exponentially with the order of the cumulant. In practice, we have only reached the computation of the fourth cumulant with this program. However, a number of optimizations can be implemented in order to reach higher order cumulants, which we summarize in what follows.

### a.1 Translational symmetry

Our lattice is symmetric under translations, implying that all geometrical factors are proportional to , the number of links. Fixing, e.g., the first link of the product, one obtains the same , but divided by a common factor . The same factor is gained in the overall speed of the program. In addition to that, the degree of the polynomials is also reduced by one, and it suffices with (instead of ) different sizes in order to recover the dependence. One can go even further by realizing that the geometrical factor corresponding to non-neighboring links, , is the only one with maximum degree . This allows us to express it in terms of the remaining factors,

 1NdG(n,n) = (Nd)n−1 (33) −1Nd∑{(n1,n2)}∖(n,n)G(n1,n2),

which are only of order or less. This means that it is enough to run the program for lattice sizes, compute all the geometrical factors but via the Lagrange interpolator, and then with the previous expression find the dependence of this last factor.

### a.2 From permutations to combinations

The product of links commutes, so its contribution to the geometrical factors is the same regardless of the order. Then, we can change the main loop over permutations with repetition to a loop over combinations with repetition, by taking into account the multiplicity of each combination. Schematically, we perform

 ∑i,j,…,k contrib(lilj⋯lk) (34) →∑i≤j≤⋯≤kmult×% contrib(lilj⋯lk),

where contrib represents a function in our program that takes a product of links and returns the contribution to the geometrical factors. If there are different links, each one appearing times, the multiplicity of the combination is given by

 mult=n!k1!⋯kr!. (35)

### a.3 Blocks - Grouping links together

Many of the link products have few, if any, repeated spins, and their contributions to the geometrical factors can be counted without having to analyze one by one each of them. This is possible by grouping them in sets of links that we will call in what follows blocks, and replacing the loop over link products by a loop over block products. When the blocks in a product are not neighbors (i.e., they do not have any common spin), we do not need to perform the computation link by link and the contribution can be summed up trivially. Let and be two non-neighboring blocks, each one composed by links, and let us denote the contributions to the geometrical factors by , where is an integer counting how many products of links have () spins in the first (second) sublattice. Then we have

 contrib(b1b3)=N2b(2,2), (36)

or in general, for the product of non-neighboring blocks, . Following this strategy, we divide our lattice into unidimensional blocks of links, in a way that the th block, , contains all links the first spin of which belongs to the th column. As a consequence, is a neighbor of blocks and , and, taking into account the boundary conditions, and are neighbors too.

When we have a product of neighboring blocks, we proceed as before, analyzing the link products one by one, and there is no computational saving. But when the blocks are not neighbors, we move from iterations to a single one.

### a.4 Clusters of blocks

The block method, as defined above, fails to save any computation time if two or more blocks are neighbors in a given block product. However, we can extend the method by dividing each block product into several subproducts, which we will denote as clusters. In each cluster, one can always connect one block to another by the equivalence relation of being neighbors (sharing spins). And in the same way, in each product different clusters never share any spin. This allows us to compute the contributions of each cluster separately, and then compose them with the following law,

 λ1(a,b)⊕λ2(c,d)=λ1λ2(a+c,b+d). (37)

If the contributions of the clusters involve more than one geometrical factor, linearity applies,

 ∑abλab(a,b) ⊕ ∑cdλcd(c,d)= (38) ∑ab,cdλabλcd(a+c,b+d).

Processing one cluster with blocks takes a computing time proportional to . So dividing the whole block product in smaller clusters implies for almost every block product a significant amount of time saved. Only when all the blocks are part of the same cluster there is no speed up.

Another major optimization can be performed by realizing that translational invariance can also be applied here, since a given cluster, say , and any of its translations, , have the same contribution to the geometrical factors. Then, when a cluster is going to be computed, we can express it in terms of its equivalence class, compute its contribution, and store it in memory. Every time one of its translations appears, we just take the value from the memory, saving a lot of computing time. In addition to that, once we have computed the factors for the first size , we know in advance all the cluster contributions for any lattice (the blocks keep its size constant). Since almost all the computing time is spent in figuring out the cluster contributions, we reduce in this way the full problem of computing the geometrical factors in lattices of different sizes to only one size, the smallest one, . In practice, the time spent by the rest of the sizes needed is barely the of that of the first size.

### a.5 Computation of a cluster

The last optimization concerns the computation of the clusters themselves. Until now it is done simply by performing a loop over each possible permutation of links belonging to each of the blocks in the cluster. However, one can go one step further and divide the blocks composing the cluster into smaller sets, that we will call sites. A site is simply the set of two links the first spin of which lies in the site , that is,

 site(i,j)≡{sijsi+1,j,sijsi,j+1}. (39)

With this new subdivision, we can apply in the same way the techniques described above. In order to compute the cluster , we start a loop over every permutation of sites , with . Each site product is divided into clusters, the contributions of which can be summed with Eq. (38) and are calculated by performing another loop over each link product ( iterations for a site product of elements). Finally, by summing up each site product contribution, we obtain the whole cluster contribution.

All the described optimizations do not remove the exponential dependence on of the algorithm. However, they allow us to reach the eighth cumulant, which takes about three days of computing time in a modern laptop.

## Appendix B Numerical tables

In this appendix we present some of the data corresponding to the figures in Sec. IV. In addition to that, we provide numerical results for several observables at and .

## References

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters