Canonical fermion determinants in lattice QCD - Numerical evaluation and properties

We analyze canonical fermion determinants, i.e., fermion determinants projected to a fixed quark number q. The canonical determinants are computed using a dimensional reduction formula and are studied for pure SU(3) gauge configurations in a wide range of temperatures. It is demonstrated that the center sectors of the Polyakov loop very strongly manifest themselves in the behavior of the canonical determinants in the deconfined phase, and we discuss physical implications of this finding. Furthermore the distribution of the quark sectors is studied as a function of the temperature.


Introductory remarks
With the already running and planned experiments at RHIC, LHC and FAIR, Quantum Chromodynamics with temperature and non-vanishing density has recently shifted into the center of attention. Of particular interest are nonperturbative results focusing on various aspects, such as obtaining information about the phase diagram, characterization of the different phases, and a possible understanding of the mechanisms that drive the various transitions. In principle lattice QCD is such a non-perturbative approach, but numerical simulations at finite density are plagued by serious phase cancellation problems.
There exist two alternative approaches to QCD with finite density: One can work grand-canonically using a chemical potential, or alternatively consider the canonical ensemble, i.e., a fixed number of quarks. The majority of lattice calculations at finite density uses a chemical potential, while the canonical approach has seen comparatively little attention [1] - [6].
In this article we revisit the canonical approach. Using a recently proposed dimensional reduction for the fermion determinant [6] we can numerically evaluate canonical determinants much faster, i.e., fermion determinants describing a fixed number of quarks. We combine the dimensional reduction with different strategies for computing canonical determinants and compare their efficiency.
Canonical determinants are not only objects that might provide new and more efficient strategies for lattice QCD with finite density, but they also have interesting properties themselves. In particular they show simple transformation properties under center transformations -a symmetry that is of crucial importance for understanding the high temperature phase transition of pure gauge theory. Thus the numerical analysis of canonical determinants might provide important insight into mechanisms responsible for transitions in the QCD phase diagram. In this article we demonstrate that the influence of center symmetry and its breaking are very strongly manifest in the properties of the canonical determinants in the high temperature phase and we discuss the physical implications of this finding. Furthermore, we show that in the high temperature phase the canonical determinants with odd quark numbers become suppressed.

Canonical partition sums on the lattice
In this section we present our lattice Dirac operator and compile formulas for canonical partition sums and canonical determinants. This serves to set our notational conventions and to summarize the formulas which we later need for the analysis of the canonical determinants.

Wilson fermions with chemical potential
In this article we use Wilson's lattice Dirac operator. The chemical potential is introduced as a boundary condition, such that the lattice Dirac operator reads The lattice indices x = ( x, x 4 ), y = ( y, y 4 ) run over the L 3 × β lattice, byĵ we denote the shift vector in j-direction, and we use the convention U −j (x) = U j (x −ĵ) † . β is the inverse temperature in lattice units, i.e., β is a positive integer. The temporal lattice index x 4 runs from −β/2 + 1 to β/2. In our notation the hopping parameter κ is related to the bare quark mass m via κ = 1/(4 + m). The chemical potential µ is attached to the link terms connecting x 4 = β/2 and x 4 = −β/2 + 1 in the usual exponential form. These terms also come with an additional relative sign which implements the anti-periodic temporal boundary conditions for the fermions. The link variables U ν (x) are elements of the gauge group SU (3). The corresponding color indices, as well as the Dirac indices are not displayed explicitly and we use vector/matrix notation for them. For later use we recollect two important properties of the lattice Dirac operator D(µ) and its determinant: The first one is its behavior under a similarity transformation with γ 5 (µ is real here), which often is referred to as generalized γ 5 -hermiticity. It implies for the determinant of the Dirac operator the identity where the asterisk denotes complex conjugation.
The second property reflects the behavior under a time reflection of the fields, and is conveniently expressed by the identity (4) On the left-hand side the determinant is evaluated for the time reflected gauge field T U (see, e.g., [7] for a detailed definition of the time reflection), while on the right-hand side the original gauge field configuration U is used. Combining (3) and (4) gives rise to

Definition of canonical partition sums
This subsection gives a short summary of the definition of canonical partition sums to set our notation. For N mass degenerate flavors of fermions, the grand canonical partition sum Z GC (µ) at chemical potential µ is given by In the following we will often refer to the fermion determinant det[D(µ)] or its N -th power as the grand canonical determinant. By X G we denote the path integral expression for some observable X in pure gauge theory, where S G [U ] is the gauge action and D[U ] denotes the path integral measure over all gauge fields, which on the lattice is given by the product of the Haar measures for all link variables U ν (x). For the grand canonical partition sum one may perform a fugacity expansion, i.e., write it as a sum over canonical partition sums Z (q) with a fixed quark number q ∈ Z, As a matter of fact, for a finite L 3 × β lattice the sum in the fugacity expansion runs only over the finite interval q ∈ [−q max , q max ], where, e.g., for one flavor of Wilson fermions one has q max = 6L 3 . The canonical partition sums Z (q) with a fixed quark number q may themselves be written as path integrals, where the canonical fermion determinant with a fixed quark number q is denoted as det N [D] (q) . We attach the number of flavors N as a lower index and omit the argument µ of the Dirac operator, since it has disappeared from the canonical expressions and enters only through the expansion coefficients exp(µqβ) in the fugacity expansion (8).
We stress at this point, that the canonical N -flavor determinant det N [D] (q) with fixed quark number q does not factorize into one-flavor canonical determinants. This is different from the grand canonical determinant which, as seen in, e.g., (6), is simply obtained by taking the N -th power of the grand canonical one-flavor determinant det[D(µ)].
We remark that once an efficient way of computing the fermion determinant det N [D] (q) for a fixed quark number q is available, the expression for the corresponding canonical partition sum Z (q) is rather simple. Eq. (9) shows that up to a normalization it can be written as an expectation value ... G in pure gauge theory.
A possible way of computing det N [D] (q) is through a Fourier transformation with respect to imaginary chemical potential, The explicit form (1) of the lattice Dirac operator shows that also imaginary chemical potential is just a boundary condition which multiplies the forward running temporal link variables of the last time slice with the phase e iϕ , while the backward running temporal link variables on the last time slice are multiplied with e −iϕ . Combining Eqs. (6) -(10), one obtains the fugacity expansion of the grand canonical determinant (det[D(µ)]) N in terms of the canonical determinants det N [D] (q) . Equation (11) may now be used to obtain an expression of the canonical determinants with N > 1 flavors in terms of the canonical one-flavor determi-nants. Exploring (11) for N = 1 one finds where δ(q 1 + ... + q N , q) denotes a Kronecker delta which enforces q 1 + q 2 + ... + q N = q. Comparison with (11) allows one to read off the expression which we will later use to compute the canonical determinants det N [D] (q) for the case of N > 1 flavors from the one-flavor canonical determinants det 1 [D] (q) . We conclude this subsection with exploring the implications of γ 5 -hermiticity (2) and the time reflection symmetry (4) for the canonical determinants: Combining γ 5 -hermiticity in the form of (3) with the fugacity expansion (11), one finds the relation between canonical determinants with positive and negative quark numbers q. For vanishing quark number q = 0 the canonical determinant is real, i.e., In a similar way we use (4) together with the fugacity expansion (11) and arrive at det which shows that under a time reflection of the gauge configuration canonical determinants turn into their complex conjugate. Since in the path integral ... G over all gauge fields both, a configuration U and its time reflected image T U appear with the same weight, we find i.e., the canonical partition sums Z (q) are real. Via the fugacity expansion (8) this implies that also the grand canonical partition sum Z GC is real.

The fermion loop picture for canonical determinants
The canonical fermion determinants det N [D] (q) have an interesting interpretation as collections of fermion loops with a fixed total winding number q: As a gauge invariant quantity on the lattice, the grand canonical fermion determinant (det [D]) N can be decomposed into sums over products of closed loops on the lattice dressed with the gauge variables on the links the loops consist of. One finds In our notation s runs over sets S of closed loops, c s is a complex valued coefficient, and l runs over the individual loops in the set s. The second product is over the links (x, µ) the loop l is made of, and the trace is over color indices. One may, e.g., use hopping expansion to compute some terms in (17) explicitly.
Since the loops l are closed, one can determine their winding number q l around compactified time with q l ∈ Z (negative values correspond to backward winding). An individual term in the sum (17), which consists of a set s of loops, may then be assigned a total winding number q s = l∈s q l .
When one uses the definition (10) of the canonical determinants det N [D] (q) , a boundary condition phase e ±iϕ is attached to the temporal link variables on the last time slice. Thus a loop l that winds q l times picks up a factor of e iq l ϕ . An individual term in the sum (17) which is made from a product of loops in a set s thus acquires a factor of e iqsϕ . The Fourier integral in (10) projects to a fixed total winding number q s = q, such that for the canonical determinants we obtain the loop representation where S (q) now contains only sets s of loops that have a total winding number q s = q. Thus the canonical determinant det N [D] (q) is made of products of loops with a total winding number of q.
2.4. Center transformation properties of the canonical determinants SU(N) gauge theory has an interesting property which plays an important role for the QCD phase transitions, namely a symmetry under center transformations.
In lattice language the center transformation consists of multiplying all temporal link variables at a fixed time slice x 4 = x * 4 with a fixed element z of the center of the gauge group, while all other link variables remain unchanged. For the case considered here, i.e., gauge group SU(3), the center is Z 3 and thus z has the three possible values z ∈ 1 , e i2π/3 , e −i2π/3 .
The gauge action S G [U ] is invariant under such a transformation, not, however, the grand canonical fermion determinant. According to Eq. (18) the canonical determinants det N [D] (q) can be written as sums over products of loops, where each term has a total combined winding number of q for the loops the term consists of. Since the transformation (19) affects only a single time slice, the total factor that is picked up by the loops is given by z q . Thus we conclude that under the center transformation (19) the canonical determinants transform as i.e., they pick up a phase which depends on the quark number q. Since z may only assume the values listed in (20), which form the center group Z 3 , the powers z q are also restricted to these values. In the last expression of (21) we have stressed this fact by introducing the concept of triality, defined as q mod 3. Equation (21) shows that only the zero triality determinants remain invariant, while for non-zero triality a phase is picked up under the center transformation.
The center symmetry has a simple physical consequence: Performing a center transformation of the integration variables U in the path integral (7) leaves both the gauge action and the integration measure invariant. Thus as long as the center symmetry is not broken spontaneously in ... G , the only non-trivial transformation of an expectation value can come from the integrand inserted in ... G . For the canonical partition sums this leads to which implies that Z (q) must vanish unless z q = 1. Thus only if the canonical partition sum Z (q) is of zero triality, i.e., q mod 3 = 0, it is non-vanishing.
To include this finding in the fugacity expansion (8), we write the canonical partition sum as where we take into account only the non-vanishing zero triality terms with q = 3b by introducing the baryon number b ∈ Z. We stress again, that this chain of arguments holds only if the center symmetry is not broken spontaneously by ... G . The impact of a spontaneous breaking will be discussed in Section 5.

Techniques for computing canonical fermion determinants
In the last section we have addressed two possible ways of computing the canonical determinants det N [D] (q) . Equation (10) defines them through a Fourier transformation, while (18) formally expresses them as a sum over loops with a fixed total winding number. Both equations are not very suitable for an immediate application in a numerical evaluation, and in this section we discuss various strategies for an efficient evaluation of canonical determinants. We begin this discussion with a presentation of a dimensional reduction formula [6] for the determinant that will be a powerful tool for all approaches which we subsequently discuss. We will restrict ourselves to the simple case of only one flavor, since the case of more than one flavor, i.e., N > 1, can be reduced to N = 1 with the help of (13).

Dimensional reduction of the fermion determinant
For later use we summarize the essential formulas for an exact dimensional reduction of the fermion determinant [6]. For that purpose the lattice is decomposed into 4 domains, Λ (i) , i = 1, 2, 3, 4, which are distinguished by the values of the temporal lattice index x 4 . In particular Λ (1) consists of all lattice points in the range from The domains Λ (2) and Λ (4) consist of single time slices with x 4 = 0 and x 4 = β/2, respectively. The domain decomposition of the lattice is illustrated in Fig. 1. Based on the domain decomposition the Dirac operator is split into several pieces: The terms D (i) , i = 1, 2, 3, 4, live on the individual domains Λ (i) , i = 1, 2, 3, 4. The terms D (1,2) , D (2,3) , D (3,4) and D (4,1) contain the forward hopping temporal terms that connect the domains, and D (2,1) , D (3,2) , D (4,3) and D (1,4) are the corresponding backward hopping terms. The fermion determinant may be factorized exactly [6] in the form where the terms H 0 , H ±1 live on only the single time slice that constitutes the domain Λ (4) . They are made out of the Dirac operator terms D (i) , D (i,j) and the corresponding propagators, where The subscript 0 attached to the terms D indicates that there the chemical potential was set to 0, and instead is taken onto account by the exponential factors attached to the terms H ±1 in Eq. (24). The overall factor A 0 in (24) is independent of the chemical potential µ and given by We stress again that the formula (24) is an exact factorization that does not contain any approximation and holds for arbitrary background gauge fields. It provides a powerful tool for the analysis of lattice QCD with chemical potential: • Eq. (24) factorizes the piece that depends on the chemical potential µ and reduces it to a determinant which lives on a single time slice. Due to this dimensional reduction the numerical evaluation of the µ-dependent piece is considerably less demanding in a numerical analysis.
• The µ-independent and real factor A 0 was shown [6] to carry most of the scale, i.e., the overall size of the determinant which rapidly grows with the volume and shows a strong quark mass dependence. The second, µ-dependent determinant in (24) remains small compared to A 0 , and in isolated form contains the terms that are problematic in lattice QCD with chemical potential, i.e., the complex phase.
• Finally the factorization (24) allows one to compute the canonical determinants using a generalized hopping expansion with excellent convergence properties. We will elaborate on this approach below.

Fourier projection method
The most straightforward method for the determination of the canonical determinants is the Fourier projection as shown in Eq. (10). As a matter of fact it is sufficient to apply this formula to the one-flavor case N = 1, since the determinants for N > 1 may be built up using (13). A problem of the Fourier projection method is that for larger values of q the Fourier integral has to be computed very precisely. For the numerical evaluation this implies that the determinant has to be computed for many values of ϕ in the interval [−π, π]. Since the numerical cost for the evaluation of a determinant increases with the size N of an N × N matrix like N 3 , this is a rather expensive approach.
The situation improves considerably when the factorization formula discussed in the last subsection is applied [6]. Using (24) together with (10) one finds In this expression the determinant is for matrices that live on only a single time slice. Thus the size N of the matrix is reduced by a factor of 1/β and the cost for computing the determinant goes down by a factor of 1/β 3 . For lattice sizes typically used in finite-µ calculations one has β = 4 -6, such that the speed up of the calculation is of order O(100). Thus, at a fixed numerical cost a much larger number of intermediate values of ϕ becomes available and Fourier projection becomes an attractive approach also for larger values of q.

Generalized hopping expansion
It may be shown that the terms H 0 , H ±1 consist of paths of at least β links, dressed with link variables [6]. Consequently a product of n terms H 0 , H ±1 is made of paths with a length of at least nβ links. Since products of link variables along such paths decrease in amplitude roughly exponentially with the length of the path, this property can be used to set up an efficient expansion [6] of the determinant, which may be viewed as a generalization of the usual hopping expansion. Using det[M ] = exp(Tr ln M ) one finds where we have introduced For a perturbative expansion the sum in the exponent of Eq. (29) has to be cut off such that q is restricted to q ∈ [−Q, Q] (see last expression in (29)). The exponential in the last line of (29) may then be expanded and the resulting series can be organized with respect to powers of exp(±µβ). Comparison with (11) for the case of N = 1 allows one to read off the canonical determinants det 1 [D] (q) as the coefficients of the terms exp(µqβ). In a practical application also the infinite sum (30) has to be truncated at some finite value of n. In [6] it was demonstrated that the contributions decrease exponentially with increasing n, such that the series (30) converges quickly and the effect of the truncation can be controlled systematically. One finds that the results for det 1 [D] (q) from the expansion are reliable up to |q| ≤ Q (see also Section 4).

Generalized Fourier method
Yet another method for the evaluation of the canonical determinants was suggested in [2]. The basic idea is to expand the logarithm of the determinant into its Fourier components and then to expand the exponential of this series. In our presentation here we directly apply the method to the dimensionally reduced determinant (24). Evaluating this expression at imaginary chemical potential, we define the expansion coefficients c n , where c 0 is real, while the c n with n > 0 are complex and the asterisk again denotes complex conjugation. These coefficients may, e.g., be computed by a numerical evaluation of the corresponding Fourier transforms of ln det An alternative strategy for the determination of the c n is provided by the generalized hopping expansion discussed in the last section. A comparison of (31) with the last expression of (29) shows that the coefficients c n can be identified as where the T (n) can be computed using (30). Inserting (31) into the basic definition (10) at N = 1 gives rise to where we have cut off the sum over all Fourier coefficients at some maximal value Q. Describing the coefficients by their absolute value and phase, introduced as we find upon inserting this form in (33) Writing the cosines in the last factors again as exponentials, one can solve the remaining ϕ-integrals explicitly using where the I m denote the Bessel functions of the first kind. As for the generalized hopping expansion, the necessary truncations used in the generalized Fourier method may be controlled systematically due to the decrease of the coefficients c n (i.e., the −T (n) ) with increasing n.

Numerical evaluation of canonical fermion determinants
Having presented various strategies for the numerical evaluation of the canonical determinants in the last section, we now briefly discuss our experience with the different approaches. The tests were done for the free case, as well as for SU (3) configurations from the pure gauge ensemble (7). We begin with briefly describing these gauge ensembles in the next subsection and subsequently present the outcome of our comparison.

The gauge ensembles
It was already stressed in Section 2.2 that in the canonical approach the canonical partition sums Z (q) may be expressed (9) using the expectation value ... G of pure gauge theory as given in (7). In this paper we perform various numerical tests using such pure gauge theory ensembles. These tests are done on three different lattice volumes, 6 3 × 4, 8 3 × 4 and 10 3 × 4. For the former two we use ensemble sizes of 500 (1000 for two ensembles) configurations, while for the largest lattice 100 configurations were evaluated. We work with the Lüscher-Weisz gauge action [8] and for an estimate of the scale we in some places use the lattice spacing determined in [9] from the Sommer parameter. The update was implemented using a mix of overrelaxation and Metropolis sweeps. Our update furthermore contains a random center rotation of the configurations. As discussed, this is always a symmetry of the gauge theory as long as the system size is finite. Since we are interested in the behavior of the canonical determinants in all three Polyakov loop sectors (compare Section 2.4), we decided to include a random center rotation in the update in order to generate an equal share of configurations in each sector. For the two smaller lattice sizes we use several different values of the gauge coupling β G 1 to scan the low and high temperature regimes. For our lattices with time extent β = 4, the critical value of the gauge coupling is near β G = 7.8. We use the analysis of the critical temperature T c for the Lüscher Weisz gauge action in [10] to determine the ratio T /T c for our ensembles, which we list in Table 1, together with other properties of the gauge configurations.
For later use and in order to illustrate properties of our configurations, in Fig. 2 we show a distribution of the values of the Polyakov loop in the complex plane. It is obvious that for low temperature (i.e., the smaller values of the gauge coupling β G ) the Polyakov loop values scatter isotropically near the origin in the complex plane. For high temperature the well known anisotropical distribution starts to form: The data points move away from the origin and the phases θ P of the Polyakov loop are close to the three center values, i.e., θ P ≃ 0, θ P ≃ ±2π/3. As long as the lattice is finite, the Polyakov loop sectors corresponding to these three values are completely equivalent. In the infinite volume limit, the underlying center symmetry (see Section 2.4) is broken spontaneously in pure gauge theory at high temperature and the system selects one of the three sectors.
In most cases the quark mass m in the Dirac operator was set such that it corresponds to 100 MeV on the scale of the pure gauge theory. For two of our ensembles we also studied additional quark masses of 50 and 200 MeV. The parameters of our gauge ensembles are summarized in Table 1.
The determinant calculations and matrix inversions necessary for computing   configurations could dominate the result. The situation might be improved by implementing, e.g., a reweighting approach in a fully dynamical simulation. However, already the approximation allows for interesting insights into the behavior of the canonical determinants and we expect that in particular the center symmetry properties will carry over to simulations with more advanced sampling techniques.

Comparison of the numerical approaches to canonical determinants
In order to assess the three different strategies for computing the canonical determinants which we discussed in Sections 3.2 -3.4, we compared the three methods in the free case and for pure gauge theory ensembles on the 6 3 × 4 and 8 3 ×4 lattices. In all three approaches we made use of the dimensional reduction formula (24). Since this is an exact expression which may be evaluated very precisely [6], this step does not introduce any approximation, but considerably speeds up all three approaches. Using the standard Fourier method of Section 3.2 one can in principle evaluate the canonical determinants for arbitrary large q, as long as the numerical integration scheme used for the Fourier transformation is accurate enough. In our implementation we experimented with the trapezoidal integration scheme, as well as Simpson's rule, and found almost no difference for the cases we studied. Using a sufficiently large number N ϕ of intermediate values of ϕ ∈ [−π, π], one may compute the canonical determinants up to q = 10 with almost arbitrary precision (on 8 3 × 4). We considered N ϕ up to N ϕ = 4096, which for the smaller values of q is certainly a very conservative number. On the other hand, for the larger values of q such an order of magnitude for N ϕ is required for sufficient accuracy. By increasing the number N ϕ of points used for the integration one can systematically improve the accuracy for each volume and set of parameters β G , m, q. Thus we consider the results from the Fourier approach as our standard values where we compare the other results to.
An independent check of the accuracy for the numerical evaluation of the canonical determinants is obtained by summing the fugacity expansion (11). As we will discuss in Section 5, for the free case and also for our ensembles at the lower values of β G , we found an excellent agreement between the sum and the exact expression for the grand canonical determinant. For larger β G higher orders of q are needed (compare Fig. 8) which were not calculated in this work.
The other two methods, the generalized hopping expansion (Section 3.3) and the generalized Fourier method (Section 3.4) have a cutoff parameter which plays the same role in both approaches. This is the number Q in Equations (29) and (33), respectively. The parameter Q sets the cutoff up to which winding number the terms in the logarithm of the determinant are taken into account. The actual calculation of the terms is, however, different in the two approaches: In the generalized hopping expansion they are obtained as sums over traces of products of H 0 , H ±1 , while for the generalized Fourier method they are again obtained through a numerical Fourier integral.
In order to obtain the canonical determinants from the expansion of the logarithm of the grand canonical determinant, also the exponential has to be expanded. For the two methods described in Sections 3.3 and 3.4 this second expansion is done in different ways, but in principle, both methods allow to formally calculate canonical determinants det N [D] (q) for arbitrary large q in the second expansion, irrespective of the cutoff value Q. However, in our comparison we found that for both methods sizable deviations from the results of the standard Fourier method appear as one increases q beyond the cutoff Q. This finding does not really come as a surprise, since when expanding the exponential, the leading term for det N [D] (q) comes from the term T (q) in the generalized hopping method and from the term c q in the generalized Fourier method. Thus for q > Q the leading terms are missing and both methods cannot be used reliably for such values of q. On the other hand, as long as the bound q ≤ Q is respected we found that both methods agree well with the results from the standard Fourier method.

Properties of canonical fermion determinants
Having discussed and compared the possible strategies for the evaluation of the canonical fermion determinants, we now analyze their properties numerically. Except for the discussion of the results for free fermions in the next subsection, throughout this section the canonical determinants were evaluated on our pure gauge ensembles using the generalized hopping expansion method for quark numbers ranging from q = −6 to q = 6.

The case of free fermions
We start our discussion of canonical determinants with the case of free fermions. The study of the free case allows for a precise assessment of the three approaches discussed in Sections 3.2 -3.4, and furthermore provides an interesting illustration of the properties of canonical determinants. In this situation it is rather straightforward to evaluate the grand canonical determinant and thus also the canonical determinants. For the free case Fourier transformation allows one to compute explicitly all eigenvalues of the Dirac operator (1) for arbitrary boundary conditions. The product of the eigenvalues gives the grand canonical determinant and the canonical determinants are obtained easily using the Fourier transformation (10).
It is interesting to analyze how the canonical determinants det 1 [D] (q) are distributed as a function of q. For that purpose we consider for one flavor the ratio det 1 [D] (q) / det 1 [D] (0) , i.e., we normalize with respect to the trivial sector q = 0. More generally, in the sum for the fugacity expansion the canonical determinants are weighted with exp(µqβ). In order to directly illustrate the behavior of these terms, in Fig. 3   The µ = 0 data in Fig. 3 show that the distribution of the canonical de- Changing the value of the mass parameter has only a relatively small effect, with the distribution becoming wider as the quark mass in decreased. Increasing the spatial volume at fixed β leads to a wider distribution. If both the spatial volume and the temporal extent β are increased at a fixed aspect ratio L/β and fixed mass m, the width of the distribution remains essentially invariant.

Scatter plots of the canonical determinants
After illustrating the properties of the canonical determinants in the free case, we now switch to our gauge ensembles. We begin our analysis  Figure 4 shows the scatter plots for the ensemble at β G = 7.4, i.e., in the center symmetric phase. For q = 0 the canonical determinant is real (see Section 2.2), while for q > 0 the values scatter in the complex plane. For the two cases with vanishing triality, q = 3 and q = 6, one expects a nonvanishing expectation value det 1 [D] (q) G according to the discussion of the center symmetry in Section 2.4. Thus there must be an anisotropy in the distribution in the complex plane. However, such an anisotropy is hardly visible for q = 3 and q = 6, which illustrates that computing, e.g., the free energies for baryonic states in the confining phase from the corresponding canonical determinants is heavily plagued by phase cancellations.
In Fig. 5 we now show the scatter plots for an ensemble in the phase with broken center symmetry (β G = 8.2). The breaking of the center symmetry manifests itself in much more organized scatter plots. Again the q = 0 canonical determinant is exactly real, but now the other two cases with vanishing triality, q = 3 and q = 6, also have their data points concentrated near the real axis. The other cases, q = 1, 2, 4, 5, show the angular distribution familiar from the Polyakov loops (compare Fig. 2), where the data points have phases close to the center values 0, 2π/3 and −2π/3. The pattern is most pronounced for small values of q and becomes less clean as q is increased. This behavior can be understood due to the fact that the number of windings of the loops building up the canonical determinants is given by q: Although the canonical determinants always transform according to (21), loops that wind many times in order to give rise to large q amplify the noise introduced by the individual gauge configuration, and the pattern reflecting the center transformations becomes more blurred with increasing q. Changing the mass has essentially no influence on the center pattern.
It would be very interesting to analyze how the sharpness of the center pattern changes as one goes to larger and finer lattices (a question beyond the scope of this study). The experience from the Polyakov loop leads one to expect that the pattern becomes sharper as the volume is increased and the resolution of the lattice becomes finer. Anyway, it is surprising that for the smaller values of q the center pattern is so clearly pronounced, much more pronounced than for the Polyakov loop as may be seen by comparison with Fig. 2.
In Section 2.2 we discussed the relation between canonical one-flavor determinants and the canonical determinants for N > 1 flavors. For completeness in Fig. 6 we show a few scatter plots of canonical determinants in the two-flavor case as computed from (13). In particular we show det 2 [D] (q) for q = 0, 1, 2 (left to right). The top row is for an ensemble in the center symmetric phase (β G = 7.4), while in the bottom row we use a high temperature ensemble (β G = 8.2).
It is obvious that for N = 2 qualitatively the behavior is the same as in the one-flavor case. In both phases the q = 0 determinants are real. In addition for the case of two flavors the values are strictly non-negative as may be seen from (13) and (14). In the center symmetric phase the values for q = 1 and q = 2 scatter isotropically in the complex plane, while for the center broken phase we again find the familiar center pattern. The amplitudes are roughly the squares of the amplitudes for the one-flavor case, as is expected from (13) where the two-flavor case is expressed as a sum of squares of one-flavor determinants.

Implications of the center properties
In the last section we have illustrated that the center pattern is surprisingly strongly pronounced in the high temperature phase. As a consistency check whether we understand the connection between the sectors of the Polyakov loop as shown in Fig. 2 and the center pattern in Fig. 5, we divide the 1000 configurations in our β G = 8.2 ensemble into three subsets according to the phase θ P of the Polyakov loop. In each of the three sectors, the real (θ P ≃ 0), and the two complex sectors (θ P ≃ ±2π/3), we have roughly a third of the total statistics.
In Fig. 7 we show scatter plots for the canonical q = 1 (top row) and q = 2 (bottom) determinants. In the lhs. column we display the sector with real Polyakov loop, followed by two columns showing the two complex sectors.
The figure nicely demonstrates, that for q = 1 the center angle follows exactly the angles determining the sectors of the Polyakov loop. For q = 2 the phase of det 1 [D] (2) is complex conjugate compared to the Polyakov loop phase because the two-flavor determinant transforms under center rotations with z 2 = z * . Fig. 7 thus confirms that the center pattern is manifest as expected from Eq. (21).
Let us now discuss an important consequence of the strong center pattern observed in the canonical determinants. The transition from the real Polyakov loop sector to the sector with θ P ≃ 2π/3 may be implemented by a center rotation of the gauge configuration with z = exp(i2π/3) ≃ exp(iθ P ), and equivalently for the other complex Polyakov loop sector. According to Eq.   (37) is a sum of essentially real and positive terms. The fact that for θ P = 0 the det 1 [D] (q) are essentially real and positive follows from the scatter plots in Fig. 5 and 7. The complex Polyakov loop sectors, on the other hand, are characterized by θ P ≃ 2π/3. As discussed in the last paragraph, these phases can also be obtained from the real sector by a center rotation with z = exp(±i2π/3).
Inserting these values of z into (37) leads to relative complex phases in the fugacity expansion and thus to cancellations. Consequently one expects that the values of the canonical determinants are smaller for the complex Polyakov loop sectors.
In order to study this effect we analyzed the expectation value of the modulus of the grand canonical determinant, | det[D(µ)]| G , again dividing the gauge configurations into the three Polyakov loop sectors. Table 2 shows the corresponding results for our β G = 8.2 ensemble. It is obvious from the table, that for the complex Polyakov sectors the average size of the grand canonical determinant is much smaller than for the real Polyakov loop sector. Thus the latter sector receives a much larger weight in the path integral and thus is selected by the system in the high temperature phase. When analyzing | det[D(µ)]| G in the low temperature phase, we found essentially no discrepancy between the three Polyakov loop sectors 2 . We conclude that in the high temperature phase the selection of the real Polyakov loop sector in the dynamical case can be understood as a consequence of the center symmetry properties of the canonical determinants.

Center symmetry and observables
Having demonstrated the importance of the center symmetry, we briefly discuss its relevance for observables. An observable is a gauge invariant object and as such may be expressed as a sum over products of closed loops dressed with link variables 3 , similar to Eq. (17) which gives the loop representation for the fermion determinant. Each term in this sum may be classified with respect to the total winding number r ∈ Z of its loops. Organizing the loop representation of an observable O with respect to the winding numbers, we find Since the terms O (r) consist of loops with a total winding number of r, they transform under the center transformation (19) as Let us now analyze the vacuum expectation value of an observable. We find In the second step we inserted the fugacity expansion (11) for the determinant and the sum (39) for the observable. The individual terms in the sum transform under a center transformation as As long as the center symmetry is unbroken, this implies that the contribution vanishes unless (q + r) mod 3 = 0 , which follows from the same reasoning as used in (22). In other words, only those terms, where the combined triality of the canonical determinant and the observable term O (r) is trivial, give a non-vanishing contribution in vacuum expectation values. For high temperature, where the center symmetry becomes broken, the three center elements are not summed over in the path integral and the selection rule (42) does no longer apply, and O (r) receives contributions from all quark sectors in the center broken phase. Particularly simple are observables where the sum (39) reduces to a single term with a fixed value of r. Examples are the Polyakov loop (r = 1), the Wilson loop (r = 0), the plaquette (r = 0), the dual chiral condensate [11] (r = 1) and similar projected observables [12] (r = 1). The whole observable in the low temperature phase then receives only contributions from the quark sectors q that have (q + r) mod 3 = 0. For example the Polyakov loop has its leading contribution in the low temperature phase from the sector q = −1, while in the high temperature phase the leading contribution comes from q = 0 (see also [4,13] for this example).

Size distribution properties of canonical determinants
Analyzing scatter plots in Section 5.3 we established interesting properties for the angular distribution of the canonical determinants at high temperature, which, as we illustrated, is related to center symmetry. In this section we will show that also the modulus of the canonical determinants shows an interesting behavior.
We begin our discussion with Fig. 8, where we plot the average over the modulus of the canonical one flavor determinants | det 1 [D] (q) / det 1 [D] (0) | as a function of q, normalized relative to the q = 0 determinant. In other words we study the average size of the canonical determinant for different q. In the top panel of Fig. 8 we show for the 8 3 × 4 ensembles the results for the low temperature phase, i.e., from β G = 7.4 up to β G = 7.8 (which is roughly the gauge coupling for the transition). The bottom panel is for high temperature (β G = 7.8, ... 8.2).
The low temperature results presented in the top plot of Fig. 8 show a similar behavior as the results for the free case (compare the µ = 0 data in Fig. 3). The distribution is symmetrical around q = 0 (as it must be since det 1 [D] (−q) = (det 1 [D] (q) ) * ) and has a Gaussian-like distribution as a function of q. The distribution is rather narrow for the lowest temperature (β G = 7.4) and widens as β G and thus the temperature are increased.
As a consistency check we summed up the fugacity expansion (11) with the available terms, q ∈ [−6, 6], and compared it to the exact result for the grand canonical determinant. For the lowest values of β G , where, as Fig. 8 shows, the canonical determinants become small very quickly, we found excellent agreement, and the fugacity expansion with q ∈ [−6, 6] reproduces the grand canonical µ = 0 determinant with an error smaller than 10 −6 (on 8 3 ×4). When  increasing µ, which shifts the distribution towards larger q (compare Fig. 3), we found that the error remains below 1% up to µ = 0.35 (in lattice units). This finding is an interesting result which shows the fast convergence of the fugacity expansion at low temperatures and moderate lattice volumes. Beyond that it is also an important consistency and accuracy check for our evaluation of the canonical determinants. We stress that the achieved accuracy was not only established for the ensemble average, but holds also for individual configurations.
Let us now come to the bottom plot in Fig. 8, which shows our results in the high temperature phase. It is obvious that the trend towards a wider distribution at increasing temperature (increasing β G ) continues, although for the largest values of β G the widening slows down. Most interesting, however, is the observation that the shape of the distribution is no longer Gaussian, but a relative enhancement of the determinant values for even q becomes visible for the highest values of β G . In order to study this enhancement of even quark numbers further, in Fig. 9 we compare the distribution for the low and high temperature phases now using the two ensembles on the larger 10 3 × 4 lattices. While the canonical determinants in the low temperature phase (β G = 7.5, T /T c ∼ 0.77) again show a Gaussian-like distribution, for the high temperature phase (β G = 8.1, T /T c ∼ 1.32) the enhancement of the even quark numbers is much more pronounced for the larger 10 3 × 4 volumes shown in this figure. On the other hand, we found that on our smaller 6 3 × 4 volumes the effect vanishes. Thus the enhancement of even quark numbers is an effect that needs a minimal spatial volume, and probably would become more pronounced if the volume was increased further.
The relative enhancement of the canonical determinants for even quark numbers may be interpreted as an enhancement of meson states, of diquark states, et cetera in the high temperature phase, while single quark states, baryon states, et cetera are suppressed 4 . A possible formation of diquark correlations at not too large temperatures has been suggested [14] and preliminary lattice studies may be found in the literature [15]. Although it certainly needs to be improved (larger volumes, more systematical scanning of temperatures and masses), our study provides further ab-initio support for such a picture and generalizes it to all even quark number sectors. Of course a crucial question, which will have to be addressed in future work, is whether the effect persists when a better sampling beyond the pure gauge ensembles is being used. Such an analysis, as well as the question what happens at larger values of T /T c (here we have T /T c ∼ 1.32) is planned for the future.
We stress again at this point, that in the low temperature phase also a suppression of certain quark sectors is manifest -the sectors with non-trivial triality. At low temperature, however, the mechanism is a different one, as we discussed in Section 2.4: The unbroken center symmetry leads to the survival of only the sectors with vanishing triality, i.e., mesons, baryons et cetera. This happens through the averaging over the phases of the canonical determinants. Since in our observable | det 1 [D] (q) / det 1 [D] (0) | the phases are neglected, the suppression of the sectors with non-vanishing triality is not manifest in the corresponding low temperature data in Figures 8 and 9.

Summary and outlook
In this article we discuss the numerical evaluation and properties of canonical fermion determinants using Wilson fermions on pure SU(3) configurations. We obtain a considerable speedup in the evaluation of the canonical determinants by implementing a domain decomposition approach which leads to a dimensional reduction of the fermion determinant [6].
Based on this dimensional reduction we compare three different strategies for evaluating the canonical determinants: The conventional Fourier transformation approach, a generalized hopping expansion, and a generalized Fourier method, which combines a Fourier transform of the logarithm of the determinant with expansion techniques. Our comparison shows that all three methods work, are sufficiently accurate for the volumes, quark numbers q and couplings we considered, and allow for a systematical improvement of the accuracy if needed.
We evaluate the canonical determinants with quark numbers q ∈ [−6, 6], which for our lower temperatures and not too large chemical potential is sufficient for a very accurate representation of the grand canonical determinant through the corresponding fugacity expansion. The obtained accuracy of the fugacity expansion is an important test for the correctness of our implementation. For larger chemical potential and high temperatures more values of the quark number q would be needed.
Subsequently we explore the properties of the canonical determinants. Scatter plots of the determinant values show that in the low temperature phase the determinant values are distributed essentially isotropically in the complex plane. This holds also for the sectors with vanishing triality (except for the zero quark sector), illustrating that the evaluation of the free energy of, e.g., a baryon from the corresponding q = 3 canonical determinant is heavily plagued by phase cancellations.
In the high temperature phase we observe a very pronounced center pattern of the canonical determinants and we show that this pattern behaves exactly as expected from the center transformation properties of the canonical determinants.
Based on the observation that the center symmetry is so strongly manifest in the canonical determinants, we explain and numerically confirm that the grand canonical determinant is much smaller for high temperature gauge configurations with complex phases of the Polyakov loop. Thus the fact that the Polyakov loop comes out real in dynamical simulations may be understood as a consequence of center symmetry and its breaking. The discussion of the role of the center symmetry is furthermore extended to observables, giving rise to selection rules for the relevant quark sectors in the low temperature phase.
Analyzing the distribution of the canonical determinants as a function of the quark number, we find for low temperatures a Gaussian-like distribution which widens as the temperature is increased. At high temperatures an interesting effect is observed: Canonical fermion determinants with even quark numbers q are enhanced in the high temperature phase. This might indicate an enhancement of diquark states et cetera.
We conclude with remarking again that an important caveat must be kept in mind: The numerical results were obtained using ensembles of pure gauge theory. Although this is formally correct, in a numerical simulation it could be the case that the distribution of the generated gauge configurations has only a small overlap with the main contributions to the canonical determinants, and very large ensembles might be needed for accurate results. This is an open issue which might be answered by using dynamical ensembles and suitable reweighting -a question that we plan to address in the future.