A Geometric Monte Carlo Algorithm for the Antiferromagnetic Ising model with"Topological"Term at $\theta=\pi$

In this work we study the two and three-dimensional antiferromagnetic Ising model with an imaginary magnetic field $i\theta$ at $\theta=\pi$. In order to perform numerical simulations of the system we introduce a new geometric algorithm not affected by the sign problem. Our results for the $2D$ model are in agreement with the analytical solutions. We also present new results for the $3D$ model which are qualitatively in agreement with mean-field predictions.


Introduction
Since its introduction many years ago, the Ising model has been a prototype statistical system for studying phase transitions and critical phenomena [1]. With the advent of the epoch of computer numerical simulations to study statistical systems, this model has become even more important as a test bench to develop new algorithms. There are many interesting physical systems for which, due to the sign problem, we do not have efficient numerical algorithms. Some examples include QCD at finite density or with a non-vanishing θ term. This situation has hindered progress in such fields for a long time, and it is thus of great interest to study novel simulation algorithms. In the present work we develop and test a new algorithm, which belongs to a class of "geometric" algorithms [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], and which is applicable to the D ě 2 dimensional antiferromagnetic Ising model with an imaginary magnetic field iθ (see [21] and [22]) at θ " π, with which we are able to solve the sign problem that afflicts this model when using standard algorithms. Preliminary results were presented in [24]. This paper is organized as follows. In Section 2 we introduce the Antiferromagnetic Ising model. In Section 3 we derive our geometric representation of the model. Section 4 is devoted to the construction of the numerical algorithm. In Section 5 we present the numerical results and finally Section 6 contains our conclusions. Some technical details on the ergodicity of the algorithm as well as on the numerical analysis are contained in the appendices.

The Antiferromagnetic Ising model with a topological term
We consider the Ising model in D ě 2 dimensions, defined on a hypercubic lattice Λ with an even number of sites L " 2n in each direction, and with either open or periodic boundary conditions. The Hamiltonian of the model is (1) Here the spin variables are s x "˘1, and the sum ř px,yqPB is over the pairs of sites px, yq that are nearest neighbors; we denote the set of all such pairs by B. Moreover, J is the coupling between nearest neighbors, and B is an external magnetic field. The reduced Hamiltonian H " H{pk B T q, where T is the temperature and k B the Boltzmann constant, is written as with F " J{pk B T q, h " 2B{pk B T q. As the total number of spins is L D " p2nq D , and therefore even, the quantity Q " 1 2 ř x s x is an integer number, taking values between´L D {2 and L D {2. Q can then be thought of as playing the role of a topological charge. It is then worth studying what happens for imaginary values of the reduced magnetic field h, i.e., for h " iθ. The topological charge Q is odd under the Z 2 transformation s x Ñ´s x @x: while at θ " 0 the system is symmetric under this transformation, for θ ‰ 0, π the Z 2 symmetry of the system is explicitly broken. At θ " π the contribution of the topological charge to the Boltzmann factor amounts to i.e., this contribution is Z 2 invariant, and therefore the Z 2 symmetry is restored; it has to be checked if it is spontaneously broken or not.
3 Geometric representation of the model 3

.1 Partition function
We will now introduce a geometric representation for the model at h " iθ " iπ. Let us rewrite the partition function of the system, where we have taken into account that the volume of the system, V " L D , is a multiple of 4, and that coshp´xq " coshpxq, sinhp´xq " sinhpxq. Note now the following symmetry of this partition function. Define the two staggered lattices Λ p1,2q as follows: Λ p1q " tx " pi 1 , . . . , i D q P Λ | pi 1`¨¨¨`iD qmod 2 " 0u , Λ p2q " tx " pi 1 , . . . , i D q P Λ | pi 1`¨¨¨`iD qmod 2 " 1u .
Nearest neighbors always belong to different staggered lattices; as a consequence, if we change variables in the sum in Eq. (4) by changing the sign of all the spins in only one of the two staggered lattices, say, s x Ñ´s x @x P Λ p2q , Eq. (4) becomes (6) Therefore, ZpF, θ " πq " Zp´F, θ " πq, i.e., at θ " π the ferromagnetic and antiferromagnetic models are essentially equivalent. In conclusion, we can write ZpF, θ " πq " Let us denote by B! the power set of B, i.e., A subset b can be seen as a configuration of "active bonds" (we will sometimes also refer to active bonds as "dimers") between neighboring sites. Let N rbs be the number of elements of b, i.e., the number of active bonds; clearly,N rbs " N rBs´N rbs is the number of inactive bonds. Finally, we define the quantity π x rpy, zqs " and let π x rbs be the number of bonds in b that "touch" x, π x rbs " ÿ py,zqPb π x rpy, zqs .
Armed with this notation, we can rewrite the product over pairs of neighboring sites in Eq. (7) as follows, ź px,yqPB rcoshp|F |q`sinhp|F |qs x s y s " The sum over spin configurations in Z vanishes unless each spin appears an even number of times in the product above, and gives a factor of 2 per spin otherwise, that is, Summarizing, This is the geometric representation that we will use in our algorithm.

Observables
It is useful to generalize the partition function to the case of variable couplings, i.e, H rts x u, tF xy u, θ " πs "´ÿ px,yqPB This allows us to calculate all the correlation functions for an even number of spins xs x 1 s x 2 . . . s x 2k y, 1 by taking derivatives with respect to F xy for an appropriate set of px, yq. Indeed, choosing a set of paths C 1 , . . . C k connecting the spins pairwise (there are no restrictions on these paths), and then performing derivatives with respect to all the pairs px, yq appearing in those paths (if a pair appears mpx, yq times, one has to take the mpx, yq-th derivative with respect to the corresponding coupling), , . -ˇˇˇˇt The geometric representation for the partition function with variable couplings is similar to the one obtained for constant coupling, and the final result takes the form Taking derivatives with respect to F x 1 y 1 , . . . , F x l y l , and carrying them inside the summation over b, one obtains an extra factor cothpF where δ b px, yq " 1 if px, yq P b, and 0 otherwise. Noting that (18) and setting F xy " F , @px, yq, we finally obtain Finally, denoting by ∆rb; tx i , y i us " ř l i"1 δ b px i , y i q, and dividing by the partition function, we obtain where the symbol xx. . .yy indicates the average taken with the probability distribution 2 P pbq " tanhpF q N rbs ÿ bPB!, tπxrbsmod 2"1@xu tanhpF q N rbs .
As an example, let us write down the two-point correlation function xs x s y y for x and y lying on the same lattice axis, e.g., y " x`l1. Let C be a path connecting x and y on the lattice; the simplest choice is a straight-line path C " Y l i"1 px i , y i q with y i " x i`1 , x i " x`pi´1q1 (i " 1, l`1), x 1 " x and y l " y. We have Cpl, F q " xs x s x`l1 y " xxtanhpF q l´2∆rb;tx i ,y i us yy , where ∆rb; tx i , y i us is the number of active bonds in configuration b along the straight-line path connecting x and x`l1. However, we stress the fact that the specific choice of the path is irrelevant, as they are all equivalent, as long as the endpoints are fixed. We can also obtain expressions for bulk observables, such as the energy density and the specific heat, by taking derivatives of the partition function with respect to the coupling constant. Using the geometric representation it is easy to see that such observables are directly related to the average number of active bonds and its fluctuations. We denote by B " N rbs the total number of active bonds, H 0 is defined as and we denote the fluctuations in these quantities by ∆H 0 " H 0x H 0 y, ∆B " B´xxByy. Then we have the following relations for the energy density ε and the specific heat c V (we only show for simplicity the relation in the case of periodic boundary conditions): where D is the dimensionality of the system. For small values of the coupling F the energy density and the specific heat are approximately equal to the average occupation number of the dimers and to its fluctuations:

Monte Carlo algorithm
In order to perform calculations by means of Monte Carlo methods, we need an efficient algorithm to explore the space of configurations, that as we have seen is given, in the geometric representation, bỹ B " tb P B! | π x rbs mod 2 " 1 @xu , or, in words, by those configurations for which the number of active bonds touching any site x of the lattice is odd. We will call the setB the set of admissible configurations. For simplicity we will describe in detail the algorithm only for D " 2, as it can be easily generalized to D ą 2.
As far as numerical simulations are concerned, it is enough if we know at least one admissible configuration, and a set of updating rules that take us from one admissible configuration to another and that satisfy detailed balance and ergodicity.
To describe the updating rules we found convenient to work with the dual lattice (equivalently, the set of squares of the original lattice). Let us denote a point in the dual lattice by x˚, and let us enumerate the four links of the corresponding square in anticlockwise order, as shown in Fig. (1). To describe a configuration of bonds, we introduce for each site of the dual lattice a vector with four components, Apx˚q " pA 1 px˚q, A 2 px˚q, A 3 px˚q, A 4 px˚qq, defined such that A i px˚q " 1 if the corresponding bond is active and A i px˚q " 0 if it is inactive. The bond configuration is completely specified by the dual lattice field Apx˚q, although this description is redundant: indeed, one has that A 3 px˚q " A 1 px˚`1q and A 2 px˚q " A 4 px˚`2q. In Table 1 we show all the possible configurations of a given square on the lattice, i.e., a point x˚in the dual lattice. The graphs Spx˚q corresponding to these configurations, and the number of active bonds, wpx˚q " Suppose now that we are given an admissible configuration, and we want to update it to a new admissible configuration. This requires updating some bonds by changing their state, i.e., A i px˚q Ñ 1´A i px˚q.  sequence, the parity of the number of active bonds touching this site would change from odd to even under the update, and the resulting configuration would not be admissible. The most general admissible update consists therefore in changing the state of bonds belonging to a closed (possibly disconnected) path. This is easily seen to be equivalent to perform consecutive updates on the set of elementary squares that cover that part of the lattice enclosed by the path, 3 Table 2: Transformation of configurations under conjugation C. Acted upon a graph, we denote the transformation with a hat. The quantity ∆wpx˚q is the variation in the number of active bonds when passing from Spx˚q toĈSpx˚q. AsĈ is an involution,ĈrĈSpx˚qs " Spx˚q, so the conjugate of a configuration in column 3 is the corresponding configuration in column 1, and the variation of the number of active bonds changes sign when passing from a configuration in column 3 to the corresponding configuration in column 1.
where C stands for conjugation, and where we have introduced the vector I " p1, 1, 1, 1q. Under conjugation, A i px˚q Ñ 1´A i px˚q, i " 1, . . . 4; clearly, C 2 " I, where I is the identity, IApx˚q " Apx˚q. The variation ∆wpx˚q in the number of active bonds under conjugation is given by ∆wpx˚q " (30) In Table 2 we show the pairs of configurations of an elementary square connected by conjugation, together with ∆wpx˚q. These updating steps can be applied independently to all the sites of the dual lattice, and are clearly reversible. We have now the ingredients to set up the first version of a Metropolis algorithm (Algorithm 1).

Algorithm 1
1. At a given site x˚of the dual lattice, compute ∆wpx˚q corresponding to a conjugation step.

Repeat the procedure for all the dual lattice sites.
It is easy to see that this algorithm satisfies detailed balance. The only question remaining is that of the ergodicity of the algorithm.

Ergodicity
Concerning ergodicity the situation is slightly different for open or periodic boundary conditions.
The simplest case is that of open boundaries. In this case we can prove (A.1) that all the admissible configurations can be transformed to the same configuration through a sequence of conjugation moves. As these transformations are reversible, any admissible configuration is connected to all the others. Therefore in this case algorithm 1 is ergodic and can be used as is to simulate the system.
The case of periodic boundary conditions is slightly more complicated. One can see that in this case, the total number of vertical (respectively horizontal) active bonds modulo 2 on any given row (respectively, column) of the dual lattice is conserved under the updating moves, and moreover is the same for any row (column). Calling these numbers vertical and horizontal parity, P V and P H , respectively, this defines four different "sectors" of admissible configurations, classified by parities (P V , P H ) P tp0, 0q, p1, 0q, p0, 1q, p1, 1qu. It can be shown that the updating moves are ergodic within each sector separately (A.2). Therefore we have to modify slightly algorithm 1 to obtain an ergodic algorithm (algorithm 2).
The case of higher dimensionality is a straightforward generalization of the algorithms presented here, obtained by applying algorithm 1. Start from a configuration of parity, say, p0, 0q.
2. After a certain number of sweeps through the whole lattice, using algorithm 1, propose a change of the parity P H of the configuration, by proposing the inversion of the state of all the horizontal bonds on a row i " i 0 , j " 1, . . . L, of the direct lattice, chosen randomly, so (possibly) passing to a configuration of parity p0, 1q with probability tanhp|F |q ∆w , where ∆w " and we have denoted B H pi, jq " A 1 pi, jq, and CB H pi, jq " 1´B H pi, jq.
3. After the same number of sweeps through the whole lattice, (try to) change the parity P V of the configuration by proposing the inversion of the state of all the vertical bonds on a column j " j 0 , i " 1, . . . L, of the direct lattice, chosen randomly, so (possibly) passing to a configuration of parity p1, 1q, if the proposal in the previous point has been accepted, or p1, 0q, if the proposal in the previous point has been rejected, with probability tanhp|F |q ∆w , where now ∆w " and we have denoted B V pi, jq " A 4 pi, jq, and CB V pi, jq " 1´B V pi, jq.

Iterate the procedure.
First of all, we tested the algorithms in the case of periodic boundary conditions. The acceptance rate of the global update varies widely with both the coupling and the size of the system. For example, for the 2D model at F "´1.0 the acceptance rate drops from about 60% at L " 16 to about 30% at L " 64; for a smaller coupling, F "´0.6, the acceptance rate drops from 27% at L " 16 to 3% at L " 64. 5 Then we checked that the bounds on the number of active bonds in each configuration of our system ( B) and also on the average of this quantity ( C) were respected. For simplicity we only show in Figs. (2) the average number of active bonds.
Then we calculated the correlation functions (22) as well as the energy density (24) and the specific heat (25) both for the two-dimensional and three-dimensional systems. We have evaluated these quantities for different values of the coupling F and various volumes, with lattice sizes ranging from L " 16 to L " 1024. Simulations were done collecting 100k measurements for each value of F . We discarded between 10k and 20k configurations at the beginning of each run in order to ensure thermalization. The data analysis was done using the jackknife method over bins at different blocking levels.
In Figs. (3), (4) and (5) we show how the correlation functions depend on the distance d, both for the 2D and 3D models, choosing different antiferromagnetic couplings F ă 0 and varying the lattice volume V . In Fig. (6) we show the staggered and the standard magnetization squared in the 2D model for two values of the coupling F , together with a solid line indicating the analytical result. As can be seen the staggered magnetization squared is always different from zero while the standard magnetization vanishes for all values of the coupling F , and the results are in perfect agreement with the analytical solution of Refs. [21][22][23]. In Fig. (7) we show the corresponding results for the 3D model at the same values of the coupling, as well as the mean-field prediction for the staggered magnetization obtained in Ref. [25]. Again the standard magnetization vanishes, whereas the staggered one does not, and its value is close to the mean-field prediction for large values of the coupling |F |. Therefore, despite the vanishing of the standard magnetization, the Z 2 symmetry of the sys- tem is spontaneously broken both in the 2D and in the 3D models. We can also notice an apparent decrease of Cpd, F q at large d for small values of the coupling F both for the 2D and 3D models. This is due to the heavy-tailed distributions of the correlators, which are also responsible for the noisy behavior seen at large d. In Fig. (8) we show the probability distributions of the logarithm of the correlators for a lattice size L " 64 and for F "´0.4 and F "´2.0. Clearly we notice that for a small coupling |F | the values are spread in a wider range than for F "´2.0; also a long tail develops for large distances d, thus making more difficult a precise evaluation of the correlators. 6 Finally, in Figs. (9) and (10) we show the energy density ε and the specific heat c V as a function of |F | for various lattice sizes L. The solid lines present in the figures for the 2D model are the analytical results for these quantities [22]. We can see that the energy density and the specific heat do not show any sign of singular behavior in F for the range of couplings studied.

Conclusions
In this paper we studied the 2D and 3D antiferromagnetic Ising model with a "topological" θ-term at θ " π. For this model we introduced a new geometric algorithm free from the sign problem. The numerical part of the work has been devoted to testing the algorithm for the two-dimensional model against known analytical results, with which we obtain perfect agreement, and then afterwards to study the three-dimensional system. Our findings strongly support the scenario that, despite the vanishing of the standard magnetization, the staggered magnetization is non-zero for all D ě 2, and therefore the Z 2 symmetry is spontaneously broken for all values of F at θ " π.
It would be interesting to study whether it is possible to introduce a "worm" in our algorithm, in the spirit of [2]. This could change the dynamics of the system; in particular an implementation that allows the worm to wind through the lattice might be able to tunnel between the different parity sectors. We leave the study of such a possibility for a future work.      Table 3, starting from the upper-right corner of the dual lattice, proceeding downward along a column of the dual lattice, and then moving to the column to the left. We call this transformation "reduction", which we denote with R. It is easily seen from Table 3 that it coincides with the identity if A 2 px˚q " 0, and with conjugation if A 2 px˚q " 1, i.e., In the case of open boundary conditions, the sites of the dual lattice are x˚pi, jq with i, j P t1, . . . , N´1u. By construction, after R has been applied to the right-most column of the dual lattice, the resulting configuration will not have any vertical bond on the right-hand side of this column, i.e., A 2 px˚q " 0 for x˚pi, N´1q, i P t1, . . . , N´1u, and only horizontal bonds will be present. As R transforms an admissible configuration into another admissible configuration, the only possibility is that all the horizontal bonds of the right-most column of the dual lattice are active, as the rightmost sites of the direct lattice need to be touched by at least one active bond and they cannot have more than one (see Fig. (11)). We now apply R to the following column: as the vertical bonds move to the left, and all the sites in the before-last column of the direct lattice already have an active bond, the horizontal bonds in the before-last column of the dual lattice must be inactive    (see Fig. (12)). If we now repeat the procedure, we find ourselves as after the first step: there is a column of sites of the direct lattice with no vertical bonds, and no horizontal bonds to the right, exactly as the right boundary of the lattice. Therefore, the result iterates for pairs of columns of the dual lattice, until we reach the left boundary. As all the sites are already connected horizontally to the right, and as the uppermost site can have at most a single vertical bond, it is easy to see that no vertical bond can appear. We have then reduced the initial configuration to the reduced configuration of Fig. (13). As no reference has been made to the specific form of the initial configuration, the procedure applies equally to any admissible configuration, which completes the proof of ergodicity for open boundary conditions. This also provides an admissible configuration, so completing the construction.  Figure 13: Reduced configuration, after completing the reduction process.

A.2 Periodic boundary conditions for the 2D model
Consider a 2D NˆN square lattice with periodic boundary conditions. In this case, the sites of the dual lattice are x˚pi, jq with i, j P t1, . . . , N u, and are closed as well by periodic boundary conditions. We start by applying R to column N of the dual lattice. Again, after this procedure there are no vertical bonds active on the righthand side of the column; however, this time the right-most sites of the direct lattice can have horizontal active bonds entering from the left or from the right, due to the periodic boundary conditions (see Fig. (14)). There are 2 N possibilities, according to the position of the horizontal bonds (left or right). Repeating the procedure for column N´1, one immediately sees that the active horizontal links must be the same as in column 1, as there are no other possibilities due to the absence of vertical bonds in column N´1 of the direct lattice (see Fig. (15)). Repeating again, we find for column N´2 the same configuration of horizontal bonds of column N , and so on. Since there is now an even number of columns on the dual lattice, the final pattern will be N {2 identical pairs of columns on the dual lattice, with all the odd columns equal to each other, and all the even columns equal to each other.
At first sight, one could think that there are therefore 2 N reduced configurations instead of one (see Fig. (16)). At second sight, one could think that there are even more, as it is possible that all the vertical bonds in column 1 of the direct lattice are active (see Fig. (17)): in this case all the sites in this column have two more active bonds (vertical bonds are however absent in the rest of the lattice, by construction). However, one can show that many of these configurations, that we call almost-reduced, can be transformed into one another: by direct inspection, one can show that a configuration with two horizontal bonds on the right of two adjacent sites can be transformed into the configuration with two bonds on the left of the same two adjacent sites, all the rest unchanged (and vice-versa); and that a configuration with one bond on the right and one on the left of two adjacent sites can be transformed into the configuration with the first bond on the left and the second one on the right of the same two adjacent sites, all the rest unchanged. In the case without vertical bonds, it suffices to apply a conjugation to a site of the dual lattice belonging to the row corresponding to the two given sites, and then repeat the reduction procedure; in the case with vertical bonds, it is enough to repeat the reduction procedure along the whole lattice (see Fig. (18)).
Exploiting these equivalences, one can "swap" pairs of bonds, finally reducing to one of the four configurations of Fig. (19). However, it is not possible to transform one of these configurations into one another by means of our admissible moves. To see this, define the number of active vertical bonds on row i of the dual lattice, and the number of active horizontal bonds on column j of the dual lattice, Observing that the admissible updates always change A 1 pi, jq`A 3 pi, jq and A 4 pi, jq`A 2 pi, jq by 0 or˘2, one has that the vertical parity P V piq " N V piq mod 2 and the horizontal parity P H pjq " N H pjq mod 2 are conserved under admissible moves. It is evident that for each of the four configurations of Fig. (19), P V piq " P V @i and of P H pjq " P H @j, due to the horizontal periodicity of the configuration, and moreover that these configurations have different values of (P V , P H ). As these quantities are conserved, we have that P V piq " P V @i and P H pjq " P H @j also for a generic configuration, which under R will be transformed to the reduced configuration with the same pair pP V , P H q.

A.3 Ergodicity in the 3D model
Consider a 3D cubic NˆNˆN lattice, with N even, and a configuration of bonds satisfying the constraint that every site is touched by an odd number of bonds. We define plaquette update (PU) the process of inverting the state (active/inactive) of all the bonds surrounding an elementary square of the lattice. We denote by P µν pnq the square formed by the links pn, n`μq, pn`μ, n`μ`νq, pn`μ`ν, n`νq and pn`ν, nq, with µ, ν " 1, 2, 3 and n " pn 1 , n 2 , n 3 q where n i " 1, 2, . . . , N .
We now show that by means of PU's it is possible to transform any admissible configuration into any other in the case of open boundary conditions (obc), while in the case of periodic boundary conditions (pbc) it is necessary to supplement these transformations with a few global transformations. In this way we can construct an ergodic algorithm. The strategy is to reduce any admissible configuration to one and the same, "elementary" configuration.
The proof is as follows. Consider the squares lying in the p1, 2q planes of the lattice. Fix n 1 " N´1, and @ n 2 , n 3 apply a PU on      Figure 19: The four inequivalent reduced configurations for periodic boundary conditions. The pair pP V , P H q is reported above each configuration.
As a result, we have now a configuration where there are no active bonds in the directions 2 and 3, except possibly at n 1 " 1. Therefore, in the bulk of the lattice the constraint on the number of bonds touching a site has to be enforced by means of bonds in direction 1, and therefore only one bond can touch a site. This immediately implies the following relation: B 1 pn 1 , n 2 , n 3 q " 1´B 1 pn 1`1 , n 2 , n 3 q , n 1 " 1, . . . , N´2 , (36) for n 2 , n 3 " 1, . . . , N . For obc, since B 1 pN´1, n 2 , n 3 q " 1 @ n 2 , n 3 in order for the sites at n 1 " N to satisfy the constraint, this implies obc : B 1 pn 1 , n 2 , n 3 q " mod pn 1 , 2q , n 1 " 1, . . . , N´1 , @ n 2 , n 3 .
(37) Since all sites at n 1 " 1 are touched by at least one bond, active bonds in directions 2 and 3 at n 1 " 1 must form closed paths, so that they contribute an even number of active bonds to all sites. Any non-selfintersecting path of active bonds can be made inactive by performing a PU on all the plaquettes contained in the path; self-intersecting paths of active bonds can be decomposed in non-self-intersecting paths that have no link in common, and so also in this case all bonds in direction 2 and 3 can be made inactive, i.e., we obtain B 2,3 pn 1 , n 2 , n 3 q " 0 @ n 1 , n 2 , n 3 . This is the sought-after "elementary configuration", to which all other configurations can be reduced in the case of obc.
For pbc, in order for the sites at n 1 " N to satisfy the constraint, one has furthermore that B 1 pN´1, n 2 , n 3 q " 1´B 1 pN, n 2 , n 3 q, @ n 2 , n 3 , and so B 1 p1, n 2 , n 3 q " 1´B 1 pN, n 2 , n 3 q, which implies that a single bond in direction 1 touches the sites at n 1 " 1. This again implies that the active bonds in directions 2 and 3 at n 1 " 1 must form closed paths. While for closed paths that do not wind around the lattice the considerations made above apply, so that the corresponding bonds can be made inactive, this is not true for winding paths. However, we can basically repeat the same strategy used above: starting from n 2 " N´1, perform a PU in the plaquettes P 23 p1, n 2 , n 3 q if the rightmost bond (pn`2, n`2`3q) is active, for all n 3 , and then repeat the procedure for all n 2 until we reach n 2 " 1. At this point, there can be only bonds along direction 2, and possibly bonds along direction 3 at n 2 " 1. Enforcing the constraint that every site has to be touched by an even number of bonds lying in the plane n 1 " 1, we conclude that B 2 p1, n 2 , n 3 q " B 2 p1, 1, n 3 q @ n 2 , i.e., active bonds (if any) in direction 2 form closed straight-line paths winding around the lattice. As a consequence, the number of 2-bonds touching a site at n 2 " 1 is even (either zero or 2), so that the same must apply for the 3-bonds at n 2 " 1, and so again B 3 p1, 1, n 3 q " B 3 p1, 1, 1q @ n 3 . Furthermore, it is easy to see that the straight lines formed by the 2-bonds can be parallelly shifted, and that a pair of such straight lines at n 3 and n 3`1 can be made inactive, by means of PU's. In conclusion, it is always possible to reduce the configuration of 2-and 3-bonds on the plane n 1 " 1 to one of the following cases: 1 : B 2 p1, n 2 , n 3 q " 0 , @ n 2 , n 3 , B 3 p1, n 2 , n 3 q " 0 , @ n 2 , n 3 ; 2 : B 2 p1, n 2 , 1q " 1 @ n 2 , B 2 p1, n 2 , n 3 q " 0 @ n 2 , @ n 3 ‰ 1 , B 3 p1, n 2 , n 3 q " 0 , @ n 2 , n 3 ; 3 : B 3 p1, 1, n 3 q " 1 @ n 3 , B 3 p1, n 2 , n 3 q " 0 @ n 2 ‰ 1, @ n 3 , B 2 p1, n 2 , n 3 q " 0 , @ n 2 , n 3 ; 4 : B 2 p1, n 2 , 1q " 1 @ n 2 , B 2 p1, n 2 , n 3 q " 0 @ n 2 , @ n 3 ‰ 1 , B 3 p1, 1, n 3 q " 1 @ n 3 , B 3 p1, n 2 , n 3 q " 0 @ n 2 ‰ 1, @ n 3 .
The configurations 1´4 are characterized by the number (0 or 1) of 2-bonds and 3-bonds in a strip at fixed n 2 and n 3 , respectively. Notice that these numbers do not depend on which strip we choose. Since there are no other 2-bonds and 3-bonds in the rest of the lattice, these numbers are the same if we count the 2-bonds in a slice (i.e., all n 1 and n 3 ) at fixed n 2 , and if we count the 3-bonds in a slice (i.e., all n 1 and n 2 ) at fixed n 3 . Again, these numbers do not depend on the chosen slice. Since a PU in the pµ, νq plane does not change the parity of the number of µ-or ν-bonds in a slice at fixed n µ or n ν (although changing possibly their number), we can determine to which of the above configurations in the n 1 " 1 plane a given generic configuration can be reduced, by simply computing these parities, which again do not depend on which slice we choose. Defining P 2 " mod p ÿ n 1 ,n 3 B 2 pn 1 , n 2 , n 3 q, 2q , P 3 " mod p ÿ n 1 ,n 2 B 2 pn 1 , n 2 , n 3 q, 2q , (39) one can easily classify the configurations 1´4 (see Table 4).
The last step is to simplify as much as possible the configuration of 1-bonds. For each pn 2 , n 3 q, the configuration of 1-bonds is entirely determined by the value of B 1 p1, n 2 , n 3 q. It is easy to see that by means of PU's, we can simultaneously change the configurations at pn 2 , n 3 q and pn 2`1 , n 3 q or pn 2 , n 3`1 q, i.e., we can go from B 1 p1, n 2 , n 3 q " b 1 , B 1 p1, n 2`1 , n 3 q " b 2 (b i " 0, 1) to B 1 p1, n 2 , n 3 q " 1´b 1 , B 1 p1, n 21 , n 3 q " 1´b 2 , or from B 1 p1, n 2 , n 3 q " b 1 , B 1 p1, n 2 , n 3`1 q " b 2 to B 1 p1, n 2 , n 3 q " 1´b 1 , B 1 p1, n 2 , n 3`1 q " 1´b 2 . Using this observation, we can change the position of those rows characterized by B 1 p1, n 2 , n 3 q " 0, and trade a pair of such rows for a pair with B 1 p1, n 2 , n 3 q " 1. For definiteness, we move them first towards n 2 " 1 at fixed n 3 , removing them when two show up at neighboring sites: as a result, B 1 p1, n 2 , n 3 q " 1 @ n 2 ‰ 1, n 3 . Next, we move them towards n 3 " 1 at fixed n 2 " 1, again removing them when two show up at neighboring sites: as a result, B 1 p1, n 2 , n 3 q " 1 @ n 2 ‰ 1, n 3 ‰ 1. Then, two possibilities remain: either B 1 p1, 1, 1q " 1 or B 1 p1, 1, 1q " 0. This means that in the first case there is an even number of 1-bonds in any slice at fixed n 1 , while in the second case this number is odd. Since, as mentioned above, the parity of the number of 1-bonds in a slice at fixed n 1 is not changed by a PU, the "reduced configuration" of 1-bonds corresponding to any given generic configuration is determined by the value of the quantity P 1 " mod p ÿ n 2 ,n 3 B 1 pn 1 , n 2 , n 3 q, 2q .

B Determination of the upper and lower number of bonds permitted in a given configuration (pbc case)
Due to the constraints, the number of dimers touching a given site must be odd for an admissible configuration. In dimension D and using periodic boundary conditions, this means that this number is 1, 3 . . . , 2D´1. We denote by V k the number of vertices with k dimers in a given configuration, and by V the total number of sites. One has the two following relations: where B is the total number of dimers. The first equation above can be rewritten as and substituting this into the second equation we get Since the terms under the summation signs are positive, we get the inequalities V ď 2B ď p2D´1qV and dividing by the total number of links DV we obtain