text:
# Solution of the Complex Action Problem in the Potts Model for Dense QCD 11footnote 1This work is supported in part by funds provided by the U.S. Department of Energy (D.O.E.) under cooperative research agreements DE-FC02-94ER40818 and DE-FG02-96ER40945.

## 1 Introduction and Summary

Non-perturbative dense QCD can presently not be studied from first principles because Monte Carlo simulations of lattice QCD with non-zero baryon chemical potential $`\mu `$ suffer from a severe complex action problem. The Boltzmann factor in the path integral can then not be interpreted as a probability and standard importance sampling methods fail. In particular, when the $`\mu `$-dependent part of the Boltzmann factor is included in the measured observables, due to severe cancellations the required statistics is exponentially large in the space-time volume .

The complex action problem prevents the numerical simulation of a large class of interesting physical systems including other field theories at non-zero chemical potential or non-zero $`\theta `$-vacuum angle as well as some fermionic field theories with an odd number of flavors. A special case of the complex action problem is the so-called fermion sign problem which arises for fermionic path integrals formulated in a Fock state basis. The problem is due to paths that correspond to an odd permutation of fermion positions which contribute negatively to the path integral. There are numerous condensed matter systems ranging from the repulsive Hubbard model away from half-filling to antiferromagnetic quantum spin systems in an external magnetic field that cannot be simulated with standard Monte Carlo algorithms. Meron-cluster algorithms have been used to solve the sign or complex action problems in several of these cases. For example, the first meron-cluster algorithm has led to a solution of the complex action problem in the 2-d $`O(3)`$ symmetric field theory at non-zero $`\theta `$-vacuum angle . In this model, some of the clusters are half-instantons, so they are called meron-clusters. The complex action problem also arises in the 2-d $`O(3)`$ model at non-zero chemical potential. When formulated as a D-theory  — i.e. in terms of discrete variables that undergo dimensional reduction — the complex action problem has also been solved with a meron-cluster algorithm . Recently, the meron concept has been generalized to fermions . Meron-cluster algorithms have led to a complete solution of the fermion sign problem in a variety of models including non-relativistic spinless fermions , relativistic staggered fermions  and some models in the Hubbard model family .222The models investigated so far only show s-wave superconductivity. Recently, a meron-cluster algorithm has been used to solve the sign problem that arises for quantum antiferromagnets in an external magnetic field . For a review of these recent developments and a preliminary version of the present results see .

In the conventional formulation of lattice QCD the quarks are represented by Grassmann fields. When the quarks are integrated out, they leave behind a fermion determinant that acts as a non-local effective action for the gluons. At zero chemical potential and for an even number of flavors, the fermion determinant is real and positive and can thus be interpreted as a probability for generating gluon field configurations. Despite the fact that standard importance sampling techniques apply, the non-local nature of the effective gluon action makes lattice QCD simulations with dynamical fermions very time consuming. With a non-zero chemical potential for the baryon number, the fermion determinant becomes complex and standard importance sampling techniques fail completely . This is the reason why non-perturbative QCD at non-zero baryon density can presently not be studied from first principles.

It is natural to ask if a meron-cluster algorithm could be used to solve the complex action problem in QCD. When one integrates out light quarks, one obtains a non-local effective action for the gluons and it appears unlikely that the meron concept will apply. On the other hand, when one describes the quarks in a Fock state basis, the complex action problem is still present, in the form of a fermion sign problem. Our hope is that this problem will eventually be solved by a meron-cluster algorithm applied to the D-theory formulation of QCD , since the quark and gluon degrees of freedom are then discrete and are much easier to handle. In this paper we address a simpler situation first. We consider QCD in the limit of very heavy quarks with a large chemical potential. These can be integrated out, introducing Polyakov loops into the effective gluon action. When quarks are integrated out at non-zero chemical potential $`\mu `$ we expect a complex action, and in this case it arises because a Polyakov loop $`\mathrm{\Phi }`$ and its charge conjugate $`\mathrm{\Phi }^{}`$ get different weights when $`\mu 0`$.

Polyakov loops are only non-local in the Euclidean time direction, so this effective gluon action is more manageable than the one that arises for a general fermion determinant. Indeed, Blum, Hetrick and Toussaint have simulated the theory in this form on lattices of moderate size where the complex action problem is less severe . Recently, Engels, Kaczmarek, Karsch and Laermann have studied QCD with heavy quarks at fixed baryon number. Again, for moderate baryon density and moderate volumes the complex action problem is not too severe and simulations are possible . Ultimately one would like to be able to solve the complex action problem for this gluon action completely. At the moment, we still cannot apply a meron-cluster algorithm to solve the problem, because the construction of efficient cluster algorithms for non-Abelian gauge theories seems to be impossible for Wilson’s formulation of lattice field theory. Here we will simplify the problem further by replacing the gauge dynamics by that of the $`(3)`$ Potts model representing the Polyakov loops . We have found a cluster algorithm that solves this complex action problem in the Potts model approximation to QCD.

The 3-d $`(3)`$-symmetric Potts model has often been used as an approximation to QCD with static quarks. In particular, the phase transition to a broken $`(3)`$ symmetry phase at high temperature corresponds to the first order deconfining phase transition in QCD. As has been noted by Condella and DeTar, a term that corresponds to a chemical potential can also be included in the Potts model, explicitly breaking the $`(3)`$ symmetry . As the coefficient of this term grows, the first order deconfinement phase transition persists but it becomes weaker and ultimately disappears in a critical endpoint. This point is expected to be in the universality class of the 3-d Ising model. In this paper we will confirm this expectation with numerical simulations.

In principle, one can imagine deriving an effective 3-d 3-state Potts model directly from QCD by integrating out all degrees of freedom except for the $`(3)`$ phase of the Polyakov loop. However, the resulting Potts model action would be very complicated and cannot be derived in practice, except in the strong coupling limit. Here we approximate QCD with heavy quarks by a 3-d $`(3)`$-symmetric Potts model with a standard nearest-neighbor interaction. Universal features like the nature of the critical endpoint of the deconfinement phase transition are correctly reproduced in this approximation. Figure 1 contains the phase diagram of the

3-d 3-state Potts model in the $`(h,\kappa )`$-plane. The parameter $`h`$ represents $`\mathrm{exp}(\beta (\mu M))`$ in QCD with quarks of mass $`M`$ at chemical potential $`\mu `$. We study the limit $`M,\mu \mathrm{}`$ for any given $`\mu M`$. Large $`h`$ corresponds to $`\mu >M`$ and small $`h`$ to $`\mu <M`$. Because $`\mu MM,\mu `$ we are always, for any $`h`$, in the immediate neighborhood of the onset of non-zero density for the heavy quarks. This means that it does not matter whether they are fermions or bosons, since they never move. The difference would only become apparent above the onset, where either a Fermi surface or a degenerate Bose gas would occur, and our order of limits is such that we never get that far from the onset. The parameter $`\kappa `$ is the standard Potts model coupling, which corresponds roughly to the temperature $`T=1/\beta `$. The ordinary first-order deconfinement phase transition at $`h=0`$ (point $`T`$ in Figure 1) extends into a line of first order transitions that terminates in the critical endpoint $`E`$. This endpoint occurs at such a low value of $`h`$ that the complex action problem is not very severe there, and we found that the most efficient way to locate and study it was to employ a reweighted Metropolis algorithm, which can in this case be applied at volumes large enough to show the critical behavior. Similar methods were used recently by Karsch and Stickan  in a version of the 3-d 3-state Potts model where the action is real, and the endpoint was found to have the critical exponents of the 3-d Ising model. We find that in the Potts model with complex action the endpoint has the same critical properties. Furthermore its position is barely shifted in comparison to the model with real action. In this paper we do not limit our attention to the endpoint, but develop a method that solves the complex action problem everywhere in the phase diagram.

We also calculate the potentials between static quarks and anti-quarks in the Potts model approximation to QCD. In the confined phase at $`\mu =0`$ the static quark-anti-quark potential is linearly rising with the distance as a manifestation of confinement. For the same reason the quark-quark and anti-quark-anti-quark potentials are infinite at all distances. In the deconfined phase the quark-anti-quark potential reaches a plateau at twice the (now finite) free energy of a quark. Similarly, the quark-quark and anti-quark-anti-quark potentials are no longer infinite. It should be noted that quark-quark and anti-quark-anti-quark potentials are usually not calculated in lattice simulations. This is because—as a consequence of the $`(3)`$ Gauss law—quark or anti-quark pairs cannot exist in a finite spatial volume with periodic boundary conditions . Interestingly, this changes for $`\mu 0`$ because then there are compensating background charges in the medium that can absorb the $`(3)`$ flux of an external quark. Since the chemical potential explicitly breaks the $`(3)`$ symmetry, there is no longer a clear distinction between confinement and deconfinement for $`\mu 0`$. This manifests itself in the phase diagram by the fact that confined and deconfined phases are analytically connected.

Figure 2 shows the quark-anti-quark, quark-quark and anti-quark-anti-quark potentials on the confined side (a) and on the deconfined side (b) of the crossover. Note that at $`\mu 0`$ even in the confined phase the quark-anti-quark potential now reaches a plateau. The plateau height corresponds to the sum of the free energies $`F_Q`$ of an external quark and $`F_{\overline{Q}}`$ of an external anti-quark. For $`\mu >0`$ quarks are favored in the medium while anti-quarks are suppressed. As a consequence, the free energy of an external static quark is larger than that of an external static anti-quark. While an external static anti-quark can bind with a single background quark from the medium and form a meson, an external static quark needs two quarks from the medium to form a baryon. Indeed, on the confined side of the transition $`F_Q`$ is clearly larger than $`F_{\overline{Q}}`$, while on the deconfined side $`F_Q`$ and $`F_{\overline{Q}}`$ are more or less the same. We have normalized the potentials such that at zero distance a static quark-anti-quark pair has zero energy. In the Potts model, two quarks at zero distance are indistinguishable from a single anti-quark, and similarly, two anti-quarks on top of each other behave like a single quark. Hence, at zero distance the quark-quark potential $`V_{QQ}(0)`$ agrees with the free energy of a single anti-quark $`F_{\overline{Q}}`$ and the anti-quark-anti-quark potential obeys $`V_{\overline{Q}\overline{Q}}(0)=F_Q`$. At asymptotic distances the potentials $`V_{Q\overline{Q}}(\mathrm{})`$, $`V_{QQ}(\mathrm{})`$ and $`V_{\overline{Q}\overline{Q}}(\mathrm{})`$ take the values $`F_Q+F_{\overline{Q}}`$, $`2F_Q`$ and $`2F_{\overline{Q}}`$, respectively. This behavior is consistent with our numerical data shown in fig.2.

In the absence of a chemical potential, the Potts model can be simulated with the original Swendsen-Wang cluster algorithm . When a chemical potential is introduced, the Potts model suffers from the complex action problem and standard importance sampling methods including the cluster algorithm fail. In this paper, we construct an improved estimator for the $`\mu `$-dependent part of the Boltzmann factor by averaging analytically over all configurations related to each other by cluster flips. In contrast to the original Boltzmann factor, the improved estimator is real and positive and can be used for importance sampling. This solves the complex action problem completely.

Although the Potts model inherits the complex action problem from QCD, it can be transformed into a “flux model” that has no complex action problem . The flux model has been simulated in  and the disappearance of the first order deconfining phase transition at large chemical potential has been observed numerically. These results may give encouragement to the hope that QCD itself could be transformed into a model without a complex action problem. In this paper we show that, at least for the Potts model, it is more efficient to leave it in its usual form and solve the complex action problem with our cluster algorithm than to transform it into a flux model and use conventional Metropolis methods.

The paper is organized as follows. Section 2 contains a derivation of the effective gluon action resulting from static quarks with large chemical potential as well as its Potts model approximation. In section 3 we describe the cluster algorithm that solves the complex action problem. Section 4 contains the derivation of the flux representation of the Potts model and a description of a Metropolis algorithm to simulate it. A comparison of the Metropolis algorithm for the flux model and the cluster algorithm for the original Potts model shows that the latter is more efficient. In section 5 we present the physical results concerning the critical endpoint $`E`$. Using finite-size scaling techniques, we are able to determine the position of the critical endpoint of the deconfinement phase transition to high accuracy. Our results are consistent with the expected universal 3-d Ising behavior. Finally, section 6 contains our conclusions.

## 2 QCD with Heavy Quarks and the 3-d 3-State Potts Model

The partition function for a pure $`SU(3)`$ Yang-Mills theory is given by

$$Z=𝒟A\mathrm{exp}(S[A]),$$
(2.1)
where

$$S[A]=_0^\beta 𝑑td^3x\frac{1}{2g^2}\text{Tr}[F_{\mu \nu }F_{\mu \nu }],$$
(2.2)
is the Euclidean action for the gluons and $`\beta `$ is the inverse temperature. The action is invariant under gauge transformations

$$g(\stackrel{}{x},0)=g(\stackrel{}{x},\beta )z,$$
(2.3)
that are periodic in Euclidean time up to an element $`z`$ of the center $`(3)=\{\mathrm{exp}(2\pi in/3),n=1,2,3\}`$ of the non-Abelian gauge group. In the presence of a single external heavy quark of bare mass $`M`$ at an undetermined position $`\stackrel{}{x}`$ the partition function turns into

$$Z_Q=𝒟A\mathrm{\Phi }[A]\mathrm{exp}(S[A])\mathrm{exp}(\beta M),$$
(2.4)
where

$$\mathrm{\Phi }[A]=d^3x\text{Tr}[𝒫\mathrm{exp}(_0^\beta 𝑑tA_4(\stackrel{}{x},t))],$$
(2.5)
is the spatial integral of the Polyakov loop. Ultimately, the mass $`M`$ will be sent to infinity. Note that while the center transformation of eq.(2.3) leaves the pure gluon action $`S[^gA]=S[A]`$ invariant, the Polyakov loop transforms into

$$\mathrm{\Phi }[^gA]=z\mathrm{\Phi }[A].$$
(2.6)
This shows that in the presence of the external quark, the $`(3)`$ symmetry is explicitly broken. The partition function for a system of gluons in the presence of a single heavy anti-quark is given by

$$Z_{\overline{Q}}=𝒟A\mathrm{\Phi }[A]^{}\mathrm{exp}(S[A])\mathrm{exp}(\beta M),$$
(2.7)
where $``$ denotes complex conjugation. Let us now consider a system of gluons in a background of $`n`$ static quarks and $`\overline{n}`$ static anti-quarks. The partition function then takes the form

$$Z_{n,\overline{n}}=𝒟A\frac{1}{n!}\mathrm{\Phi }[A]^n\frac{1}{\overline{n}!}(\mathrm{\Phi }[A]^{})^{\overline{n}}\mathrm{exp}(S[A])\mathrm{exp}(\beta M(n+\overline{n})).$$
(2.8)
The factors $`1/n!`$ and $`1/\overline{n}!`$ appear because quarks are indistinguishable, as are antiquarks. Introducing the quark chemical potential $`\mu `$ that couples to $`(n\overline{n})`$, i.e. three times the baryon number, we obtain the grand canonical partition function

$`Z(\mu )`$ $`=`$ $`{\displaystyle \underset{n,\overline{n}}{}}Z_{n,\overline{n}}\mathrm{exp}(\beta \mu (n\overline{n}))`$
$`=`$ $`{\displaystyle \underset{n,\overline{n}}{}}{\displaystyle 𝒟A\frac{1}{n!}\mathrm{\Phi }[A]^n\frac{1}{\overline{n}!}(\mathrm{\Phi }[A]^{})^{\overline{n}}\mathrm{exp}(S[A]\beta n(M\mu )\beta \overline{n}(M+\mu ))}`$
$`=`$ $`{\displaystyle 𝒟A\mathrm{exp}(S[A]+\mathrm{exp}(\beta (M\mu ))\mathrm{\Phi }[A]+\mathrm{exp}(\beta (M+\mu ))\mathrm{\Phi }[A]^{})}.`$
(2.9)
As expected, the presence of quarks and anti-quarks leads to an explicit breaking of the $`(3)`$ center symmetry. Furthermore, in the presence of a non-zero chemical potential the effective action for the gluons is complex. Note that in the $`SU(2)`$ case the action remains real because then the Polyakov loop itself is real, i.e. $`\mathrm{\Phi }[A]^{}=\mathrm{\Phi }[A]`$. The action becomes real even in the $`SU(3)`$ case if $`\mu `$ is purely imaginary. Furthermore, one can see that the chemical potential explicitly breaks the charge conjugation symmetry that replaces $`\mathrm{\Phi }[A]`$ by $`\mathrm{\Phi }[A]^{}`$. In fact, under charge conjugation the action turns into its complex conjugate. We have assumed that the quarks are static. Hence, to be consistent we must consider the limit $`M\mathrm{}`$. In order to obtain a non-trivial result, we simultaneously take the limit $`\mu \mathrm{}`$ such that $`M\mu `$ remains finite. The partition function then simplifies to

$$Z(\mu )=𝒟A\mathrm{exp}(S[A]+\mathrm{exp}(\beta (M\mu ))\mathrm{\Phi }[A]).$$
(2.10)
As discussed in  and , a similar result can be obtained by simplifying the full QCD quark determinant in the static quark limit. In general the determinant would contain all Wilson loops, but because $`M`$ is large most of them are suppressed. The only ones that survive are those for which the enhancement due to the chemical potential compensates for the suppression due to the mass. These are the Polyakov loops that progress in a straight line from Euclidean time $`t=0`$ to $`t=\beta `$ at some position $`\stackrel{}{x}`$. In the loop expansion of the quark determinant, each of these has a weight $`\mathrm{exp}(\beta (M\mu ))`$.

Up to this point we have treated QCD consistently in the static quark limit. The resulting effective action for the gluons is complex and we presently don’t know how to simulate it efficiently. For that reason we now replace the gluon system by a simple 3-d lattice 3-state Potts model. The Potts spins $`\mathrm{\Phi }_x(3)`$ replace the original Polyakov loop variables and the partition function turns into

$$Z(h)=𝒟\mathrm{\Phi }\mathrm{exp}(S[\mathrm{\Phi }]+h\underset{x}{}\mathrm{\Phi }_x),$$
(2.11)
where $`h`$ replaces $`\mathrm{exp}(\beta (M\mu ))`$. Note that the Potts model action is still complex. In principle, one can imagine integrating out all QCD degrees of freedom except for the $`(3)`$ phase of the Polyakov loop and thus derive an effective Potts model action directly from QCD. In practice this is impossible, except in the strong coupling limit. For simplicity, we therefore replace the pure gluon action $`S[A]`$ by a standard nearest-neighbor Potts model interaction

$$S[\mathrm{\Phi }]=\kappa \underset{x,i}{}\delta _{\mathrm{\Phi }_x,\mathrm{\Phi }_{x+\widehat{i}}}.$$
(2.12)
The coupling constant $`\kappa `$ is not related in a simple way to the parameters of QCD. Still, a large value of $`\kappa `$ corresponds qualitatively to the high-temperature deconfined phase, while small $`\kappa `$ values correspond to the confined phase. As mentioned in the introduction, the Potts model also retains the general features of the QCD phase diagram. At $`h=0`$ ($`M`$ infinite, $`\mu `$ finite) there is a first-order phase transition as a function of $`\kappa `$, between the disordered (confined) phase that respects the $`(3)`$ symmetry and the ordered (deconfined) phase that spontaneously breaks it. An order parameter for this transition is $`\mathrm{\Phi }`$. As $`h`$ rises from zero, the chemical potential term explicitly breaks the $`(3)`$ symmetry, the phase transition weakens, and then ends at a critical point. Correspondingly, in heavy-quark QCD the quarks begin to contribute to the partition function when $`\mu `$ gets close to $`M`$, and there is no longer an order parameter for deconfinement. The deconfining phase transition terminates at a critical endpoint.

## 3 Cluster Algorithm Solution of the Complex Action Problem

In this section we first discuss the general nature of the complex action problem and then discuss the cluster algorithm that solves this problem for the Potts model. We also construct improved estimators for various physical quantities.

### 3.1 The General Nature of the Complex Action Problem

When the action is complex the resulting Boltzmann factor cannot be interpreted as a probability and hence standard importance sampling techniques fail. When one uses just the absolute value of the Boltzmann factor for importance sampling and includes its complex phase in measured observables $`O`$, expectation values take the form

$`O`$ $`=`$ $`{\displaystyle \frac{1}{Z}}{\displaystyle 𝒟\mathrm{\Phi }O[\mathrm{\Phi }]\mathrm{exp}(S[\mathrm{\Phi }]+h\underset{x}{}\mathrm{\Phi }_x)}`$ (3.1)
$`=`$ $`{\displaystyle \frac{O\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R}{\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R}}.`$
The subscript $`R`$ refers to a modified ensemble with a real action described by the partition function

$$Z_R=𝒟\mathrm{\Phi }\mathrm{exp}(S[\mathrm{\Phi }]+h\underset{x}{}\text{Re}\mathrm{\Phi }).$$
(3.2)
By definition we have

$`\mathrm{exp}(ih{\displaystyle \underset{x}{}}\text{Im}\mathrm{\Phi }_x)_R`$ $`=`$ $`{\displaystyle \frac{1}{Z_R}}{\displaystyle 𝒟\mathrm{\Phi }\mathrm{exp}(ih\underset{x}{}\text{Im}\mathrm{\Phi }_x)\mathrm{exp}(S[\mathrm{\Phi }]+h\underset{x}{}\text{Re}\mathrm{\Phi })}`$ (3.3)
$`=`$ $`{\displaystyle \frac{Z}{Z_R}}\mathrm{exp}(V(ff_R)),`$
where $`f`$ and $`f_R`$ are the free energy densities of the original complex and the modified real action systems, respectively, and $`V`$ is the spatial volume. Hence, the denominator in eq.(3.1) becomes exponentially small as one increases the volume. The same is true for the numerator, because $`O`$ itself is not exponentially large in $`V`$.

Although, in principle, simulating the modified ensemble is correct, in practice this method fails for large volumes. The reason is that observables are obtained as ratios of exponentially small numerators and denominators which are themselves averages of quantities of order one. This leads to very severe cancellations and requires an exponentially large statistics in order to obtain accurate results. To see this, we estimate the relative statistical error in the determination of the average phase of the Boltzmann factor $`\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)`$. Since $`\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R=Z/Z_R`$ the average itself is real. When one generates $`N`$ statistically independent field configurations in a Monte Carlo simulation, the resulting error to signal ratio is given by

$`{\displaystyle \frac{\mathrm{\Delta }\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)}{\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R}}`$ $`=`$ $`{\displaystyle \frac{\sqrt{|\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R|^2_R}}{\sqrt{N}\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R}}`$ (3.4)
$`=`$ $`{\displaystyle \frac{\sqrt{1\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R^2}}{\sqrt{N}\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R}}{\displaystyle \frac{\mathrm{exp}(V(ff_R))}{\sqrt{N}}}.`$
For large $`V`$ we have used $`\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R1`$ as implied by eq.(3.3). Consequently, in order to obtain an acceptable error to signal ratio one must generate at least $`N\mathrm{exp}(2V(ff_R))`$ configurations. For large volumes this is impossible in practice.

### 3.2 The Cluster Algorithm for the Potts Model

Let us now outline the ideas that underlie the cluster algorithm that we use to solve the complex action problem. It is based on the original Swendsen-Wang cluster algorithm  for the Potts model without chemical potential. In fact, in the limit $`h=0`$ our algorithm reduces to that algorithm. The Swendsen-Wang cluster algorithm decomposes the lattice into independent clusters of connected sites. Each spin belongs to exactly one cluster and all spins within a cluster are assigned the same random $`(3)`$ element. In this paper, we construct an improved estimator for the $`h`$-dependent part $`\mathrm{exp}(h_x\mathrm{\Phi }_x)`$ of the Boltzmann factor by analytically averaging it over all configurations related to each other by cluster flips. Although, for an individual configuration $`\mathrm{exp}(h_x\mathrm{\Phi }_x)`$ is in general complex, its improved estimator is always real and positive and can thus be used for importance sampling. This completely solves the complex action problem.

Let us first describe the original Swendsen-Wang algorithm for $`h=0`$. In this method one introduces variables $`b=0,1`$ for each bond connecting neighboring lattice sites $`x`$ and $`y=x+\widehat{i}`$ and one writes the nearest neighbor Boltzmann factor as

$$\mathrm{exp}(\kappa \delta _{\mathrm{\Phi }_x,\mathrm{\Phi }_y})=\underset{b=0,1}{}[\delta _{b,1}\delta _{\mathrm{\Phi }_x,\mathrm{\Phi }_y}(e^\kappa 1)+\delta _{b,0}].$$
(3.5)
In the enlarged configuration space of spin and bond variables, the bond variables impose constraints between the spin variables. When a bond is put (i.e. when $`b=1`$), the spin Boltzmann factor is $`\delta _{\mathrm{\Phi }_x,\mathrm{\Phi }_y}(e^\kappa 1)`$ and hence the spin variables $`\mathrm{\Phi }_x`$ and $`\mathrm{\Phi }_y`$ at the two ends of the bond must be identical. On the other hand, when the bond is not put ($`b=0`$), the spin Boltzmann factor is $`1`$ and thus the variables $`\mathrm{\Phi }_x`$ and $`\mathrm{\Phi }_y`$ fluctuate independently. The spin variables, in turn, determine the probability to put a bond. When the spins $`\mathrm{\Phi }_x`$ and $`\mathrm{\Phi }_y`$ are different, the bond Boltzmann factor is $`\delta _{b,0}`$ and thus the bond is not put. On the other hand, when $`\mathrm{\Phi }_x`$ and $`\mathrm{\Phi }_y`$ are the same, the bond Boltzmann factor is $`[\delta _{b,1}(e^\kappa 1)+\delta _{b,0}]`$. Consequently, a bond between parallel spins is put with probability $`p=1e^\kappa `$. Note that for $`\kappa =0`$ no bonds are put, while for $`\kappa =\mathrm{}`$ parallel spins are always connected by a bond.

The Swendsen-Wang cluster algorithm updates bond and spin variables in alternating order. First, for a given spin configuration, bonds are put with probability $`p`$ between parallel neighboring spins. No bonds are put between non-parallel spins. Then the spins are updated according to the constraints represented by the resulting bond configuration. Spins connected by bonds must remain parallel, while spins not connected by bonds fluctuate independently. Hence, to update the spins, one must identify clusters, i.e. sets of spins that are connected by bonds. All spins in a cluster are parallel and are assigned the same random $`(3)`$ element in the spin update. All spins belong to exactly one cluster. It should be noted that a cluster may consist of a single spin. A configuration consisting of $`N_C`$ clusters can be viewed as a member of a sub-ensemble of $`3^{N_C}`$ equally probable configurations which result by assigning $`(3)`$ elements to the various clusters in all possible ways. As was already pointed out by Swendsen and Wang, one can construct improved estimators for various physical quantities by averaging analytically over all $`3^{N_C}`$ configurations in a sub-ensemble. Since the number of clusters is proportional to the volume, this effectively increases the statistics by a factor that is exponentially large in $`V`$.

Let us construct an improved estimator for the $`h`$-dependent part $`\mathrm{exp}(h_x\mathrm{\Phi }_x)`$ of the Boltzmann factor. Although for an individual configuration this term is in general complex, its average over a sub-ensemble of $`3^{N_C}`$ configurations is always real and positive. This results from the following observations. The $`h`$-dependent part of the Boltzmann factor is a product of cluster contributions

$$\mathrm{exp}(h\underset{x}{}\mathrm{\Phi }_x)=\underset{C}{}\mathrm{exp}(h\underset{xC}{}\mathrm{\Phi }_x).$$
(3.6)
Since the clusters are independent, the sub-ensemble average is a product

$$\mathrm{exp}(h\underset{x}{}\mathrm{\Phi }_x)_{3^{N_C}}=\underset{C}{}\mathrm{exp}(h\underset{xC}{}\mathrm{\Phi }_x)_3,$$
(3.7)
of 3-state averages for the individual clusters

$`\mathrm{exp}(h{\displaystyle \underset{xC}{}}\mathrm{\Phi }_x)_3`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{exp}(h|C|\mathrm{\Phi })`$ (3.8)
$`=`$ $`{\displaystyle \frac{1}{3}}[\mathrm{exp}(h|C|)+2\mathrm{exp}(h|C|/2)\mathrm{cos}(\sqrt{3}h|C|/2)]`$
$`=`$ $`W(C),`$
which defines a weight $`W(C)`$ for each cluster. We have used the fact that all spins $`\mathrm{\Phi }_x`$ in a given cluster $`C`$ take the same value $`\mathrm{\Phi }(3)`$ so that $`_{xC}\mathrm{\Phi }_x=|C|\mathrm{\Phi }`$ where $`|C|=_{xC}1`$ is the cluster size. It is easy to show that the expression in eq.(3.8) is always positive and can hence be used for importance sampling. This is crucial for a complete solution of the complex action problem.

For a given bond configuration one can integrate out the spin variables and one obtains

$`Z`$ $`=`$ $`{\displaystyle 𝒟b(e^\kappa 1)^{N_b}3^{N_C}\underset{C}{}W(C)}`$ (3.9)
$`=`$ $`{\displaystyle 𝒟b(e^\kappa 1)^{N_b}\underset{C}{}[\mathrm{exp}(h|C|)+2\mathrm{exp}(h|C|/2)\mathrm{cos}(\sqrt{3}h|C|/2)]}.`$
Here $`N_b`$ is the number of bonds that are put (i.e. have $`b=1`$). The factor $`3^{N_C}`$ represents the number of allowed spin configurations for a given bond configuration and the factors $`W(C)`$ come from the improved estimator. The effective action for the bond variables depends only on the sizes $`|C|`$ of the clusters corresponding to a given bond configuration. Note that the factor $`1/3`$ per cluster in eq.(3.8) cancels against the factor $`3^{N_C}`$.

Our algorithm directly updates the partition function of eq.(3.9), i.e. it only operates on the bond variables while the spins are already integrated out analytically.333B. Scarlet was first to realize that the spin variables need not even be simulated. The bond variables that define the clusters are updated with a local algorithm. A bond whose value does not affect the cluster sizes is put with probability $`p=1e^\kappa `$. This happens when the two sites at its ends belong to the same cluster because they are connected indirectly through other bonds. A bond whose value affects the cluster sizes is put with a probability that depends on the sizes of the corresponding clusters. When the bond is not put ($`b=0`$), its endpoints $`x`$ and $`y`$ belong to two different clusters $`C_1`$ and $`C_2`$ of sizes $`|C_1|`$ and $`|C_2|`$ and the corresponding Boltzmann weight is $`3^2W(C_1)W(C_2)`$. On the other hand, when the bond is put ($`b=1`$), its endpoints belong to the combined cluster $`C_1C_2`$ of size $`|C_1|+|C_2|`$. In that case, the Boltzmann weight is $`3W(C_1C_2)(e^\kappa 1)`$. Hence, the bond is put with probability

$$q=\frac{W(C_1C_2)(e^\kappa 1)}{3W(C_1)W(C_2)+W(C_1C_2)(e^\kappa 1)}.$$
(3.10)
### 3.3 Improved Estimators for Physical Quantities

In order to measure physical observables, it is crucial to construct improved estimators for them as well. Here we construct improved estimators for the Polyakov loop $`\mathrm{\Phi }_x`$, its charge conjugate $`\mathrm{\Phi }_x^{}`$, as well as for the correlators $`\mathrm{\Phi }_x\mathrm{\Phi }_y^{}`$, $`\mathrm{\Phi }_x\mathrm{\Phi }_y`$ and $`\mathrm{\Phi }_x^{}\mathrm{\Phi }_y^{}`$. The expectation values

$$\mathrm{\Phi }_x=\mathrm{exp}(\beta F_Q),\mathrm{\Phi }_x^{}=\mathrm{exp}(\beta F_{\overline{Q}}),$$
(3.11)
determine the free energies $`F_Q`$ of a quark and $`F_{\overline{Q}}`$ of an anti-quark. The Polyakov loop correlators determine the quark-anti-quark potential $`V_{Q\overline{Q}}(xy)`$, the quark-quark potential $`V_{QQ}(xy)`$ and the anti-quark-anti-quark potential $`V_{\overline{Q}\overline{Q}}(xy)`$ via

$`\mathrm{exp}(\beta V_{Q\overline{Q}}(xy))=\mathrm{\Phi }_x\mathrm{\Phi }_y^{},`$
$`\mathrm{exp}(\beta V_{QQ}(xy))=\mathrm{\Phi }_x\mathrm{\Phi }_y,`$
$`\mathrm{exp}(\beta V_{\overline{Q}\overline{Q}}(xy))=\mathrm{\Phi }_x^{}\mathrm{\Phi }_y^{}.`$ (3.12)
The improved estimator for the Polyakov loop is given by the sub-ensemble average

$$\mathrm{\Phi }_x\mathrm{exp}(h\underset{z}{}\mathrm{\Phi }_z)_{3^{N_C}}=\frac{1}{3}\underset{\mathrm{\Phi }(3)}{}\mathrm{\Phi }\mathrm{exp}(h|C_x|\mathrm{\Phi })\underset{CC_x}{}W(C),$$
(3.13)
where $`C_x`$ is the cluster that contains the point $`x`$. Hence, we obtain

$$\mathrm{\Phi }_x=\frac{1}{Z(h)}𝒟b\frac{1}{3W(C_x)}\underset{\mathrm{\Phi }(3)}{}\mathrm{\Phi }\mathrm{exp}(h|C_x|\mathrm{\Phi })(e^\kappa 1)^{N_b}3^{N_C}\underset{C}{}W(C),$$
(3.14)
i.e. after integrating out the spin variables, the Polyakov loop is represented by

$$\mathrm{\Phi }_x=\frac{1}{3W(C_x)}\underset{\mathrm{\Phi }(3)}{}\mathrm{\Phi }\mathrm{exp}(h|C_x|\mathrm{\Phi }).$$
(3.15)
Similarly, the operator representing the charge conjugate Polyakov loop is

$$\mathrm{\Phi }_x^{}=\frac{1}{3W(C_x)}\underset{\mathrm{\Phi }(3)}{}\mathrm{\Phi }^{}\mathrm{exp}(h|C_x|\mathrm{\Phi }).$$
(3.16)
The improved estimator for the Polyakov loop correlator $`\mathrm{\Phi }_x\mathrm{\Phi }_y^{}`$ is given by the sub-ensemble average

$`\mathrm{\Phi }_x\mathrm{\Phi }_y^{}\mathrm{exp}(h{\displaystyle \underset{z}{}}\mathrm{\Phi }_z)_{3^{N_C}}`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }\mathrm{exp}(h|C_x|\mathrm{\Phi }){\displaystyle \frac{1}{3}}{\displaystyle \underset{\mathrm{\Phi }^{}(3)}{}}\mathrm{\Phi }^{}\mathrm{exp}(h|C_y|\mathrm{\Phi })`$ (3.17)
$`\times `$ $`{\displaystyle \underset{CC_x,C_y}{}}W(C),`$
if the points $`x`$ and $`y`$ belong to two different clusters $`C_x`$ and $`C_y`$. If the points $`x`$ and $`y`$ belong to the same cluster (i.e. if $`C_x=C_y`$), the improved estimator is simply given by

$$\mathrm{\Phi }_x\mathrm{\Phi }_y^{}\mathrm{exp}(h\underset{z}{}\mathrm{\Phi }_z)_{3^{N_C}}=\underset{C}{}W(C),$$
(3.18)
because then $`\mathrm{\Phi }_x\mathrm{\Phi }_y^{}=1`$. Hence, the operator representing the correlator is

$`\mathrm{\Phi }_x\mathrm{\Phi }_y^{}`$ $`=`$ $`{\displaystyle \frac{1}{9W(C_x)W(C_y)}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }\mathrm{exp}(h|C_x|\mathrm{\Phi }){\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }^{}\mathrm{exp}(h|C_y|\mathrm{\Phi }),\text{if}C_xC_y,`$
$`\mathrm{\Phi }_x\mathrm{\Phi }_y^{}`$ $`=`$ $`1,\text{if}C_x=C_y.`$ (3.19)
Similarly, we have

$`\mathrm{\Phi }_x\mathrm{\Phi }_y`$ $`=`$ $`{\displaystyle \frac{1}{9W(C_x)W(C_y)}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }\mathrm{exp}(h|C_x|\mathrm{\Phi }){\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }\mathrm{exp}(h|C_y|\mathrm{\Phi }),\text{if}C_xC_y,`$
$`\mathrm{\Phi }_x\mathrm{\Phi }_y`$ $`=`$ $`{\displaystyle \frac{1}{3W(C_x)}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }^2\mathrm{exp}(h|C_x|\mathrm{\Phi }),\text{if}C_x=C_y,`$
$`\mathrm{\Phi }_x^{}\mathrm{\Phi }_y^{}`$ $`=`$ $`{\displaystyle \frac{1}{9W(C_x)W(C_y)}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }^{}\mathrm{exp}(h|C_x|\mathrm{\Phi }){\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }^{}\mathrm{exp}(h|C_y|\mathrm{\Phi }),\text{if}C_xC_y,`$
$`\mathrm{\Phi }_x^{}\mathrm{\Phi }_y^{}`$ $`=`$ $`{\displaystyle \frac{1}{3W(C_x)}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}(\mathrm{\Phi }^{})^2\mathrm{exp}(h|C_x|\mathrm{\Phi }),\text{if}C_x=C_y.`$ (3.20)
### 3.4 Severity of the Complex Action Problem

In order to estimate the severity of the complex action problem, we like to determine the expectation value of the complex phase of the Boltzmann factor in the modified real action ensemble

$$\mathrm{exp}(ih\underset{x}{}\text{Im}\mathrm{\Phi }_x)_R=\frac{Z}{Z_R}.$$
(3.21)
Rather than implementing this directly in a simulation that uses the absolute value of the Boltzmann factor for importance sampling, one can measure $`Z_R/Z`$ with the cluster algorithm. In fact, an improved estimator for this quantity is given by $`_CW_R(C)/W(C)`$, where

$$W_R(C)=\mathrm{exp}(h\underset{xC}{}\text{Re}\mathrm{\Phi }_x)_3=\frac{1}{3}[\mathrm{exp}(h|C|)+2\mathrm{exp}(h|C|/2)],$$
(3.22)
is the weight that replaces $`W(C)`$ in the real action ensemble. Alternatively, one can construct a cluster algorithm that simulates the real action ensemble. In that case, one needs to measure $`_CW(C)/W_R(C)`$ in order to obtain $`Z/Z_R`$.

Figure 3 shows $`\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R`$ as a function of the volume $`V=L^3`$ at the critical endpoint of the transition line (point $`E`$ in figure 1). Indeed, one finds an exponentially small signal, as expected from eq.(3.3). Defining a scale parameter $`L_0`$ that measures the severity of the complex action problem by

$$\mathrm{exp}(ih\underset{x}{}\text{Im}\mathrm{\Phi }_x)_R\mathrm{exp}\left(\frac{L^3}{L_0^3}\right)$$
(3.23)
we get $`L_080`$ at the endpoint $`E`$. This means that the complex action problem at the endpoint is extremely mild. In fact, in practice it is not a problem at all up to volumes as big as $`100^3`$. Since the current computer hardware restricts the system size to a couple of million degrees of freedom, one can also study the point $`E`$ with an algorithm that does not solve the complex action problem. On very large lattices the meron-cluster algorithm will become superior also at the point $`E`$ since its computational effort is polynomial in the system size as opposed to exponential for the reweighted Metropolis algorithm.

It should be noted that the complex action problem is most severe for intermediate values of $`h`$. While it is obvious that there is no complex action problem at $`h=0`$, it is perhaps less obvious that there is also no problem for large $`h`$. This is because

$$\frac{W_R(C)}{W(C)}=\frac{\mathrm{exp}(h|C|)+2\mathrm{exp}(h|C|/2)}{\mathrm{exp}(h|C|)+2\mathrm{exp}(h|C|/2)\mathrm{cos}(\sqrt{3}h|C|/2)}$$
(3.24)
approaches $`1`$ in the limit $`h\mathrm{}`$ so that $`\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R_R=Z_R/Z1`$. Figure 4 shows the complex action problem scale parameter $`L_0`$ as a function of $`h`$ for fixed $`\kappa =0.5495`$. It has a minimum at $`h1`$ meaning that the complex action problem ist most severe in that region. Defining the “practical complex action problem” as being present when $`\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R<0.01`$ on a $`100^3`$ lattice (i.e. $`L_0<60`$) we see that at $`\kappa =0.5495`$ there is no practical complex action problem for $`h<0.003`$ as well as for $`h>6`$. A more physical definition of a practical complex action problem would compare $`L_0`$ with the correlation length $`\xi `$ but we have not measured the correlation length.

## 4 Flux Representation of the Potts Model

In this section we map the Potts model to an equivalent flux model that does not suffer from the complex action problem. Then we describe a Metropolis algorithm to update the flux model and compare its efficiency with the cluster algorithm.

### 4.1 Mapping the Potts Model to a Flux Model

As pointed out in  the Potts model can be rewritten as a flux model that does not suffer from the complex action problem. As we have seen, the complex action problem can also be solved in the cluster formulation. Hence, the question arises if the flux or the cluster formulation leads to more efficient numerical simulations. Let us first match the flux model to the original Potts model. The flux model is formulated in terms of “electric” charges $`Q_x\{0,\pm 1\}`$ defined on the lattice sites $`x`$ and electric flux variables $`E_{x,i}\{0,\pm 1\}`$ living on the links. The charge and flux variables are related by the $`(3)`$ Gauss law constraint

$$Q_x=\underset{i}{}(E_{x,i}E_{x\widehat{i},i})\text{mod}3.$$
(4.1)
The action of the flux model takes the form

$$S[E,Q]=\frac{g^2}{2}\underset{x,i}{}E_{x,i}^2+\beta \underset{x}{}(MQ_x^2\mu Q_x).$$
(4.2)
The mass $`M`$ and the chemical potential $`\mu `$ of the $`(3)`$ charges are not directly related to the mass and chemical potential of quarks in QCD (which we also denoted by $`M`$ and $`\mu `$ in section 2) but qualitatively they play the same role. The partition function of the flux model takes the form

$$Z=\underset{x}{}\underset{Q_x\{0,\pm 1\}}{}\underset{x,i}{}\underset{E_{x,i}\{0,\pm 1\}}{}\underset{x}{}\delta _x\mathrm{exp}(S[E,Q]).$$
(4.3)
The $`\delta `$-function $`\delta _x`$ imposes the Gauss law of eq.(4.1) at the point $`x`$ and can be written as

$$\delta _x=\frac{1}{3}\underset{\mathrm{\Phi }_x(3)}{}\mathrm{\Phi }_x^{Q_x_i(E_{x,i}E_{x\widehat{i},i})}.$$
(4.4)
Inserting this as well as eq.(4.2) for the action in eq.(4.3) one can integrate out the $`E_{x,i}`$ and $`Q_x`$ variables. The result of the $`E_{x,i}`$ integration is

$$\underset{E_{x,i}\{0,\pm 1\}}{}(\mathrm{\Phi }_x^{}\mathrm{\Phi }_{x+\widehat{i}})^{E_{x,i}}\mathrm{exp}\left(\frac{g^2}{2}E_{x,i}^2\right)=1+2\mathrm{exp}\left(\frac{g^2}{2}\right)\text{Re}(\mathrm{\Phi }_x^{}\mathrm{\Phi }_{x+\widehat{i}}).$$
(4.5)
In the Potts model (up to an overall factor) the corresponding term is $`\mathrm{exp}(\kappa \delta _{\mathrm{\Phi }_x,\mathrm{\Phi }_{x+\widehat{i}}})`$. Thus, the flux model matches the Potts model if

$$\mathrm{exp}(\kappa )=\frac{1+2\mathrm{exp}(g^2/2)}{1\mathrm{exp}(g^2/2)}.$$
(4.6)
Note that the $`g0`$ limit of the flux model corresponds to the $`\kappa \mathrm{}`$ limit of the Potts model. When one integrates out the charges $`Q_x`$ one obtains

$$\underset{Q_x\{0,\pm 1\}}{}\mathrm{\Phi }_x^{Q_x}\mathrm{exp}(MQ_x^2\mu Q_x)=1+\mathrm{exp}(\beta (M\mu ))\mathrm{\Phi }_x+\mathrm{exp}(\beta (M+\mu ))\mathrm{\Phi }_x^{}.$$
(4.7)
In the original Potts model (up to an overall factor $`A`$) the corresponding term is $`\mathrm{exp}(h\mathrm{\Phi }_x)`$. Hence, the flux model matches the Potts model if

$`{\displaystyle \frac{1}{3}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{exp}(h\mathrm{\Phi })=A,`$
$`{\displaystyle \frac{1}{3}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }\mathrm{exp}(h\mathrm{\Phi })=A\mathrm{exp}(\beta (M+\mu )),`$
$`{\displaystyle \frac{1}{3}}{\displaystyle \underset{\mathrm{\Phi }(3)}{}}\mathrm{\Phi }^{}\mathrm{exp}(h\mathrm{\Phi })=A\mathrm{exp}(\beta (M\mu )).`$ (4.8)
These relations can be used to determine the parameters $`\mathrm{exp}(\beta (M\mu ))`$ and $`\mathrm{exp}(\beta (M+\mu ))`$ of the flux model in terms of the Potts model parameter $`h`$.

### 4.2 Metropolis Algorithm for the Flux Model and Comparison with the Cluster Algorithm

As first described in , the flux model can be updated with a simple Metropolis algorithm. One basic move in the algorithm creates or annihilates a nearest-neighbor charge-anti-charge pair across a given link. The other basic move creates or annihilates an electric flux loop around an elementary plaquette. These basic moves are proposed on every link and plaquette and are accepted or rejected in a Metropolis step. We have implemented both the Metropolis algorithm for the flux model and the cluster algorithm for the Potts model and we have verified that physical observables obtained with the two algorithms agree with each other.

The question arises which of the two algorithms is more efficient. The Metropolis algorithm is expected to suffer from critical slowing down at the endpoint of the first order phase transition with a dynamical critical exponent $`z2`$. Cluster algorithms are known to drastically reduce critical slowing down, in some cases even to $`z0`$. However, our cluster algorithm cannot eliminate critical slowing down completely because the decision to put a bond is more time-consuming than the one in the Swendsen-Wang algorithm. For example, when a bond is removed, one must check if an old cluster decomposes into two new clusters. To minimize the computational effort, we simultaneously grow two clusters from the two ends of the bond. Once they touch each other, we know that the old cluster did not decay. Still, in the less likely event that the old cluster does decay, one must completely grow the smaller of the two clusters in order to decide if the bond can be deleted. Here we do not attempt to determine the dynamical critical exponent $`z`$ of our cluster algorithm at the critical endpoint.

An even more severe super-critical slowing down is expected close to the first order phase transition line. Along that line, the deconfined phase coexists with the confined phase, and the Monte Carlo simulation must tunnel between the two phases. In order to tunnel between the confined and the deconfined phase, a local algorithm must go through configurations containing both phases simultaneously. Since the interface that separates the two phases has non-zero interface tension, such configurations are exponentially suppressed. Hence, close to a first order phase transition, a local algorithm like the Metropolis algorithm necessarily suffers from exponential slowing down. This super-critical slowing down is even more severe than the power-law critical slowing down at a second order phase transition. Hence, we expect that the Metropolis algorithm for the flux model is not very well suited to study the first order phase transition line. In some cases, cluster algorithms can even eliminate super-critical slowing down. For example, in the broken phase of the Potts model three distinct deconfined phases coexist with each other. The Swendsen-Wang algorithm can efficiently tunnel from one deconfined phase to another because it assigns the same random $`(3)`$ element to all spins in a cluster in a non-local spin update. Still, even the Swendsen-Wang algorithm suffers from super-critical slowing down at the first order phase transition that separates the confined from the deconfined phase. Although cluster flips can naturally lead to tunneling between distinct deconfined phases, they do not lead directly from a deconfined to the confined phase. To cure this problem, Rummukainen  has combined the cluster algorithm with the multi-canonical methods of Berg and Neuhaus  which can reduce the exponential super-critical slowing down to a power-law behavior. Although this may well be possible, we have not yet attempted to combine our algorithm with multi-canonical methods. Hence, we expect that our cluster algorithm still suffers from super-critical slowing down close to the deconfinement phase transition.

We have compared the efficiency of the Metropolis algorithm for the flux model, the cluster algorithm for the complex action Potts model, and the reweighted Metropolis algorithm for the complex action Potts model at several points in the phase diagram. In figure 5 we compare the computer time histories of the Polyakov loop for the three algorithms mentioned above at $`(h,\kappa )=(0.01,0.5)`$ which is in the confined region of the phase diagram. Obviously, the flux model Metropolis algorithm decorrelates a lot worse than the other two algorithms. The flux algorithm performs even worse when one approaches the first order transition line. The reweighted Metropolis algorithm for the Potts model suffers from the complex action problem. What is plotted is the time evolution of $`\mathrm{Re}[_x\mathrm{\Phi }_x\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)]/V\mathrm{exp}(ih_x\text{Im}\mathrm{\Phi }_x)_R`$. Its statistical fluctuations are much larger than the expectation value which is a manifestation of the complex action problem.

## 5 Universality Class of the Critical Endpoint

In this section we present the results of our numerical simulations at the critical endpoint. The complex action problem turned out to be very weak in the vicinity of the critical endpoint (see subsection 3.4). Therefore it is possible to use a simple reweighted Metropolis algorithm even though that does not solve the complex action problem. It is usable up to sufficiently large lattices so that critical exponents can be extracted from a finite size scaling analysis. Of course, on even larger lattices the meron-cluster algorithm will eventually be superior to the reweighted Metropolis algorithm. But since at the endpoint $`E`$ the complex action problem sets in only at volumes $`100^3`$, simulations at $`E`$ are not limited by the complex action problem but by the ability to simulate large lattices on today’s computers.

Figure 1 shows the phase diagram of the model defined by equation (2.11). For $`h=0`$ our model reduces to the standard 3-d 3-state Potts model which has been studied extensively in Monte Carlo simulations . The model is known to have a weak first order phase transition. The value of the coupling $`\kappa `$ where the phase transition occurs (point $`T`$ in fig. 1) has been determined with high precision. In , the phase transition was found to occur at $`\kappa _T=0.550565(10)`$. Above this value the $`(3)`$ symmetry is spontaneously broken, i.e. for $`\kappa >\kappa _T`$ three distinct deconfined phases coexist. When we switch on the parameter $`h`$, the $`(3)`$ symmetry gets explicitly broken. Positive values of $`h`$ favor the deconfined phase with a real value of $`\mathrm{\Phi }`$. Hence, the line $`\kappa >\kappa _T`$ at $`h=0`$ is a line of first order phase transitions which cannot terminate in the deconfinement transition at the point $`T`$. In fact, $`T`$ is a triple point because two other first order transition lines emerge from it. For $`h>0`$ a line of first order transitions extends into the $`(h,\kappa )`$-plane and terminates in a critical endpoint ($`E`$ in fig. 1). Negative values of $`h`$ favor the two deconfined phases with complex values of $`\mathrm{\Phi }`$. Negative $`h`$ are unphysical in the QCD interpretation of the Potts model because $`h`$ represents $`\mathrm{exp}(\beta (M\mu ))`$ in QCD. Still, the Potts model at $`h<0`$ makes perfect sense as a statistical mechanics system (unrelated to QCD) and it has another first order transition line emerging from the point $`T`$. Interestingly, with our method the complex action problem can only be solved for $`h0`$ since otherwise the improved estimator of eq.(3.8) is not necessarily positive. It should be noted that for $`h<0`$ also the flux model suffers from the complex action problem.

The line of first order phase transitions $`\kappa _t(h)`$ is determined by the condition that the free energy densities of the confined and deconfined phases are equal, i.e. $`f_c(h,\kappa _t(h))=f_d(h,\kappa _t(h))`$. Close to the point $`T=(0,\kappa _T)`$ the free energy density of the confined phase is given by

$$f_c(h,\kappa )=f_{c,T}+e_{c,T}(\kappa \kappa _T),$$
(5.1)
where $`f_{c,T}=f_c(0,\kappa _T)`$ and $`e_{c,T}=df_c/d\kappa (0,\kappa _T)`$ is the energy density of the confined phase at the point $`T`$. Note that to leading order $`f_c(h,\kappa )`$ is independent of $`h`$ because $`\mathrm{\Phi }=0`$ in the confined phase at $`h=0`$. On the other hand, for the deconfined phase one obtains

$$f_d(h,\kappa )=f_{d,T}+e_{d,T}(\kappa \kappa _T)h\mathrm{\Phi }_T,$$
(5.2)
where $`\mathrm{\Phi }_T`$ is the value of the Polyakov loop at the point $`T`$ in the deconfined phase that is favored at $`h>0`$. Using the condition $`f_{c,T}=f_{d,T}`$ for the deconfinement phase transition at $`h=0`$, one finds

$$\kappa _t(h)=\kappa _T\frac{\mathrm{\Phi }_T}{e_{c,T}e_{d,T}}h=\kappa _T\frac{h}{r}.$$
(5.3)
The Monte Carlo data of  and  imply $`r=0.41(1).`$ Our data are consistent with the first order transition line being a straight line. Fitting the values of $`\kappa _t(h)`$ obtained from the infinite volume extrapolation described below yields $`r=0.430(6)`$ in reasonable agreement with the number from above. Similarly, one can determine the angle at which the third transition line leaves the point $`T`$ in the direction of negative $`h`$. A similar argument can be applied to the Potts model with real action that was studied in . Also in that case the first order phase transition line is consistent with a straight line and again the predicted position for the transition line agrees with the numerical data. Interestingly, for the Potts model with both real and complex action, information at $`h=0`$ is sufficient to predict the position of the transition line for $`h>0`$. This is because the transition at $`h=0`$ is rather weak and the line ends already at small values of $`h`$. If the transition would extend deep into the $`(h,\kappa )`$-plane one would expect deviations from a straight line that would be hard to predict based on data at $`h=0`$.

To determine the location of the transition line numerically we perform for given values of $`h`$ and the volume $`V`$ simulations at 3 to 5 different values of $`\kappa `$. These simulations are then combined with Ferrenberg Swendsen multi-histogram reweighting . To estimate the position of the transition line we use the specific heat

$$C=\frac{1}{V}\left((S[\mathrm{\Phi }]h\underset{x}{}\mathrm{\Phi }_x)^2(S[\mathrm{\Phi }]h\underset{x}{}\mathrm{\Phi }_x)^2\right)$$
(5.4)
and determine the position of its maximum $`\kappa _t(h,V)`$ for a given $`h`$ and $`V`$. The transition point $`\kappa _t(h)`$ is determined in the infinite volume limit using

$$\kappa _t(h,V)=\kappa _t(h)+\frac{A(h)}{V}$$
(5.5)
where $`\kappa _t(h)=\kappa _t(h,V=\mathrm{})`$. This ansatz is used sucessfully in the whole $`h`$-range, i.e. for the first order region as well as in the crossover region. The values for $`\kappa _t(h)`$ are plotted with errorbars into the phase diagram (figure 1). On the first order transition line they are consistent with a straight line which intersects the $`\kappa `$ axis exactly at the point $`T`$. The crossover line has a slight curvature.

After we have determined the tangent to the first order transition line close to the endpoint we still have to find the exact location of the endpoint on that line. Also we want to extract critical exponents. For that purpose we consider the Polyakov loop susceptibility

$$\chi =\frac{1}{V}\left((\underset{x}{}\mathrm{\Phi }_x)^2\underset{x}{}\mathrm{\Phi }_x^2\right)$$
(5.6)
along the transition line, i.e. at the points $`(\kappa _t(h),h)`$. One could consider a “rotated susceptibility” instead, where an admixture of the kinetic energy term is added to the Polyakov loop to diagonalize the fluctuation matrix. Nevertheless, this is not necessary, since the “magnetic field direction” is the dominant one. We explicitly checked that the results do not depend on the admixture we chose, unless one comes close to the linear combination where the discontinuities of the Polyakov loop and the kinetic energy cancel. This linear combination would correspond to the kinetic energy term in the Ising model. Of course, if one would be interested in observables related to the kinetic energy of the Ising model, one would have to exactly choose the linear combination that corresponds to the kinetic energy direction.

Close to the critical point the following scaling ansatz describes the susceptibility (see e.g. )

$$\chi =L^{\gamma /\nu }f(x),x=(hh_c)L^{1/\nu }.$$
(5.7)
For the fit, the function $`f(x)`$ is expanded in a Taylor series around $`x=0`$ up to second order, $`f(x)=f_0+f_1x+f_2x^2/2`$. We perform two fits:

* Six parameter fit: It results in $`1/\nu =1.532(57)`$, $`\gamma /\nu =2.064(70)`$, $`h_c=0.000445(18)`$, $`f_0=0.70(16)`$, $`f_1=3.99(95)`$, and $`f_2=16.4(60)`$ with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.=1.33`$. The results for the exponents agree with the estimates for the 3d-Ising universality class $`1/\nu =1.587(2)`$ and $`\gamma /\nu =1.963(3)`$  almost within errorbars.
* Four parameter fit: The critical exponents are fixed to the Ising model values. The result is $`h_c=0.000470(2)`$, $`f_0=0.9775(46)`$, $`f_1=4.567(55)`$, and $`f_2=16.5(17)`$ with $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.=1.35`$. The good value for $`\chi ^2`$ supports again the universality class of the 3-d Ising model.

The susceptibility and the four parameter fit are shown in Figure 6. Figure 7 shows the function $`f(x)`$.

## 6 Conclusions

We have used a cluster algorithm to solve the notorious complex action problem in the Potts model approximation to QCD with heavy quarks at large chemical potential. We use a simple analytically constructed improved estimator that gives an exponential reduction in the required statistics. Since the improved estimator is real and positive, importance sampling techniques that fail for complex actions then become applicable. This makes it possible to study the whole $`h>0`$ parameter range of the Potts model, not just the $`h=0`$ axis. (Recall that $`h`$ corresponds to $`\mathrm{exp}(\beta (\mu M))`$ in QCD in the limit $`M,\mu \mathrm{}`$ at any given $`\mu M`$).

We compared our cluster algorithm with a flux model reformulation, and with the reweighted Monte Carlo algorithm. We found that the cluster algorithm was more efficient than using the flux model reformulation. In the large volume limit the cluster algorithm will always be superior to Monte Carlo reweighting. However, at very small $`h`$ it is sufficient to use reweighting techniques to obtain physically relevant results. This turned out to be the case for the endpoint of the first order line, which occurs at a very small $`h`$ because the 3-d 3-state Potts model phase transition is rather weak at $`h=0`$. We therefore used reweighted Monte Carlo to locate the first-order line and its endpoint. However, we emphasize that as computer power rises, and the maximum attainable volume becomes bigger, the meron-cluster algorithm will eventually become superior at any $`h>0`$.

We also calculated quark-quark, quark-antiquark, and antiquark-antiquark potentials, in the confined and deconfined regions of the phase diagram. We found the expected behavior: the background density of heavy quarks screens color fields, so that all potentials reach plateaux at long distances, whose values are simply related to the free energies of external static quarks and antiquarks.

The algorithm that we have developed for the Potts model belongs to the class of meron-cluster algorithms that has recently been used to solve a large variety of sign and complex action problems. Of course, the ultimate goal is to construct a similar algorithm for QCD at non-zero chemical potential and investigate the phase structure of QCD at $`\mu 0`$ from first principles. The complex action problem in full QCD is more complicated than the one in the Potts model. So far, meron-cluster algorithms have led to solutions of fermion sign problems as well as complex action problems in bosonic theories, but have not yet solved complex action problems in theories with fermions. We believe that this may ultimately become possible when one uses the D-theory formulation of QCD.

## Acknowledgements

We like to thank the INT in Seattle, where this work was initiated, for its hospitality. U.-J. W. thanks F. Karsch and S. Stickan for helpful discussions and acknowledges the support of the A. P. Sloan foundation.