Solution of the sign problem in the Potts model at fixed fermion number

We consider the heavy-dense limit of QCD at finite fermion density in the canonical formulation and approximate it by a 3-state Potts model. In the strong coupling limit, the model is free of the sign problem. Away from the strong coupling, the sign problem is solved by employing a cluster algorithm which allows to average each cluster over the Z(3) sectors. Improved estimators for physical quantities can be constructed by taking into account the triality of the clusters, that is, their transformation properties with respect to Z(3) transformations.


Introduction
At present, the non-perturbative properties of QCD at finite baryon density can not be studied from first principles, because the theory suffers from a fermion sign problem as soon as the baryon chemical potential µ is nonzero. More precisely, the fluctuating sign of the fermion determinant in the path integral at nonzero µ prevents the interpretation of the Boltzmann factor as a probability measure and hence renders importance sampling methods such as Monte Carlo (MC) simulations of lattice QCD inapplicable. One common approach to the problem is to include the fluctuating sign of the Boltzmann factor in the measured observables. While this is in principle a valid procedure, in practice it is bound to fail because the severe cancellations require the statistics to grow exponentially with the space-time volume of the lattice.
Nevertheless, it turns out that in some limiting cases the fermion sign problem is mild or even absent. One such situation concerns lattice QCD in the so-called heavy-dense limit at infinitely strong coupling. In the heavy-dense limit, the fermion contributions to the path integral are encoded in Polyakov loops associated with static quarks and antiquarks propagating in time. In the strong coupling limit, the Polyakov loops decouple in the effective action and the global Z(3) symmetry of QCD at zero density is promoted to a local one. This in turn allows to project the propagating quark states in the canonical formulation exactly onto mesonic and baryonic states localized on single sites. The contributions of those states to the path integral have a phase which is zero or extremely small, hence allowing MC calculations. Despite this happy state of affairs in the strong coupling limit, as soon as the coupling is tuned away from this limit, the sign problem strikes back with vengeanceit swamps away any coherent contribution from the colour neutral states and consequently renders MC simulations impossible.
However, the physical picture emerging in the strong coupling limit indicates how the relevant fermionic contributions to the path integral, being physical and having small or even zero phases, may look like away from infinitely strong coupling. In fact, one may speculate that whenever the quarks and antiquarks form mesonic and baryonic clusters, within which the Polyakov loops fluctuate in conjuction, a coherent physical contribution should arise, very similar to what happens in the strong coupling limit. In order to substantiate this speculation and study the mechanism in more detail, in this paper we simplify the problem further and consider a system in which the gauge dynamics is replaced by that of the Z(3) Potts model and the Polyakov loops are represented by the Potts spins z ∈ Z(3) [1,2,3]. For the canonical formulation of this model, where the quarks are represented in terms of quark occupation numbers, we are able to devise a cluster algorithm that solves the fermion sign problem completely, in a similar way as the cluster algorithm in [4] for the grand-canonical formulation with a finite chemical potential. However, in the canonical situation described in this paper, the solution has a clear and very intriguing physical interpretation in terms of the quark occupation numbers, very much in line with the physical picture that emerges for QCD in the heavy-dense and strong coupling limit as described above.
In the cluster formulation it turns out that the only physical, nonzero contributions to the canonical partition functions are exclusively from clusters which contain exactly a multiple of three quarks. The weights of these clusters are invariant under Z(3) transformations and we denote the corresponding clusters as triality-0 clusters. Clusters with one or two additional quarks are denoted as triality-1 and triality-2 clusters, respectively, refering to the transformation properties of their weights with respect to Z(3) transformations. Their total contributions to the canonical partition functions exactly vanish, but they contribute to observables which consist of nonzero-triality quark fields.
The clusters which are formed in the canonical picture allow for a straightforward physical interpretation. The triality-0 clusters containing a nonzero number of quarks represent baryon or multi-baryon states. While the quarks can move freely within the cluster, of course always respecting the Pauli exclusion principle, they are nevertheless confined within the clusters. It is therefore obvious to interpret these clusters as baryon bags. Empty clusters on the other hand represent vacuum fluctuations of the underlying gauge degrees of freedom. In case we also allow antiquarks in the system, the triality-0 clusters may also contain mesons and represent all kinds of multi-baryon-multi-meson bags.
This physical picture is also useful to interpret the transition from the confined phase to the deconfined one. At low temperature and density, the clusters tend to be very small and form a dilute gas of hadrons, representing the system in the confined phase. When the temperature is increased, the clusters grow larger and start to percolate, at which point the system undergoes a first order phase transition into the deconfined phase. At even higher temperatures, one single cluster essentially fills the whole volume and the quarks can hence move freely within the whole volume.
Before moving on to work this picture out in full detail, it is useful to discuss the symmetry properties of the Z(3) Potts model at finite density. While the model at zero fermion density is invariant under global Z(3)-symmetry transformations, this is no longer true at finite density. In that case the introduction of the chemical potential breaks the Z(3) symmetry explicitely. In contrast, in the canonical formulation the contributions from the static charges have definite transformation properties under Z(3) transformations, such that the contributions vanish unless the number of charges is a multiple of 3. As a consequence, the non-vanishing contributions are manifestly invariant under global Z(3)-symmetry transformations.
The paper is organised as follows. Section 2 contains the definition of the Z(3) Potts model and its canonical formulation. In section 3 we introduce the bond formulation and describe the cluster algorithm for the canonical formulation, while in section 4 we construct improved estimators for various physical quantities such as the quark-antiquark, the quarkquark, the antiquark-antiquark and the 3-quark correlator, as well as the free energy of a single quark and of a single antiquark. A practical algorithm to sample the various sectors in configuration space based on the calculation of the multiplicities of quark configurations with subsequent reweighting is presented in section 5. Next, in section 6 we discuss the severity of the sign problem, while in section 7 we present a selection of simulation results, including an illustration of how the expected phase transition at very low fermion density manifests itself in the canonical setup. Finally, section 8 contains our conclusions, while some technical details are collected in the appendices.

The 3-dimensional Z(3) Potts model
In order to introduce the 3-dimensional Z(3) Potts model we follow [4] and approximate the grand-canonical partition function of QCD in the heavy-dense limit by where the Potts spins z x ∈ Z(3) are defined on discrete lattice points x and stand for the Polyakov loop variables generated by the infinitely heavy quarks. The parameter h = (2κe µ ) β refers to the hopping parameter κ, the chemical potential µ and the inverse temperature β of the original theory. In the canonical formulation the partition function for N Q quarks can be written in the generic form Here the quark occupation numbers n x ∈ {0, . . . , n max } are collected in the quark configuration n = {n x } and the sum is over all quark configurations which fulfill n ≡ ∑ x n x = N Q . The finite maximal number of quarks n max on each site implements the Pauli exclusion principle. The fermionic weights g[z x , n x ] can in principle be derived from QCD and in general also involve the anti-Polyakov loops z * x and corresponding occupation numbersn x for the antiquarks. Here we replace the fermionic weights by For this case, the connection to the grand-canonical partition function in eq.(1) can be made explicit and it can be shown that the correspondence is exact for h ≪ 1, i.e. at low densities. We refer to Appendix A for further details.
For the action S[z] representing the pure gluon action we choose the standard nearestneighbour Potts model interaction where the sum is over all nearest-neighbour lattice sites. Note that the total action including the fermionic weights is manifestly complex and suffers from a sign problem just as in the grand-canonical formulation. However, the properties of the action under global Z(3) transformations ensure that In the limit γ → 0, corresponding to the strong coupling limit in QCD, the Potts spins fluctuate independently on each site and the global Z(3) transformation is promoted to a local one which hence enforces n x = 0 mod 3, ∀x ('strong coupling limit' γ → 0).
The fermionic weights then become trivial and the partition function simply counts the number of allowed quark configurations, In case the quark occupation number is restricted to n x ≤ n max = 3, the number of allowed quark configurations yields G(N Q 3) = V N Q 3 .

Bond formulation and cluster algorithm
For finite gauge coupling the canonical partition function reads and with the help of e γ⋅δz x,zy = 1 bxy=0 δ zx,zy δ bxy,1 (e γ − 1) + δ bxy,0 it can be written as Dz ⟨xy⟩ δ zx,zy δ bxy,1 (e γ − 1) + δ bxy,0 where the sum is over all bond configurations b = {b xy }. For future reference we introduce a double bracket notation to abbreviate the sum over {b}, {n} and Dz weighted by the bond factors. For a generic quark number N Q = ∑ x n x and a generic function i.e., eq.(10) becomes Let us now consider for a moment the case N Q = 0 when there are no fermionic contributions through the Polyakov loops. From the expression above we read off that for a given spin configuration a bond b xy is occupied with probability p(b xy = 1) = (1 − e −γ ) if the two neighbouring spins z x and z y are aligned, else it is empty (b xy = 0). The weight of such a bond configuration is then simply W ({b}) = (e γ − 1) N b , where N b = ∑ ⟨xy⟩ b xy is the total number of occupied bonds, and it is independent of the orientation of the spins within the connected clusters of bonds. Hence, all spins within a cluster can be flipped into any of the three states independently of the spins in the other clusters. Up to here, this is just the well-known Swendsen-Wang cluster reformulation [5] of the Potts model. One can now integrate out the spins for a given bond configuration by summing over all possible combinations of cluster orientations, yielding where N C is the number of clusters for the given bond configuration. This is just the random cluster model representation of the Potts model as derived by Fortuin and Kasteleyn in [6]. Note that the integration over the spins is possible because the bond variables decouple the spin orientations of the clusters from each other, so that the integration over the spin degrees of freedom can be factorized over the clusters. This factorization is in fact the basis for the solution of the sign problem when the fermionic contributions are included.
In order to simulate the system according to the weights given in eq.(13), the bonds can be updated with a local algorithm as follows. A bond whose value does not change the number of clusters is activated with probability p = 1 − e −γ . A bond which connects two clusters C 1 and C 2 , which would otherwise be separate, is a bridging bond and therefore called a bridge in short. Activating a bridge decreases the number of clusters N C by one, and it is therefore done with a probability while the bridge is deactivated with the probability Hence, the update of a bond requires one to check whether the bond under consideration is a bridge or not. This is a well known, difficult problem in computational complexity theory and goes under the name fully-dynamic connectivity problem. For further details and our strategy to deal with the problem we refer to Appendix B.
For N Q > 0 we can include the fermionic contribution in our considerations by constructing an improved estimator for it. For an individual configuration the contribution ∏ x z nx is in general complex, but its average over the subensemble of the 3 N C configurations related by the Z(3) transformations of the individual clusters is rather simple. In the following we denote this subensemble average by ⟨ ⋅ ⟩ 3 N C . We first observe that the total weight can be factorized into individual cluster weights W 0 (C), where by ⟨ ⋅ ⟩ 3 we denote the average over the three Z(3) orientations of a given cluster. This average can be calculated due to the fact that all spins in the cluster are aligned, Denoting the fermion number content of each cluster modulo 3 by the cluster weight simply becomes W 0 (C) = δ n C ,0 . In the following we refer to n C as the triality of the cluster C. Consequently, all bond clusters can be classified according to their fermion content n C , i.e. their triality. Clusters with nonzero triality n C = 1 or 2 have a vanishing contribution to the partition function, 1 while clusters with zero triality n C = 0 contribute positively with weight 1. Hence, the partition function can now be written as where the sums are with respect to all bond and occupation number configurations (with ∑ x n x = N Q ). The fermionic contribution ∏ x z nx x simply projects onto the sector of configurations containing only triality-0 clusters. Note that although the contributions to the partition function are now all positive, the expectation value of the fermionic contribution with respect to the full, unrestricted ensemble with fixed N Q , can still be exponentially small, because the probability to find configurations without nonzerotriality clusters among all possible cluster and fermion number configurations can be exponentially small in the volume. In fact, it is easy to see that eq.(20) is just the expectation value of the complex phase in the modified real action ensemble and represents a measure for the severity of the sign problem. We refer to Sec. 6 for a more detailed discussion.
Nevertheless, the partition function can in principle be updated by directly operating on the bond variables and the fermion occupation numbers. Updates need to ensure that the system stays in the configuration space containing only triality-0 clusters. The bond update can be implemented as before with the cluster break-up probability in eq.(15) modified to where the δ-functions make sure that no clusters with nonzero triality are generated. Note that this step fulfills detailed balance since the situation involving the reverse step, i.e. connecting two clusters with at least one having nonzero triality, can never occur. Once the bond configuration is updated, the fermion occupation numbers can be updated locally by proposing to shift a fermion from site x to one of its nearest-neighbour sites y. The shift is accepted with probability 1 as long as the fermion remains in the same cluster. Hops to sites which are already saturated are forbidden, as well as hops into a different cluster, since these would leave behind two clusters with nonzero triality. Hence, for two neighbouring sites x and y we have where C x and C y denote the clusters containing the sites x and y, respectively. Additional update steps can of course be envisaged and might be favourable with respect to efficiency of the algorithm. For example, one could let three fermions, i.e. a baryon, hop together. Since this update does not change the triality of the clusters, the hop is also allowed across the clusters.
In a sense, the algorithm realizes a kind of a bag model in which the fermions can move freely within bags defined by the individual clusters, i.e. the fermions are confined to the clusters, but deconfined within the clusters. In the confined phase we expect the clusters to be small, so the quarks within a cluster are confined within the small regions of the clusters. In the limit γ → 0 no bonds are allowed and the clusters consist of single sites only, so the quarks are bound into baryons confined to single sites, just as it happens for QCD in the strong coupling limit. Towards the deconfined phase the clusters proliferate, so the quarks can move freely within large regions and hence form baryons of extended size. In the deconfined phase, the clusters percolate the whole volume, so the quarks can move unrestricted throughout the whole volume and the baryon looses its meaning as a bound state of three quarks. Hence, what we have here is a set of degrees of freedom which captures and describes the physics of the underlying system very naturally.

Improved estimators for physical quantities
The reformulation in terms of bond variables and fermion occupation numbers allows to construct improved estimators for physical observables. First we note that the expectation values of observables with nonzero triality vanish in the ensembles with N Q = 0 mod 3. Hence, for the partition function of a single (anti-)Polyakov loop in the restricted ensemble without nonzero-triality clusters we have which reflects the fact that the free energy of a single quark or antiquark is infinite in the background of integer baryon number, F q N Q 3∈N = Fq N Q 3∈N → ∞ due to Gauss' law. The same holds true for the quark-quark or antiquark-antiquark correlation functions in the background of integer baryon number. In contrast, a single quark or antiquark can exist as a probe in a background with non-integer baryon number. Indeed, the partition function including a single quark or antiquark as a source picks up contributions exclusively from the sectors of configurations which contain nonzero-triality clusters, and similarly for the partition functions including a quark-quark or an antiquark-antiquark correlator.
To be more precise, let us construct an improved estimator for the partition function of a single quark at position x in the background of N Q = 2 mod 3 quarks. For a given cluster configuration and fermion number configuration (with ∑ y n y = N Q − 1 and N Q = 0 mod 3) we denote by C x the cluster containing the site x and calculate the subensemble average according to As before the weights W 0 (C) = δ n C ,0 project onto the triality-0 clusters, while the weight i.e. it is nonzero only if C x is a cluster with triality n Cx = 2. Hence, after integrating out the spin variables, the observable for a single quark at position x in the background of N Q = 2 mod 3 quarks is just The observable simply counts the number of configurations with all clusters having triality n C = 0 except the one containing the site x which has triality n Cx = 2. Similarly, the observable for an antiquark in the background of N Q = 1 mod 3 quarks can be written as which counts the number of configurations with all clusters having triality n C = 0 except the cluster C x with triality n Cx = 1. The first δ-function in eq.(28) is just the weight W 1 (C x ) defined in analogy to eq.(26) by Of course, while the absolute free energy can not be measured directly in a Monte Carlo simulation, the free energy difference can be determined from appropriate ratios of partition functions. The free energy difference between N Q = 0 mod 3 quarks and a quark in the background of N Q − 1 quarks for example can be calculated from corresponding to ⟨z x ⟩ ρ in the thermodynamic limit V → ∞ with ρ = N Q V fixed. The crucial step to simulate this system with the spins already integrated out is the observation that the bond cluster configurations contributing to the numerator and the denominator are the same -the only difference is in the fermion numbers and possible fermion configurations. Hence, we can measure the ratio by simulating both the sector where all clusters have zero triality n C = 0 and the sector with one quark less and exactly one cluster containing the site x with triality n Cx = 2. This can be achieved by setting up a Monte Carlo process which probes the configuration space of clusters and quark configurations with N Q and N Q − 1, the latter with exactly one triality-2 cluster connected to site x. More precisely, we add two update steps which try to remove or add a quark from a random site y depending on whether the system is in the sector with N Q = 0 mod 3 or N Q − 1 quarks. The removal is only accepted when the site y is in the same cluster as x, so for a random site y we have while the balancing acceptance probability for adding a quark is Moreover, in the sector N Q − 1 the triality-2 cluster can only be a broken up in two if the triality-2 cluster remains connected to site x, so the cluster break-up probability is given by We can of course make use of translational invariance of F qx and collect statistics for each site in the triality-2 cluster. One just has to make sure that in the update steps above maximally one triality-2 cluster is present in the system. Denoting by N 1 and N 2 the number of clusters with triality 1 and 2, respectively, this means that the allowed configuration space is extended to include configurations with (N 1 , N 2 ) = (0, 0) and (0, 1) nonzero-triality clusters.
It is straightforward to construct similar update steps for simulating both the sector where all clusters have zero triality n C = 0 and the sector with one additional quark and hence exactly one cluster with triality n C = 1. The allowed configuration space is then extended to configurations with (N 1 , N 2 ) = (0, 0) and (1, 0) nonzero-triality clusters. The relative occurrence of configurations in the two sectors eventually yields the free energy difference between N Q = 0 mod 3 quarks and a single antiquark in the background of N Q + 1 quarks, corresponding to ⟨z * x ⟩ ρ in the thermodynamic limit V → ∞ with ρ = N Q V fixed. It is noteworthy that the probability for the transition between the sectors (1, 0) → (0, 0) and (0, 0) → (0, 1) are suppressed by 1 V because the numbers of configurations in the corresponding configuration spaces scale accordingly. As a consequence, the ratios in eq. (30) and (34) are more and more difficult to measure accurately in numerical simulations towards the thermodynamic limit. This can be alleviated by artificially enhancing and suppressing the transition probabilities between the sectors followed by an appropriate a-posteriori reweighting. A more drastic approach is the following. Since the bond cluster configurations contributing to the numerator and the denominator are the same up to the difference in the fermion numbers, we can measure the ratios by simulating only in the sector (N 1 , N 2 ) = (0, 0) and reweighting to the fermion configurations with (1, 0) and (0, 1) using the corresponding multiplicities. We refer to Sec. 5 for the detailed discussion of this alternative approach.
Observables with nonvanishing expectation values in integer baryon number sectors are for example the quark-antiquark correlator ⟨z x z * y ⟩ N Q =0 mod 3 and the 3-point quark correlator ⟨z x z y z w ⟩ N Q =0 mod 3 . Let us first consider the improved estimator for the quark-antiquark correlator. For a given fermion number configuration with ∑ x n x = 0 mod 3 and fixed cluster configuration we calculate the subensemble average while keeping track whether the sites x and y belong to the same cluster, where we have used that zz * = 1. Expressing the weights of the clusters explicitely in terms of δ-functions, we can write the improved estimator for the correlator z x z * y as Hence, the quark-antiquark correlator receives contributions both from the zero triality sector and from the sector with exactly one cluster with triality n Cx = 2 and one cluster with triality n Cy = 1. 2 Therefore, in a numerical simulation we allow the breakup of a zero triality cluster with ∑ x∈C n x > 0 into the two required nonzero triality clusters, and the probability in eq. (21) to remove a bridge is modified accordingly, In addition, we allow a fermion to hop from the cluster C x with zero triality into the zero triality one C y , thereby generating the two required clusters with nonzero triality. Hence, the probability in eq.(22) for a fermion to hop from site v to site w is modified according to In this way we efficiently probe both the sectors which contribute to the correlator, and we obtain the quark-antiquark potential in the background of N Q quarks as Of course one can generalize the above steps in such a way that the allowed configuration space is extended to include configurations with just (N 1 , N 2 ) = (0, 0) and (1, 1) nonzerotriality clusters independent of x and y. By making use of translational invariance, one can then collect statistics for all possible distances x − y on a given cluster and fermion configuration. However, since the bond cluster configurations contributing to the numerator and the denominator are the same up to the difference in the fermion numbers, we can measure the ratio as in the case of ⟨z x ⟩, namely by simulating in the sector containing triality-0 clusters only and reweighting to the fermion configurations with n Cx = 2 and n Cy = 1. We refer again to Sec. 5 for a more detailed discussion of this approach.
The improved estimators for the quark-quark and antiquark-antiquark correlators can be constructed in an analogous way. Here we only note that the quark-quark correlator is defined exclusively in the background of N Q − 2 = 0 mod 3 quarks, allowing only configurations with (N 1 , N 2 ) = (1, 0) and (0, 2) nonzero-triality clusters. In contrast, the antiquark-antiquark correlator is defined only in the background of N Q + 2 = 0 mod 3 quarks, requiring exclusively configurations with (N 1 , N 2 ) = (0, 1) and (2, 0) nonzero-triality clusters. Consequently, the physical quantity of interest is the free energy difference, or potential energy, between N Q = 0 mod 3 quarks and the two (anti)quarks in the corresponding background, Finally, let us now consider the expectation value of the 3-point quark correlator ⟨z x z y z w ⟩ in the background of N Q = 0 mod 3 quarks. The derivation of the corresponding improved estimator is analogous to the derivations before. When calculating the subensemble average of the 3-point correlator z x z y z w including the canonical fermion weight ∏ v z nv v one has to keep track of which of the three spins belong to the same cluster. They can all belong to the same cluster, or to two or three different ones. The latter two possibilities generate contributions in the sectors with exactly two or three clusters with nonzero triality. More explicitly, one has and we could now write down the expression for the partition function ⟪z x z y z w ⟫ N Q in terms of δ-functions which simply single out the various contributions from the sectors containing up to three clusters with nonzero triality. As before, in the simulation we allow the breakup of a zero triality cluster with ∑ x∈C n x > 0 into two clusters with triality n C = 1 and 2, respectively, as long as the resulting clusters contain the source sites x, y, w as given above. In addition, we also allow a further breakup of the cluster with triality n C = 1 into two clusters with triality n C = 2, so as to allow contributions of the kind when all sources sit in different clusters.
By dividing ⟪z x z y z w ⟫ N Q by Z N Q one obtains an expression for the 3-quark potential in the background of N Q quarks, This expectation value can be measured by probing the extended configuration space which allows for clusters with nonzero triality. Denoting by N 1 and N 2 the number of clusters with triality 1 and 2, respectively, the allowed configuration space is restricted to configurations with (N 1 , N 2 ) = (0, 0), (1, 1) and (0, 3) nonzero-triality clusters. It is straightforward to generalize the previous update steps in such a way that only configurations in that restricted configuration space are generated.
Finally, we note that summing the 3-point quark correlator over all positions x, y, w yields the free energy difference between N Q = 0 mod 3 and N Q + 3 quarks, i.e. the baryon free energy in the background of N Q quarks, This quantity gives a direct link between the baryon density ρ B = N Q 3V defined in the canonical formulation discussed here and the quark chemical potential µ in the grand-canonical formulation through the relation

Multiplicities of quark configurations and reweighting
The classification of the clusters and fermion number configurations according to their triality guarantees the identification of those configurations, for which the contributions cancel exactly in the partition function or in the improved observables, so that we are left with positive contributions only. The positivity of configuration weights and the existence of improved estimators provide a complete solution of the complex action problem.
From a practical point of view, this is however not sufficient. While it is easy to identify the relevant configurations, one also needs a suitable algorithm to efficiently sample all the relevant configurations. The algorithm discussed in the previous section suffers for example from inefficiencies due to the fact that the transition, e.g. from the sector (N 1 , N 2 ) = (0, 0) N Q → (0, 1) N Q −1 when simulating the ratio in eq.(30), is suppressed like 1 V at low density. Similar issues may arise in the transitions (0, 0) ↔ (1, 1) when calculating the quark-antiquark correlation, or (0, 0) ↔ (1, 1) ↔ (0, 3) in the calculation of the 3-point quark correlator. In this section we describe how these difficulties can be circumvented by simulating only the (N 1 , N 2 ) = (0, 0) N Q sector with appropriate reweighting into all the other sectors relevant for the various observables. The approach is based on the fact that the bond configurations contributing to the partition function and the improved observables are the same and the contributions only differ in the multiplicities determined by the allowed quark number configurations.
To be more precise, we recall the partition function in eq. (19) and note that the sum over the quark number configurations for a fixed bond configuration is only affected by the product of constraining δ-functions. Hence, for a given cluster configuration, which is controlled by the bond configuration b, we define as the multiplicity of the bond configuration b, i.e. the number of quark configurations compatible with the constraints. The partition function for N Q quarks then becomes Similarly, the correlation function can be written as where we introduced the number of quark configurations compatible with the new constraints as N xy (N Q , b) ≡ {n} δ Cx,Cy δ n Cx ,0 + (1 − δ Cx,Cy )δ n Cx ,2 δ n Cy ,1 C≠Cx C≠Cy We then have where the latter average is over bond configurations {b} distributed according to the Boltz- The multiplicity of a bond configuration, N (N Q , b), can be computed in two steps: first determine all possible partitions {B} of the total baryon number N B = N Q 3 over the clusters C, and secondly for each of these partitions compute the number of possible arrangements of n quarks within each cluster of volume C , denoted by the multiplicity P(n, C ) in the following. Note that P depends only on the number of quarks and the volume of the cluster, while the shape of the cluster is irrelevant. The multiplicity can be computed using the recurrence relation together with the recurrence stopping condition P(n, 1) = 1 for n ≤ n max and 0 otherwise. The recurrence relation is justified by the fact that the number of ways of distributing n quarks on v sites can be computed by separating a site out and then counting how many configurations have no quark on that site, P(n, v − 1), how many have one quark on that site, P(n − 1, v − 1), etc., and then adding up all these contributions.
To generate all partitions {B} of N B baryons into N C clusters one can use the Algorithm 1 given in Appendix C. Note that the partitions only depend on the number of clusters and are not affected by the sizes of the clusters. However, since the number of partitions grows quickly with the number of clusters and the baryon number we will use stochastic estimators, as detailed further below, in order to evaluate the sum over the baryon partitions.
So for a given baryon partition B = {B C }, where B C denotes the number of baryons in the cluster C, the multiplicity, i.e. the total number of quark arrangements fulfilling the triality-0 constraint, is simply the product of quark multiplicities over all clusters. Denoting with N B the total number of baryons in partition B, that is N B ≡ ∑ C B C , we generate the partitions B such that N B = N Q 3 and we finally have the multiplicity of the bond configuration b, i.e. the number of triality-0 quark partitions for that bond configuration.
For N xy (N Q , b) we need partitions that are similar to the ones for N (N Q , b). When C x = C y we have N xy (N Q , b) = N (N Q , b). For the case when C x ≠ C y the trialities of C x and C y are 2 and 1, respectively. A partition satisfying this condition can be generated from a triality-0 partition by moving a quark from C x to C y . All triality-broken partitions with (N 1 , N 2 ) = (1, 1) and n Cx = 1, n Cy = 2 are generated from triality-0 partitions that allow this move, namely the ones with ∑ v∈Cx n v > 0 and ∑ v∈Cy n v < C y ⋅ n max . Denoting the set of triality-0 partitions with B 0,0 ≡ {B} (N 1 =0,N 2 =0) and the triality-broken ones with where M C = ∑ x∈C n x is the number of quarks in cluster C. Above we stressed that the fraction F does not depend on x and y directly, but rather on the clusters that these sites belong to. As such, for a given bond configuration b we can compute F(C ′ , C ′′ ) for all cluster combinations C ′ , C ′′ and then use it to collect statistics for correlators at all distances. For moderate number of clusters and quarks the calculation can be carried out by listing and summing over all relevant baryon partitions. Unfortunately, the number of such partition grows quickly as the number of quarks and clusters increases. It is therefore advisable to use an estimator for the fraction function, where the average if over triality zero partitions B ∈ B 0,0 distributed according to the Boltzmann factor ∏ C P(M C , C ). This can be easily generated using baryon moves from cluster to cluster, accepted or rejected based on a Metropolis process. Note that we assumed that the function P is zero when the quark number is negative or exceeds saturation.
There are two other two-point functions of interest, namely ⟨z x z y ⟩ and ⟨z * x z * y ⟩ which are defined as with The calculation proceeds along the same steps as with the quark-antiquark correlator, with a modified estimator for the multiplicity fraction where the average is again over triality zero partitions B ∈ B 0,0 as in eq.(53).
Other observables of interest are the 1-point averages ⟨z⟩ and ⟨z * ⟩ defined by with The corresponding multiplicity ratios are again calculated using the estimators Finally, the quark chemical potential µ is given by the partition function ratio This ratio can be computed using the multiplicity ratio average where the average is over bond configurations {b} distributed according to the Boltzmann factor P (b) ∝ (e γ − 1) N b 3 N C N (N Q , b). For a given baryon partition B with B = N Q 3 on a fixed bond configuration {b} we can generate a baryon partition B ′ with B ′ = (N Q + 3) 3 by adding one baryon to any of the clusters. However, this map is not one-to-one, because several B partitions can generate the same B ′ partition by adding one baryon. The exact number of such partitions B is N 0 (B ′ ) which is the number of clusters that have B ′ C ≠ 0. Therefore, we have where B C is the quark partition generated from B by adding 3 quarks to cluster C. Thus, the multiplicity ratio can be computed using the estimator 6 Severity of the sign problem A measure for the severity of the sign problem is given by the expectation value of the complex phase of the full Boltzmann weight measured in the ensemble with the absolute value of the Boltzmann weight, i.e. the so-called real action ensemble. The latter is particularly simple in the canonical formulation because the complex fermionic contributions are just a pure phase, where φ x ∈ {0, ±2π 3} is the phase of the Z(3) spin at site x. Indicating the average w.r.t. the absolute value of the measure by the subscript ⋅ , we have which we recognize as eq. (20). It is easy to construct an improved estimator for this quantity. Integrating out the spins we realize that the denominator becomes where P(N Q , V ), defined in the previous section, just counts the number of (unrestricted) quark configurations with N Q quarks on a lattice with V sites. The numerator on the other hand yields where N (N Q , b) from the previous section counts the number of quark configurations compatible with the constraints given by the bond configuration b. Essentially, the ratio simply calculates the average fraction of allowed to all fermion configurations, where the average on the r.h.s. is over all bond configurations at zero density. However, this ratio is not likely to be well estimated by generating bond configurations at N Q = 0, because if the system is in a different phase for N Q ≠ 0 the corresponding bond configurations are expected to be qualitatively different from those at N Q = 0. In order to avoid this well-known overlap problem, one can estimate the average sign from simulations at nonzero density. In that case one measures where we introduced φ = ∑ x n x φ x . Yet another way to calculate the average sign starts from the observation that which can be rewritten using ln In Fig. 1 we show the logarithm of the average phase ln⟨exp(iφ)⟩ ⋅ as a function of the baryon density ρ B , defined by ρ B = (N Q 3 + 1 2) L 3 , for different volumes and two different values of the coupling. The left plot is for the coupling γ = 0.5508 for which the system is in the deconfined phase for all quark numbers, while the right plot is for γ = 0.5480 for which the system is in the confined phase for all the quark numbers shown. We observe that at fixed volume the average sign goes to zero exponentially fast with the baryon density with a rate that depends on the volume.
In order to quantify further the severity of the sign problem as it scales with the volume at fixed baryon density, we define a scale parameter L 0 which parameterizes the average sign as exp(−V L 3 0 ) at fixed density. In Fig. 2 we show this scale parameter as a function of the baryon density for three different couplings corresponding to three different temperatures. The top curve corresponds to a high temperature for which the system is in the deconfined phase for all densities, while the bottom curve corresponds to a low temperature for which the system is confined for the range of density values shown. The middle curve is for a temperature inbetween for which the system undergoes a first order phase transition from the confined phase at low density into the deconfined phase at higher density. This phase transition will be discussed in more detail in the next section. Here we note that L 0 decreases linearly with the density, i.e. the severity of the sign problem gets worse exponentially with increasing density, independent of whether the system is in the confined or deconfined phase. However, the scale is distinctively different in the two phases and at fixed density the sign problem in the confined phase is worse than in the deconfined phase. The behaviour is consistent with the one observed in the grand-canonical formulation (cf. Fig. 4 in [4]), but the sign problem appears to be more severe in the canonical formulation.

Results
From previous studies of the phase diagram in the Potts model [7,4] one knows that the interesting physics happens at very low density which is why we concentrated our investigation of the sign problem in the previous section on that regime. We recall that the system undergoes a first order deconfinement phase transition along a critical line starting from the zero density transition point (e µ , γ) = (0, 0.550565(10)) down to the critical endpoint at (e µ , γ) = (0.000470(2), 0.549463(13)) where the system experiences a second order phase transition in the universality class of the 3-dimensional Ising model.
We first discuss our results for the quark chemical potential µ, given by the partition function ratio in eq.(60), as a function of the baryon density. The function essentially provides the translation between the canonical and grand-canonical formulations. In Fig. 3 we show the results for various volumes V = L 3 at γ = 0.5508 when the system is in the deconfined phase. For volumes V ≥ 25 3 we observe practically no finite volume effects and we find that the chemical potential grows monotonically with increasing density. The functional dependence can be well described by an ansatz for free fermions with an effective fermion mass and we show the corresponding fit as the dashed line in Fig. 3. Next, in Fig. 4 we show the results for the same quantity at γ = 0.5496. At this coupling, the system at zero density is in the deconfined phase, and our data clearly shows a first order phase transition into the deconfined   phase indicated by the non-monotonic behaviour of the chemical potential µ as a function of the density ρ B . The non-monotonicity is a typical signature for a first order phase transition, and in Fig. 5 we zoom in to inspect the transition in more detail. The left plot illustrates how the behaviour is better and better resolved towards the thermodynamic limit. Since the transition happens at rather low density ρ B ≲ 0.04 ⋅ 10 −3 one needs at least a lattice volume of V = 40 3 to see the onset of the non-monotonicity. The right plot shows the determination of the critical chemical potential corresponding to the first order phase transition on our largest volume V = 64 3 using the Maxwell construction. The construction guarantees that the two shaded regions limited by the curve µ(ρ B ) and µ = µ c (as illustrated in the plot) have the same area. The area can be related to the interface tension associated with a first order phase transition. From the Maxwell construction one can also access the values of the upper and lower densities ρ up,low at which the system enters or exits the coexistence phase. A systematic investigation using various interpolating functions and fit ranges yields µ c = −7.7895(41) at γ = 0.5496 which agrees very well with the value determined in [4]. For the upper and lower densities we obtain ρ up = 0.04386(53) ⋅ 10 −3 and ρ low = 0.00940(14) ⋅ 10 −3 .
In Fig. 6 we show once more the quark chemical potential µ as a function of the baryon density for various volumes, this time at γ = 0.5480 which is well below the critical endpoint. In this case, the system is in the confined phase for small baryon densities and runs smoothly into the deconfined phase at larger densities. The transition is no longer a phase transition, but rather a smooth cross-over which presents itself as a monotonic transition of the quark chemical potential from the baryonic gas form to the one describing effectively free quarks. As in the previous plots the dashed line is a phenomenological, effective description of the behaviour in terms of a free fermion gas.
Next, we investigate the quark-antiquark correlator ⟨z x z * y ⟩ N Q =0 at zero density as a function of the separation r = x−y for three couplings γ = 0.3, 0.4, 0.5 in the confined phase and for one coupling γ = 0.6 in the deconfined phase in Fig. 7. In the confined phase the correlator decays exponentially with the distance for large distances, with a characteristic scale that becomes longer for stronger coupling. This behaviour is consistent with the fact that at zero density the free energies of a single quark and a single antiquark are infinite because a single quark or antiquark can not exist due to Gauss' law. This reflects itself in the expectation values for a quark or an antiquark being zero, or equivalently the quark-antiquark correlator approaching zero at large distances. Above the phase transition point γ = 0.550565, the correlator approaches a constant at large distances signalling deconfinement. This is illustrated in Fig. 7 by the data at γ = 0.6. The finite asympotic value of the correlator at large distance indicates the spontaneous breaking of the global Z(3) symmetry and reflects the fact that the expectation values for a single quark and for a single antiquark are nonzero in the deconfined phase.
The situation is changed completely when the system is considered at nonzero baryon density. Now, the free energies for a single quark and a single antiquark are finite and the quark-quark, quark-antiquark and the antiquark-antiquark correlators are screened at large distance. This is illustrated in Fig. 8 where we show the correlators as defined in eqs. We first discuss the correlators in the confined phase (left plot). As before, the quarkantiquark pair at distance r = 0 annihilates and the potential energy should therefore be zero. From eqs.(36) and (39) it is clear that in our canonical cluster formulation this is exactly realized with the improved estimator. At sufficiently large distances the expectation value is expected to factorize into the expectation values of a single quark and antiquark. The potential energy should therefore flatten out and approach the sum of the free energies of the corresponding single particles, in this case F q + Fq obtained from ⟨z⟩⟨z * ⟩ as indicated in the figure by the grey line. We note that our results indeed match the expectation, not only qualitatively but also quantitatively.
For the quark-quark potential energy in the background of N B baryons we note that two quarks at the same position are equivalent to having a single antiquark on top of an additional baryon, hence the potential energy should be equal to the free energy of a single antiquark in the corresponding background of N B + 1 baryons, obtained from ⟨z * ⟩ ′ . At large separation on the other hand, the correlator is expected to factorize and the potential energy is given as twice the free energy of a single quark, i.e. 2F q obtained from ⟨z⟩⟨z⟩. Similarly, two antiquarks at the same position are equivalent to a single quark on top of an antibaryon, so the potential energy is equal to the free energy of a single quark in the corresponding background with N B − 1 baryons obtained from ⟨z⟩ ′ , while at large separation it approaches twice the free energy of a single antiquark, i.e. 2Fq as obtained from ⟨z * ⟩⟨z * ⟩. In Fig. 8 these energies, respectively the logarithms of the corresponding expectation values, are again given as grey lines, and they illustrate that our results qualitatively match the expectations. The small differences we observe between the lines and the correlator data are finite volume corrections.
We further note that the asymptotes of the correlators at large distance are approximately equidistant. This is due to the fact that in each step from bottom to top, a quark is replaced by an antiquark, which approximately costs the same amount of free energy. One last feature to point out is the flatness of the antiquark-antiquark correlator as compared to the other correlators. It suggests that screening two antiquarks is almost as effective as screening a single quark in the background with one less baryon, i.e. the system in the confined phase screens antiquarks much more effective than quarks. This is of course due to the fact that the system only contains static quarks but not antiquarks as dynamical degrees of freedom.
Let us finally discuss the correlators at finite density in the deconfined phase at γ = 0.6 (right plot). The discussions above for the confined phase essentially carry over without change, the only exception being that now all the correlators are very flat, i.e. they change by less than an order of magnitude over the whole range of distances. This indicates that in the deconfined phase a quark and an antiquark are screened equally efficient, i.e. the system does no longer discriminate between a quark and an antiquark.

Conclusions
We have used a cluster algorithm to solve the notorious sign problem in the Potts model approximation of heavy-dense QCD, i.e. QCD with heavy quarks and large chemical potential. Our approach is based on the fact that in the canonical formulation the position of the quarks determine the triality of the clusters, i.e. the transformation properties of the weights of the clusters under Z(3) transformations. This in turn allows to identify those configurations for which the contributions to the partition function cancel exactly after averaging each cluster over its Z(3) orientations, such that only configurations remain which yield positive contributions. In this sense, the concept of cluster triality is similar to the concept of meron clusters formulated in different contexts [8] and hence the algorithm formulated here belongs to the class of meron-custer algorithms [9,10].
The factorization of the weights into clusters with the subsequent subaveraging of each cluster immediately yields an increase of statistics by a factor 3 N C where N C is the number of clusters in a configuration. Since the number of clusters grows with the volume, at otherwise fixed parameters, the gain in statistics is exponential in the volume. We note that the mechanism at work here is similar to the one formulated in [11,12] using the subset method. Here it is made to work for a three-dimensional system and away from the strong coupling limit. Furthermore, the triality properties of the clusters also allow the construction of improved estimators which receive only positive contributions. This is of course crucial for the construction of efficient simulation algorithms.
We note that the solution of the sign problem for the Potts model presented here is not the first one. In [4] the authors constructed a cluster algorithm for the Potts model at finite density in the grand-canonical formulation which solved the sign problem completely. Another approach is to use the density of states method [13]. The sign problem can also be avoided by formulating the Potts model in the flux representation [14,4]. For this reason the focus of the present paper was not so much on the physics but rather on the mechanism which forms the basis for the solution of the sign problem. In that sense the solution here can be considered as a proof of a concept which can be applied in other cases.
Various extension of our solution of the sign problem in the canonical formulation are possible. In this paper, following the spirit of [4] we have only considered positive quark occupation numbers, corresponding to a system with only quarks but no antiquarks. In the context of QCD, a more natural scenario would be to also allow for antiquarks which can form antibaryons and mesons. The triality constraints for the contributions to the partition functions and improved estimators carry over to this case without modifications. The fermion number update described in Sec. 3 and 4 becomes even simpler, because a quark-antiquark pair can be generated within a cluster using the update step in eq. (22) where the first constraint is replaced by (1 − δ nx,−nmax ). The possibility of this step enhances the efficiency of the update considerably, especially at low densities. On the other hand, the reweighting with the multiplicities as described in Sec. 5 becomes more complicated because the calculation of the multiplicities is more complex.
One can also think about extending our approach to the gauge group SU (3). In this case one needs to embed cluster variables in continuous groups [15,16], in particular Z(3) into SU(3) [17]. In [18] a cluster algorithm has been constructed along these lines for SU(2) on a single time slice. Extending this construction to SU(3) and combining it with the approach outlined in this paper is currently under investigation.
which is valid for f < 1, we find that the effective action for the grand-canonical partition function is then ) . (75) For our simulations we have h of the order of 10 −3 so the higher order terms, O(h 2 ), contribute very little. If we keep only the first order terms our effective action simplifies to This is the same action as the one in the Potts model studied by Alford et al. in [4], hence we are able to reproduce the phase diagram at low density, including the transition point at zero density (h, γ) = (0, 0.550565 (10)) and the critical endpoint at (h, γ) = (0.000470 (2), 0.549463(13)), described in Fig. 1 of that paper. Note that in [4] the parameter γ is denoted as κ.

B Fully-dynamic connectivity problem
The Potts model in the bond formulation, as given by the partition function in eq.(13), belongs to the class of random cluster models. Simulating such systems with a local update algorithm requires one to check whether the bond under consideration is a bridge or not. A bridge bond connects two otherwise separate clusters when activated, or equivalently, breaks a cluster into two disjoint clusters when deactivated. The determination of the dynamically changing connectivity of bond clusters is a highly complex problem which in graph theory and computer science is well known under the name fully-dynamic connectivity problem. The issue here is to find data structures and algorithms with a minimal computational complexity for connectivity queries for the clusters, while dynamically adding and deleting bonds.
In the context of the Potts model [19], the problem has first been addressed by Sweeny in [20], based on the equivalence of the Potts model and the random cluster model described by Fortuin and Kasteleyn [6]. An alternative approach was suggested later by Swendsen and Wang [5] and is based on alternating updates of spin and bond variables which require only incremental connectivity information. It is well known that the Swendsen-Wang cluster construction can be dealt with efficiently with the Hoshen-Kopelman algorithm [21]. This is essentially a Find-Union algorithm acting on weighted tree data structures. Combining it with path compression to minimize the depth of the trees yields an efficient algorithm for which connectivity queries essentially have constant cost independent of the system size L.
Applying the same strategy to the situation of dynamically changing connectivity leads to computational critical slowing down, because whenever a bond is updated, one needs to check whether that bond acts as a bridge. Doing so requires to scan a number of bonds proportional to the typical size of the clusters. Close to the phase transition the clusters grow with the volume according to some critical exponent β and the computational cost hence grows like O(L β ) [22]. In contrast, a nearly-optimal solution for the fully-dynamic connectivity problem makes use of a spanning forest of Euler tour trees [23,24] in order to maintain a data structure which allows for an efficient checking of the bridge property. The runtime of the corresponding algorithm only grows as O(ln 2 L) and we make use of it in our implementation. We note that further improvements are possible [25] yielding asymptotically to a computational complexity O(ln 2 L ln ln L) for bond updates and O(ln L ln ln L) for connectivity queries. Additional improvements can be obtained by using randomized data structures [23]. However, one should consider that all these elaborate data structures involve some computational overhead and it is not clear a priori which of the strategies is most efficient. The fermion update in our algorithm for example only requires connectivity queries for which the Find-Union ansatz is most efficient. Finally we note that the efficiencies of some of the strategies discussed above have been investigated, and when possible compared to the one of the Swendsen-Wang algorithm, for the generic random cluster model in two dimensions in [22,26].

C Baryon partitions
Here we describe an algorithm which generates all partitions of N B baryons over N C clusters. The logic of the algorithm is to impose an ordering relation between different partitions of N B