Antiferromagnetic Ising model in an imaginary magnetic field

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 $\theta$ 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 solidstate theorists is the numerical simulation of systems with a severe sign problem. If we denote the microscopic states of a given physical system by s, and the thermodynamics of such system is described by a partition function of the form Z = s P (s), we say that the system in question presents a sign problem if the "weights" P (s) are not real and positive: This implies that we cannot interpret P (s) 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 P (s) are real but not positive definite [? ]. One can easily devise a reweighting algorithm that uses the absolute value |P (s)| as the weight of each state, and shifts the sign of P (s) 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, sign(P (s)) . If this quantity goes to zero exponentially with the volume, sign ∝ e −αV , then we would need an exponential amount (in the volume of the system V ) 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 * Corresponding author: eduroyo@unizar.es 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 twodimensional Ising model with an imaginary magnetic field h = iπ/2. 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. [? ] 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 where n, the topological charge, is an integer, and p V (n) is, up to a normalization, the probability of the topological sector n at θ = 0. 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 [1] by locating the zeros of the partition function in the complex temperaturemagnetic 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 T −8 , 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 [1], and extend the results of [2], obtained through the application of algorithms developed in [3,4], 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, [? ] 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 14 2 lattices [5], and in 12 × 13 lattices [1]. 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 [2]. 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 inter-esting 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 h = 0 and iπ/2 with the analytical solutions of [6][7][8]. 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 [1, [6][7][8][9][10][11] has been studied for a long time now, and it has known analytical solutions in the one-dimensional case at any external magnetic field h [9], and in two dimensions only for the case without magnetic field h [6] and for h = iθ/2 = iπ/2 [7,8]. 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 d ≥ 2 is non trivial [2], 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 [2] and the recent study in [12].
The partition function of the model, following the conventions of [2], is: The half magnetization is an integer taking any value between −N/2 and N/2, where N is an even number denoting the total number of spins in the lattice. It is in this sense that we identify M/2 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 F < 0, since the model with imaginary field does not define a unitary theory for arbitrary values of the ferromagnetic coupling [7,13].
As we shall see in detail in the next section, by dividing the rectangular lattice into two sublattices, introducing the respective magnetizations M 1 and M 2 , making a cumulant expansion and keeping only the first cumulant, we arrive at the following approximation to the partition function (where d denotes the dimensionality of the lattice): We recall now the mean-field analysis carried out in [2]. The resulting partition function, 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 firstcumulant expansion Z 1c as a mean-field approximation to Z, and the general expansion itself as an improvement of it, at least for small F , where the expansion is expected to converge.
Applying standard saddle-point techniques to the mean-field partition function [2], one obtains the F − θ phase diagram shown in Fig. 1. A second order critical line, separates two different phases: a staggered one, with m s = 0, for θ ≥ θ c , and a paramagnetic one, with m s = 0, for θ < θ c .

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 Ω 1 and Ω 2 in a chessboard fashion. In the two-dimensional lattice this means that if i and j index, respectively, the row and the column of a given spin, this spin will be in the first (second) sublattice if the sum i + j is even (odd). For simplicity we will require that both lengths of the lattice be even. Denoting by N the total number of points in the lattice, we define the magnetization densities m 1 and m 2 as and the density of staggered magnetization is Let us denote by g(m 1 , m 2 ) the number of microstates with magnetization densities m 1 and m 2 in sublattices Ω 1 and Ω 2 , respectively, that is, (9) A trivial computation gives: with N j+ ≡ N (1 + m j )/4 for j = 1, 2. Defining now the expected value at fixed m1, m2 as: we can rewrite the partition function (2) in the form: The θ term in Eq. (12) is just iθ (m 1 + m 2 ) N/4, and therefore constant under fixed m 1 and m 2 ; we can take it out of the expected value, arriving at 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: where the nth cumulant κ n is an nth degree polynomial in the first n noncentral moments of X, given by the following recursion formula: By expanding in cumulants in our partition function, taking t = F and X = s i s j , we obtain where now the moments are given by 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 n = 8. The results, at leading order in Now we can compute an approximation to the expected value of any observable of the form O(m 1 , m 2 ) as follows: where O depends implicitly on the number of cumulants included in the approximation, n max , and on the number of spins of the system N . Taking the limit of both n max and N 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 e, the specific heat c v and both the usual and the staggered magnetization m and m s , respectively. The precise definitions of the computed observables are the following: It must be noted that at θ = π, where the model has an analytical solution, the free energy has a singularity at F = 0 [7,8]. 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: As we have mentioned before, the complex-valued exponentials in Eq. (19) 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 m 1 and m 2 . The computational cost when computing the observables grows on one hand with N 2 due to the number of summands in (19). In addition to that, the number of digits needed grows linearly with N , increasing the cost of each multiprecision operation.

IV. RESULTS
At θ = 0 and π we know the analytical solution for the two-dimensional Ising model [6][7][8], and therefore we can compare the exact results with the approximations obtained from Eq. (19). We can see in Figs. 2 and 3 the density of free energy as a function of the coupling |F |, 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 N that we have calculated, although we have carefully checked that the finite-size effects are tiny at that value of N . We can see that the agreement with the exact result, especially for the fourth and eighth approximations, is excellent at small |F |, where we can expect the cumulant expansion to be well behaved. At |F | 0.57 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 θ = 0, 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 θ = 0 and at π suggests that the cumulant expansion can be trusted at all values of θ, as long as |F | 0.57.
We expect a nonvanishing value of m s 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 N system it is always zero [permuting m 1 and m 2 leaves Eq. (19) invariant]. However, we can compute m 2 s , which also separates the weak and strong coupling phases.
In Fig. 9 we show results for m 2 s at θ = 2. One can see how, as we approach the thermodynamic limit, m 2 s becomes a steeper function of |F |. To obtain the critical line for a given cumulant approximation, we numerically calculate the quantity d dθ m 2 s (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 N and each value of θ, F c (θ). We can see in Fig. 10 the behavior of such quantity as a function of F and N , for the specific value θ = 2, in the eight cumulant approximation. The height of the peak does not scale as N , 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 |F | = .45, where we have good evidence that the cumulant expansion works well. The change from the line corresponding to k = 1 and 4 is very large, but the results  seem to stabilize quickly with the order of the expansion, and the lines corresponding to k = 4 and 8 are quite close together. Therefore we expect the phase diagram for k = 8 to be a quite accurate approximation to the exact one. Further evidence of this is the agreement with the few maximal values for F c estimated in [1] 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 k = 8 line.
As another crosscheck we show in Fig. 12 results for the specific heat at θ = 2 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 [2], 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 θ = 0 is in very good agreement with the exact result F c = log(1 + √ 2)/2 ≈ 0.4407, 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. In order to use expressions (16) and (19), we need to compute the cumulants κ n . The nth cumulant can be calculated in terms of the first n noncentral moments µ n , by means of the recursion relation (15). The summation over < ij > 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 N spins, the magnetization of which is the sum m = i s i , and ask about the expected value of the product of n of these spins at fixed m (or fixed N + , the number of positive spins), that is, s 1 s 2 · · · s n m . One can perform this calculation by means of the microcanonical formalism, arriving at s 1 s 2 · · · s n m = 1 N N+ n k=0 (−1) k n k N − n N + − n + k .
(A2) In the above expression, k can be read as the number of negative spins in the product s 1 s 2 · · · s n . In this way, the first summand, k = 0, counts the number of states with zero negative spins in the product s 1 s 2 · · · s n and multiplies it by the expected value of the product in this case, (−1) 0 = 1. The second one, k = 1, does the same for one negative spin in s 1 · · · s n , and so on. Dividing the sum by the total number of configurations with magnetization m = 2N + /N −1, one obtains the previous expected value at fixed m. Secondly, consider an observable O(m 1 , m 2 ) in our two sublattice system, with a dependence on m 1 and m 2 such as we can write it as O 1 (m 1 )O 2 (m 2 ). In this case, from the definition (11) of the expected value at fixed m 1 and m 2 , we have (A3) This immediately applies to the spin product s 1 s 2 · · · s n . We can always divide it into two products s a · · · s b and s α · · · s β , each one containing the spins of one of the sublattices, and then s 1 s 2 · · · s n m1,m2 = s a · · · s b m1 s α · · · s β m2 . (A4) With the previous couple of results, we come back to Eq. (A1), and apply the linearity of the expected value, arriving at µ n = <ij>,<kl>,··· ,<pq> s i s j s k s l · · · s p s q m1,m2 , (A5) which is the sum of the expected values of the product of n links, running over all permutations with repetitions of these links. Then, in every summand we have the product of 2n spins, in some cases with some of them identical. Taking into account that s 2 i = 1∀i, each summand can be reduced to the expected value of the product of n 1 + n 2 different spins, n 1 and n 2 being the number of spins in each sublattice. Since by means of Eq. (A2) we already have an expression that computes s 1 · · · s n m , the problem is reduced to count how many summands in Eq. (A5) have (n 1 , n 2 ) spins. We call these numbers geometrical factors, and denote them by G(n 1 , n 2 ). Following this convention, we can write the nth central moment as where the sum runs over the couples of integers (n 1 , n 2 ) the sum of which is even and less than or equal to n.
The computation of the geometrical factors G(n 1 , n 2 ) can be done by hand for the first few cumulants. As an example, for the second noncentral moment µ 2 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, {(n 1 , n 2 )} = {(0, 0), (2, 0), (0, 2), (2, 2)}. (A7) The factors G(n 1 , n 2 ) can be computed easily in this case, even for an hypercubic lattice of arbitrary dimension d, arriving at the following expression for the second moment We can use this expression to calculate the second cumulant κ 2 , where we have taken the thermodynamic limit, keeping only the terms of order O(N ), 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 n of the cumulant, and it is quite cumbersome for just n ≥ 4. In order to get beyond this limitation, we have developed a program which computes the geometrical factors G(n 1 , n 2 ) numerically for a finite L × L bidimensional lattice. Since these factors G(n 1 , n 2 ) are polynomials in N of order ≤ n (and with integer coefficients), we can run the program for lattices of n + 1 different sizes, obtaining a set of (N ,G(N )) 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 L × M lattice, with L, M > n, n 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 n links, and perform the following steps, • We have a product of n links, or equivalently 2n spins, s 1 · · · s 2n .
• Recursively, we remove couples of equal spins from this product.
• We classify the remaining product by the number of spins in each sublattice, (n 1 , n 2 ).
• We add one to the geometric factor G(n 1 , n 2 ) and proceed to the next iteration.
When the algorithm finishes, we obtain all the G(n 1 , n 2 ) values for a given N = LM . The computational cost is associated to the number of iterations of the main loop, which grows as (LM ) n , 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.

Translational symmetry
Our lattice is symmetric under translations, implying that all geometrical factors are proportional to N d, the number of links. Fixing, e.g., the first link of the product, one obtains the same G(n 1 , n 2 ), but divided by a common factor N d. The same factor is gained in the overall speed of the program. In addition to that, the degree of the polynomials G(n 1 , n 2 ) is also reduced by one, and it suffices with n (instead of n + 1) different sizes in order to recover the N dependence. One can go even further by realizing that the geometrical factor corresponding to non-neighboring links, G(n, n), is the only one with maximum degree N n−1 . This allows us to express it in terms of the remaining factors, which are only of order n−2 or less. This means that it is enough to run the program for n−1 lattice sizes, compute all the geometrical factors but G(n, n) via the Lagrange interpolator, and then with the previous expression find the N dependence of this last factor.

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(l i l j · · · l k ) → i≤j≤···≤k mult × contrib(l i l j · · · l k ), (A11) 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 r different links, each one appearing k 1 , . . . , k r times, the multiplicity of the combination is given by

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 b 1 and b 3 be two non-neighboring blocks, each one composed by N b links, and let us denote the contributions to the geometrical factors by λ(n 1 , n 2 ), where λ is an integer counting how many products of links have n 1 (n 2 ) spins in the first (second) sublattice. Then we have or in general, for the product of k non-neighboring blocks, N k b (k, k). Following this strategy, we divide our lattice into unidimensional blocks of 2M links, in a way that the jth block, b j , contains all links the first spin of which belongs to the jth column. As a consequence, b j is a neighbor of blocks j−1 and j+1, and, taking into account the boundary conditions, b 0 and b L−1 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 n blocks are not neighbors, we move from (N d) n iterations to a single one.

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, If the contributions of the clusters involve more than one geometrical factor, linearity applies, Processing one cluster with k blocks takes a computing time proportional to (N d) k . 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 b 0 b 1 b 1 , and any of its translations, b 0+t b 1+t b 1+t , 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 G(n 1 , n 2 ) for the first size L × M , we know in advance all the cluster contributions for any L × M 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 n − 1 different sizes to only one size, the smallest one, M × M . In practice, the time spent by the rest of the sizes needed is barely the 1 − 2% of that of the first size.

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 i, j, that is, With this new subdivision, we can apply in the same way the techniques described above. In order to compute the cluster b 1 . . . b k , we start a loop over every permutation of sites s 1 . . . s k , with s i ∈ b i . Each site product is divided into clusters, the contributions of which can be summed with Eq. (A15) and are calculated by performing another loop over each link product (2 k iterations for a site product of k 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 n 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 θ = 2 and k = 8.