Oct 16, 2008 - phase transition line in the (T,Âµq) plane is interesting both from the ..... 2: Distribution of the saddle point at Î² = 3.55, 3.63, 3...

0 downloads 13 Views 684KB Size

Canonical partition function and finite density phase transition in lattice QCD Shinji Ejiri Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA (Dated: October 10, 2008)

arXiv:0804.3227v3 [hep-lat] 16 Oct 2008

We discuss the nature of the phase transition for lattice QCD at finite temperature and density. We propose a method to calculate the canonical partition function by fixing the total quark number introducing approximations allowed in the low density region. An effective potential as a function of the quark number density is discussed from the canonical partition function. We analyze data obtained in a simulation of two-flavor QCD using p4-improved staggered quarks with bare quark mass m/T = 0.4 on a 163 × 4 lattice. The results suggest that the finite density phase transition at low temperature is of first order. PACS numbers: 11.15.Ha, 12.38.Gc, 12.38.Mh

I.

INTRODUCTION

The study of the QCD phase diagram at nonzero temperature (T ) and quark chemical potential (µq ) is one of the most important topics among studies of lattice QCD. In particular, the study of the endpoint of the first order phase transition line in the (T, µq ) plane is interesting both from the experimental and theoretical point of view. The existence of such a critical point is suggested by phenomenological studies [1, 2, 3]. The appearance of the critical point in the (T, µq ) plane is closely related to hadronic fluctuations in heavy ion collisions and may be experimentally examined by an event-by-event analysis of heavy ion collisions. Many trials have been made to find the critical point by first principle calculation in lattice QCD [4, 5, 6, 7, 8, 9, 10, 11, 12]. However, no definite conclusion on this issue is obtained so far. One of the interesting approaches is to construct the canonical partition function ZC (T, N ) by fixing the total quark number (N ) or quark number density (ρ) [13, 14, 15, 16, 17, 18]. The canonical partition function is obtained from the grand canonical partition function ZGC (T, µq ) by an inverse Laplace transformation. The relation between ZGC (T, µq ) and ZC (T, N ) is given by Z X N ZGC (T, µq ) = DU (det M (µq /T )) f e−Sg = ZC (T, N )eN µq /T , (1) N

where det M is the quark determinant, Sg is the gauge action, and Nf is the number of flavors. Nf in this equation must be replaced by Nf /4 when one uses a staggered type quark action. In order to investigate the net quark number giving the largest contribution to the grand canonical partition function at (T, µq ), it is worth introducing an effective potential Veff as a function of N , Veff (N ) ≡ − ln ZC (T, N ) − N

µq f (T, N ) µq = −N , T T T

(2)

where f is the Helmholtz free energy. If there is a first order phase transition region, we expect this effective potential has minima at more than one value of N . At the minima, the derivative of Veff satisfies ∂(ln ZC ) µq ∂Veff (N, T, µq ) = − (T, N ) − = 0. ∂N ∂N T

(3)

Hence, in the first order transition region of T , we expect ∂(ln ZC )/∂N (T, N ) ≡ −µ∗q /T takes the same value at different N . Here, µ∗q (T, N ) is the chemical potential which gives a minimum of the effective potential at (T, N ). The phase structure in the (T, ρ) plane and the expected behavior of µ∗q /T are sketched in the left and right panels of Fig. 1, respectively. The thick lines in the left figure are the phase transition line. We expect that the transition is crossover at low density and becomes of first order at high density. Since two states coexist on the first order transition line, the phase transition line splits into two lines in the high density region, and the two states are mixed in the region between two lines. The expected behavior of µ∗q along the lines A and B are shown in the right figure. When the temperature is higher than the temperature at the critical point Tpc (line A), µ∗q increases monotonically as the density increases. However, for the case below Tcp (line B), this line crosses the mixed state. Because the two states of ρ1 and ρ2 are realized at the same time, µ∗q does not increase in this region between ρ1 and ρ2 . The Glasgow method [14, 15] has been a well-known method to compute the canonical partition function. Recently, such a behavior at a first order phase transition has been observed by Kratochvila and de Forcrand in 4-flavor QCD

2

*

T Tcp

µq

hot phase

T

A critical point

B

B

A

cold phase mixed state ρ1 ρ2

ρ

ρ1/T

3

ρ2/T

3

ρ/T

3

FIG. 1: Phase structure in the (T, ρ) plane and the behavior of µ∗q /T as a function of ρ.

with staggered fermions on a 63 × 4 lattice [18] calculating the quark determinant by the Glasgow algorithm. However, with present day computer resources, the study by the Glasgow method is difficult except on a small lattice. Therefore, it is important to consider a method available for a simulation on a large lattice. In this paper, we propose such a method for the calculation of the canonical partition function introducing approximations allowed in the low density region. The method proposed in this paper is based on the following ideas. We adopt a saddle point approximation for the inverse Laplace transformation from ZGC to ZC . This approximation is valid when the volume size is sufficiently large. We moreover perform a Taylor expansion of ln det M (µq ) in terms of µq around µq = 0 and calculate the expansion coefficients, as proposed in [19]. The Taylor expansion coefficients are rather easy to calculate by using the random noise method. The saddle point approximation is also based on a Taylor expansion around the saddle point. We estimate the Taylor expansion coefficients at the saddle point from the Taylor expansion at µq = 0. Although we must cut off this expansion at an appropriate order in µq , this approximation is applicable when the saddle point is not far from µq = 0, and we can estimate the application range where the approximation is valid for each analysis. For the calculation of ZC two further technical problems must be solved. The first problem is the problem of importance sampling in Monte-Carlo simulations (overlap problem). Since the method proposed in this paper is a kind of reweighting method, the configurations which give important contributions will change when the weight factor is changed by µq [20, 21]. To avoid this problem, we combine configurations generated at many simulation points β = 6/g 2 covering a wide range of the temperature using a method introduced in [22]. The second problem is the sign problem. We must deal with an expectation value of a complex number in this method. If the fluctuation of the complex phase is large, the statistical error becomes larger than the mean value. We use a technique introduced in [12]. We consider the probability distribution function in terms of the complex phase of the complex operators. We assume the distribution function is well approximated by a Gaussian function and perform the integration over the phase. Once this assumption is adopted, the sign problem is completely solved. This assumption is reasonable for sufficiently large volume and small µq /T . During the process of this calculation, we will find out why the quark number density changes sharply at the transition point and why the density approaches the value of the free quark gas in the high density limit even at low temperature. The configurations generated in µq = 0 simulations at low temperature are gradually suppressed as the density increases. In the next section, we explain the method to calculate the canonical partition function using the inverse Laplace transformation of ZGC within a saddle point approximation. The problem of the Monte-Carlo sampling is discussed in Sec. III. The sign problem is discussed in Sec. IV. We evaluate ∂(ln ZC )/∂N using data obtained with two-flavors of p4-improved staggered quarks in [6]. The result is shown in Sec. V. The behavior of ∂(ln ZC )/∂N suggests that the phase transition is of first order in the low temperature and high density region. Conclusions are given in Sec. VI.

3 II.

CANONICAL PARTITION FUNCTION

We calculate the canonical partition function using Ns3 × Nt lattice and investigate the effective potential Veff (N ). From Eq. (1), the canonical partition function can be obtained by an inverse Laplace transformation [13, 16, 17, 18], 3 ZC (T, N ) = 2π

Z

π/3

−π/3

e−N (µ0 /T +iµI /T ) ZGC (T, µ0 + iµI ) d

µ I

T

,

(4)

where µ0 is an appropriate real constant and µI is a real variable. Note that ZGC (T, µq + 2πiT /3) = ZGC (T, µq ) [23]. The grand canonical partition function can be evaluated by the calculation of the following expectation value at µq = 0. * N N + Z ZGC (T, µq ) det M (µq /T ) f 1 det M (µq /T ) f Nf −Sg = (det M (0)) e = DU . (5) ZGC (T, 0) ZGC det M (0) det M (0) (T,µq =0)

However, with present day computer resources, the exact calculation is difficult except on small lattices. We consider an approximation which is valid for large volume and low density. If we select a saddle point as µ0 in Eq. (4) when the volume is sufficiently large, the information which is needed for the integral is only the value of det M around the saddle point. Furthermore, if we restrict ourselves to study the low density region, the value of det M (µq /T ) near the saddle point can be estimated by a Taylor expansion around µq = 0. The calculations by the Taylor expansion are much cheaper than the exact calculations and the studies using large lattices are possible. First, we perform the integral in Eq. (4) by a saddle point approximation. We denote the quark number density in a lattice unit and physical unit as ρ¯ = N/Ns3 and ρ/T 3 = ρ¯Nt3 , respectively. We assume that a saddle point z0 exists in the complex µq /T plane for each configuration, which satisfies [D′ (z) − ρ¯]z=z0 = 0,

(6)

where (det M (z)/ det M (0))Nf = exp[Ns3 D(z)] and D′ (z) = dD(z)/dz. We then perform a Taylor expansion around the saddle point and obtain the canonical partition function, *Z Nf + π/3 det M (z0 + ix) 3 −N (z0 +ix) dx ZGC (T, 0) e ZC (T, ρ¯V ) = 2π det M (0) −π/3 (T,µq =0) *Z + π/3 3 1 ′′ 2 = ZGC (T, 0) dx exp V D(z0 ) − ρ¯z0 − D (z0 )x + · · · 2π 2 −π/3 (T,µq =0) s + * 3 1 . ≈ √ ZGC (T, 0) exp [V (D(z0 ) − ρ¯z0 )] e−iα/2 V |D′′ (z0 )| 2π

(7)

(T,µq =0)

Here, D′′ (z) = d2 D(z)/dz 2 , V ≡ Ns3 and D′′ (z) = |D′′ (z)|eiα . We chose a path which passes the saddle point. Higher order terms in the expansion of D(z) becomes negligible when the volume V is sufficiently large, since the saddle point approximation is a 1/V expansion. Next, we calculate the quark determinant by the Taylor expansion around µq = 0. We define n ∂ ln det M (µq /T ) 1 . (8) Dn = n!Ns3 Nt ∂(µq /T )n µq =0 The even derivatives of ln det M are real and the odd derivatives are purely imaginary [19]. The calculation of Dn is rather easy using the stochastic noise method. D(z), D′ (z) and D′′ (z) in Eq. (6) and (7) can be evaluated by D(z) = Nf Nt

∞ X

n=1

Dn z n ,

D′ (z) = Nf Nt

∞ X

n=1

nDn z n−1 ,

D′′ (z) = Nf Nt

∞ X

n=2

n(n − 1)Dn z n−2 .

(9)

Because Im(D1 ) ≪ Re(D2 ), the saddle point, i.e. the solution of Eq. (6), is distributed near the real axis and Re(z0 ) increases as ρ increases. Moreover, for the case that the saddle point z0 is on the real axis, the saddle point condition is the same as the condition where Re(D(z) − ρ¯z) is minimized. Hence, exp[V Re(D(z0 ) − ρ¯z0 )] in Eq. (7) decreases exponentially as Re(z0 ) increases.

4 In this study, we want to focus on the derivative of the effective potential with respect to N or ρ. Since the effective potential Veff (N ) is minimized in the thermodynamic limit, i.e. ∂ log ZC /∂N + µq /T = 0, we denote the derivative by µ∗q ∂ ln ZC (T, N ) 1 ∂ ln ZC (T, ρ¯V ) =− = − . T ∂N V ∂ ρ¯ Within the framework of the saddle point approximation, this quantity can be evaluated by D E q z0 exp [V (D(z0 ) − ρ¯z0 )] e−iα/2 V |D′′1(z0 )| ∗ µq (T,µq =0) E q . ≈ D 1 T −iα/2 exp [V (D(z0 ) − ρ¯z0 )] e V |D′′ (z0 )|

(10)

(11)

(T,µq =0)

This equation is similar to the formula of the reweighting method for finite µq . The operator in the denominator corresponds to a reweighting factor, and µ∗q /T is an expectation value of the saddle point calculated with this modification factor. We denote the real and imaginary parts of the logarithm of the weight factor by F and θ, ! ∞ X 1 iα F + iθ ≡ V Nf Nt Dn z0n − ρ¯z0 − ln[V |D′′ (z0 )|] − . (12) 2 2 n=1 We define the complex phase of the weight factor by θ and the absolute value of the reweighting factor is exp(F ). This weight factor plays an important role around the phase transition point at finite density. In the calculation to derive Eq. (7), we replaced the order of the path integral of gauge fields and the integral for the inverse Laplace transformation. This replacement is essentially important. If one calculates ZC from ZGC using Eq. (4) after the path integral, the equation which is satisfied at a saddle point is N = ρ¯V = ∂(ln ZGC )/∂(µq /T ).

(13)

Hence, in the thermodynamic limit, i.e. when we ignore the finite volume correction, µ∗q is just equal to the inverse function of ρ(µq ) at the saddle point. Therefore, ρ(µq ) must be a discontinuous function or a multivalued function at a first order phase transition to obtain the behavior of µ∗q (ρ) shown in Fig. 1. However, if we calculate ln ZGC by a Taylor expansion in µq at a temperature in the hadron phase, ρ(µq ) cannot be a discontinuous function. In this study, we use Eq. (11). As we will discuss in detail, we can obtain the behavior of µ∗q /T suggesting a first order phase transition, although the calculations of z0 and the weight factor in Eq. (11) are based on the Taylor expansion. The important point is that the weight factor exp(F + iθ) gives the same effect when the temperature changed and configurations which give important contributions to the calculation of µ∗q /T change gradually as ρ increases. Hence µ∗q /T does not need to increase monotonously as a function of ρ even if the saddle point z0 is a monotonous function of ρ for each configuration. III.

MONTE-CARLO ANALYSIS AND REWEIGHTING METHOD FOR β-DIRECTION

We analyze data obtained in a simulation of two-flavor QCD using p4-improved staggered quarks with bare quark mass m/T = 0.4 [6]. These data are obtained at 16 simulation points from β = 3.52 to 4.00. The corresponding temperature normalized by the pseudocritical temperature is in the range of T /Tc = 0.76 to 1.98, and the pseudocritical point (T /Tc = 1) is βpc ≈ 3.65. The ratio of pseudoscalar and vector meson masses is mPS /mV ≈ 0.7 at β = 3.65. The lattice size Nsite = Ns3 × Nt is 163 × 4. The number of configurations Nconf. is 1000 – 4000 for each β. Further details on the simulation parameters are given in [6]. We use the data of the Taylor expansion coefficients Dn up to O(µ6q ). The saddle point is found by the following procedure. (1) Because absolute values of the odd terms of Dn are much P smaller than the even terms, we first find a solution when the odd terms are neglected, i.e. a solution of Nf Nt n 2nD2n z 2n−1 − ρ¯ = 0. The odd terms of Dn are purely imaginary and the even terms are real. Although some fake solutions may appear at large z due to the truncation of the higher order terms, the solutionP is found on the real axis in the low density regime. (2) Next, in the vicinity of this solution, we calculate r2 ≡ |Nf Nt n nDn z n−1 − ρ¯|2 and find the point where r2 is zero. For the data we used in this analysis, the saddle point could be found for every configuration by this procedure except in the very low density region of ρ/T 3 < 0.37. Figure 2 shows an example of the distribution of the saddle points, obtained at β = 3.55, 3.63, 3.70 with ρ/T 3 = 2.0. For the calculation of the derivative of ln ZC , the application of the reweighting method for β-direction is crucial. Configurations in a Monte-Carlo simulation are generated with the probability in proportion to the product of the

5 1 3

ρ/T =2.0

Im[z0]

0.5

0

β=3.70

β=3.63 β=3.55

-0.5

-1

0

1

2

3

4

Re[z0] FIG. 2: Distribution of the saddle point at β = 3.55, 3.63, 3.70 with ρ/T 3 = 2.0. β=3.52 T/Tc=0.76

3.65 1.00

3.80 1.36

4.00 1.98

W(P,β)

3.60 0.90

0.78

0.8

0.82

0.84

0.86

0.88

0.9

0.92

0.94

0.78

0.8

0.82

0.84

0.86

0.88

0.9

0.92

0.94

-ln W(P,β)

4

2

0

P FIG. 3: Plaquette histogram w(P, β) at µq = 0.

weight factor (det M )Nf e−Sg and the state density of the link fields {Uµ (x)}. The expectation value is then estimated by taking an average of the operator O[Uµ ] over the generated configurations {Uµ (x)}. hOi(β) ≈

1 Nconf.

X

{Uµ (x)}

O[Uµ ].

(14)

However, if the value of O[Uµ ] changes very much during a Monte-Carlo simulation and the change of O[Uµ ] is much larger than the size of the probability p(O, β), the Monte-Carlo method is no longer valid. For example, configurations which have large O × p(O, β) are important for the evaluation of hOi(β) , but such configurations are not generated if p(O, β) is too small. Such a problem occurs in the calculation of Eq. (11). This problem is called an “overlap problem.” To clarify the problem of the importance sampling, we rewrite Eq. (11) as R µ∗q hz0 exp[F + iθ]iP w(P, β)dP = R . (15) T hexp[F + iθ]iP w(P, β)dP P is the plaquette value, h· · · iP denotes the expectation value for fixed P at µq = 0, and w(P, β) is the probability distribution of P at β, Z ′ w(P , β) = DU δ(P − P ′ )(det M (0))Nf e−Sg (β) ∝ hδ(P − P ′ )i(β,µq =0) . (16)

6 6

3

ρ/T =10.0 3 ρ/T =9.0 3 ρ/T =8.0 3 ρ/T =7.0 3 ρ/T =6.0 3 ρ/T =5.0 3 ρ/T =4.0 3 ρ/T =3.0 3 ρ/T =2.0 3 ρ/T =1.0

<|z0|>P 5

4

3

2

1

0

0.8

0.9

0.85

P FIG. 4: Expectation value of |z0 | with fixed P for each ρ/T 3 .

0

-500

-1000 3

ρ/T =1.0 3 ρ/T =2.0 3 ρ/T =3.0 3 ρ/T =4.0 3 ρ/T =5.0 3 ρ/T =6.0 3 ρ/T =7.0 3 ρ/T =8.0 3 ρ/T =9.0 3 ρ/T =10.0

-1500

-2000

-2500 0.8

0.9

0.85

P FIG. 5: Expectation vale of F with fixed P for each ρ/T 3 .

This kind of analysis is called the density of state method [24, 25, 26, 27, 28, 29]. We define the average plaquette as P = −Sg /(6Nsite β) for later discussions. This P is the plaquette value for the standard gauge action, but is a linear combination of Wilson loops for improved gauge actions. Because R hXδ(P − P ′ )i(β,µq =0) DU Xδ(P − P ′ )(det M (0))Nf , (17) = R hXiP ′ ≡ ′ hδ(P − P )i(β,µq =0) DU δ(P − P ′ )(det M (0))Nf

hXiP is independent of β for an operator X which does not depend on β explicitly. Hence, hXiP can be computed at an appropriate β. The probability distribution functions w(P ) and − ln w(P ) are given in Fig. 3. We show the values of β and corresponding T /Tc above these figures. To obtain w(P ), we grouped the configurations by the value of P into blocks and counted the number of configurations in these blocks. − ln w(P ) is normalized by the minimum value for each β. Because the transition from the hadron phase to the quark-gluon phase is a crossover for two-flavor QCD, the distribution is always of Gaussian type, and the width of the distribution becomes narrower as the volume increases. Moreover, since the suppression factor is exp(6Nsite βP ), the peak position of the distribution w(P ) moves to the right as β increases. We also calculate the expectation value of |z0 | and F when P is fixed, h|z0 |iP ′ = h|z0 |δ(P − P ′ )i/ hδ(P − P ′ )i,

hF iP ′ = hF δ(P − P ′ )i/ hδ(P − P ′ )i.

(18)

The result of h|z0 |iP is plotted in Fig. 4, and solid lines in Fig. 5 are hF iP for each ρ/T 3 . (Dashed lines will be explained in the next section.) For the calculation of these quantities, we use the delta function approximated by a Gaussian

7 20000

−ln[ZGC(β)/ZGC(3.65)] 10000

0

-10000

-20000

-30000

-40000 3.5

3.7

3.6

β

3.8

3.9

4

FIG. 6: β dependence of − ln ZGC determined by the consistency condition Eq. (21).

√ function, δ(x) ≈ 1/(∆ π) exp[−(x/∆)2 ]. If making ∆ small, this approximation becomes better but statistical errors become larger. Hence, the size of ∆ must be adjusted appropriately. In this study, we adopt ∆ = 0.0025. Let us now consider hexp[F + iθ]iP × w(P ) in Eq. (15). Since hF iP increases linearly for P < ∼ 0.85, Fig. 5 suggests that hexp[F + iθ]iP increases exponentially as P increases in this range. Moreover, the slope of hF iP is increasing as ρ/T 3 increases. Therefore, for large ρ/T 3 , this hexp[F + iθ]iP × w(P ) may not decrease even if w(P ) decreases exponentially at the tail of the distribution generated by a simulation. If the value of hexp[F + iθ]iP × w(P ) is still large even in the region where the configurations are not generated, the calculation by the Monte-Carlo method is completely wrong. Because the width of w(P ) becomes narrow for large V , it is essentially important to solve this problem if we want to use a large lattice which is required for the saddle point approximation. To avoid this problem, we combine all configurations obtained in simulations with many different β values, using the Ferrenberg and Swendsen method [22]. The expectation value hOiβ can be calculated from the data obtained by more than one simulation point, βi (i = 1, 2, · · · , Nβ ) by the following equation: hOiβ ≈

hOG(β, P )iall . hG(β, P )iall

(19)

Here, the weight factor G(β, P ) is G(β, P ) = PNβ

i=1

e6Nsite βP −1 Ni e6Nsite βi P ZGC (βi )

,

(20)

where Ni is the number of configurations at simulation points βi and h· · · iall means the average over all configurations generated at all βi . The derivation of this equation is given in Appendix A. The partition function ZGC (βi ) is determined by a consistency condition for each i, ZGC (βi ) = hG(βi , P )iall

(21)

This equation can be solved except for the normalization factor. The result of − ln[ZGC (β)/ZGC (3.65)] is plotted in Fig. 6. We should note that G(β, P ) is independent of the simulation points βi at which the operators are measured and the expectation value is simply given by the average over all configurations generated at many β. If we perform simulations at many different β and combine the data, the configurations are distributed in a wide range of P . (See Fig. 3.) Among these configurations, important configurations for each calculation are selected by O and G(β, P ) automatically. This method is particularly important when the volume is large, since the distribution w(P ) is narrow if we generate configurations on a large lattice with single β. The overlap problem is solved by this method. However, the statistical error is enlarged by the fluctuations of exp[F + iθ] when the density is increased, hence the application range of ρ is determined by the statistical error. Also, we should check that the important configurations are within the range of the plaquette distribution for each calculation. We will discuss this point in Sec. V again.

8

β=3.55

-80

-60

-40

-20

0

θ

20

40

60

β=3.63

80

-40

-20

0

θ

20

β=3.70

40

-10

-5

0

θ

5

10

FIG. 7: Histograms of the complex phase θ for ρ/T 3 = 2.0 at β = 3.55 (left), 3.63 (middle) and 3.70 (right). Dashed lines are the fit results by Gaussian functions.

IV.

SIGN PROBLEM

Next, we discuss the sign problem that also shows up in the calculation of the canonical partition function. Because z0 and the additional weight factor in Eq. (11) are complex numbers, the calculation of Eq. (11) suffers from the sign problem [21, 30]. If the weight factor changes the sign frequently, both the numerator and denominator of Eq. (11) become smaller than their statistical errors. To avoid the sign problem, we use a method proposed in [12]. As in Eq. (12), we define the complex phase of the weight factor by !# " ∞ X α n − . (22) θ ≡ Im V Nf Nt Dn z0 − ρ¯z0 2 n=1 In this definition, θ is not restricted to the range from −π to π because there is no reason that the imaginary part of Eq. (12) must be in the finite range. In fact, this quantity becomes larger as the volume increases. It has been discussed in [12] that histograms of Dn are well approximated by Gaussian functions if a simulation is performed at a point away from the critical point with sufficiently large volume. The Taylor expansion coefficients in Eq. (12) are given by combinations of traces of products of ∂ n M/∂(µq /T )n and M −1 . (See the appendix of [6].) Therefore, Dn are obtained by the sum of the diagonal elements of such matrices. When the correlation among the diagonal elements is small and the volume is sufficiently large, the distribution functions of the expansion coefficients and θ should be of Gaussian type due to the central limit theorem. For example, the diagonal element of the first coefficient, Im[∂(ln det M )/∂(µq /T )] = Im[tr[M −1 (∂M/∂(µq /T ))]], is the imaginary part of the local number density operator at µq = 0. If the spatial density correlation is not very strong, a Gaussian distribution is expected. Figure 7 is the histogram of the complex phase θ at ρ/T 3 = 2.0 for the two-flavor QCD simulations with p4-improved staggered quarks at β = 3.55 (left), 3.63 (middle), and 3.70 (right). These figures suggest that the distribution of θ is well approximated by a Gaussian function. We fit the data of the histogram by a Gaussian function. Dashed lines are the fit results. The width of the Gaussian function is different for each distribution obtained by a different parameter. If we restrict the phase to the range from −π to π by subtracting 2πn (n: integer), the complex phase distribution is almost flat for the case that the width is much larger than π, and the flatness indicates the seriousness of the sign problem. However, the measurement of the width of the Gaussian distribution is easier than the estimation of the flatness of the restricted phase distribution. Once we assume a Gaussian distribution for θ, the problem of complex weights can be avoided. We introduce the probability distribution w ¯ as a function of the plaquette P , F , and θ, Z ′ ′ ′ w(P ¯ , F , θ ) ≡ DU δ(P ′ − P )δ(F ′ − F )δ(θ′ − θ)(det M (0))Nf e6βNsite P . (23) The distribution function itself is defined as an expectation value at µq = 0, however F and θ are functions of ρ. The denominator of Eq. (11), hexp[F + iθ]i, is given by Z Z Z

F iθ 1 e e (T,µq =0) = dP dF dθ eF eiθ w(P, ¯ F, θ), (24) ZGC R R R where ZGC = dP dF dθ w(P, ¯ F, θ). Because we calculate this expectation value by the reweighting method using Eq. (19), the operator in the calculation of Eq. (24) is a function of P , F and θ.

9 Since the partition function is real even at nonzero density, the distribution function is symmetric under the change from θ to −θ. Therefore, the distribution function is a function of θ2 , e.g., w(θ) ¯ ∼ exp[−(a2 θ2 + a4 θ4 + a6 θ6 + · · · )]. And, the distribution function is expected to be well approximated by a Gaussian function when the system size is sufficiently large in comparison to the correlation length. We assume the following distribution function in terms of θ when P and F are fixed: r a2 (P, F ) ′ w(P, ¯ F, θ) ≈ w ¯ (P, F ) exp −a2 (P, F )θ2 . (25) π The coefficient a2 (P, F ) is given by

Z dθ θ2 w(P ¯ ′ , F ′ , θ) dθ w(P ¯ ′ , F ′ , θ)

2 θ δ(P ′ − P )δ(F ′ − F ) (T,µq =0)

≡ θ2 (P,F ) . = ′ ′ hδ(P − P )δ(F − F )i(T,µq =0)

1 = 2a2 (P ′ , F ′ )

Z

The integration over θ can be carried out easily and we obtain the denominator of Eq. (11), r Z Z Z

F iθ 2 a2 ′ 1 e e (T,µq =0) ≈ w ¯ (P, F )e−a2 θ eF eiθ dP dF dθ ZGC π Z Z 1 ′ = dP dF w ¯ (P, F )eF e−1/(4a2 ) ZGC Z 1 DU eF e−1/(4a2 (P,F )) (det M (0))Nf e−Sg = ZGC E D . = eF e−1/(4a2 (P,F )) (T,µq =0)

(26)

(27)

Since θ is roughly proportional to the size of the quark matrix M , the value of 1/a2 becomes larger as the volume increases. Therefore, the phase factor in this quantity decreases exponentially as a function of the volume. However, in this framework, the operator in Eq. (27) is always real and positive for each configuration. This means that the expectation value is always larger than its statistical error. Therefore, the sign problem is completely avoided if we can assume the Gaussian distribution of θ. The effect from a non-Gaussian term (a4 ) is discussed in Appendix B for the case that a4 /a2 is small. In this case, the effect from a4 is at most O(µ6q ), hence the non-Gaussian term may be neglected in the low density region even if a4 is nonzero. The complex phase distribution of the quark determinant by chiral perturbation theory has been discussed in [31]. The Gaussian distribution is suggested in the 1-loop calculation. The complex phase factor for the calculation of hz0 exp[F + iθ]i in Eq. (11) is calculated repeating the same procedure. We consider the complex phase of z0 ≡ |z0 | exp[iθ′ ], where −π < θ′ ≤ π. Replacing θ with θ + θ′ in Eq. (27), the suppression factor exp[−h(θ + θ′ )2 i/2] is estimated. We plot the result of the average of θ2 as a function of P in Fig. 8, i.e. hθ2 iP ′ ≡ hθ2 δ(P − P ′ )i/hδ(P − P ′ )i.

(28)

We adopt the approximated delta function δ(P ) used when we computed h|z0 |iP and hF iP in Sec. III. Although we need to average θ2 as a function of P and F for the calculation of Eq. (27), F dependence is not considered in Fig. 8. It is found from this figure that hθ2 iP decreases linearly as P increases in the range of P < 0.85 and the phase fluctuation is small for P > 0.85. This means that the phase factor decreases exponentially as P decreases for P < 0.85. This behavior is similar to that of hF iP in Fig. 5. Therefore, the weight factor exp[F − hθ2 i(P,F ) /2] suppresses the contribution from configurations having small P for large ρ/T 3 . Moreover, this argument implies that configurations on which the sign problem is serious do not contribute to the actual calculations of expectation values. The reason is that the fluctuations of the complex phase hθ2 i(P,F ) are large on such configurations and the configurations are suppressed by the weight factor exp[−hθ2 i(P,F ) /2]. Therefore, even if the error due to the Gaussian approximation of the complex phase distribution becomes visible when the phase fluctuations are large, the error does not affect to the practical calculations of expectation values so much. Here, we should notice that the values of P and p F are strongly correlated. We estimate the width of the distribution of F for each P and ρ/T 3 by calculating ∆F ≡ h(F − hF iP )2 iP . Dashed lines above and below the solid line for hF iP in Fig. 5 are the values of hF iP ± ∆F . Most of the configurations characterized by P and F are distributed in the narrow region between these two dashed lines. Outside this bound an accurate calculation of hθ2 i(P,F ) is difficult, since the number of configurations is not enough for the average. However, if we consider that F is approximately given as a function of P on each configuration, hθ2 i(P,F (P )) will be a function of only P .

10 5000 2

<θ >P 4000

3000

2000

3

ρ/T =10.0 3 ρ/T =9.0 3 ρ/T =8.0 3 ρ/T =7.0 3 ρ/T =6.0 3 ρ/T =5.0 3 ρ/T =4.0 3 ρ/T =3.0 3 ρ/T =2.0 3 ρ/T =1.0

1000

0

0.8

0.9

0.85

P FIG. 8: Average of the square of the complex phase as a function of P and ρ/T 3 .

0.9

0.88 0.86 0.84 3

ρq/T =6.0

0.82

3

ρq/T =4.0

0.8

3

ρq/T =2.0 3

ρq/T =0.0

0.78 3.5

3.6

β

3.7

3.8

FIG. 9: Expectation value of the plaquette for each ρ/T 3 . Dashed lines are the peak positions of the plaquette distributions at µq = 0.

V.

RESULTS AND DISCUSSIONS

We calculate the slope of ln ZC , i.e. µ∗q /T , using Eq. (11). This quantity is given by the average of the saddle point z0 multiplying the additional weight factor exp[F + iθ]. The volume V = 163 is sufficiently large, and we assume that the complex phase distribution of the reweighting factor is well approximated by a Gaussian function as discussed in Sec IV. We then replace the weight factor by exp[F − hθ2 i(F,P ) /2] to eliminate the sign problem. We find a saddle point z0 numerically for each configuration, assuming z0 exists near the real axis in the low density region of the complex µq /T plane. Before showing the result for µ∗q /T , it is worth discussing the effect of the weight factor using Eq. (15). The weight factor can be approximately estimated by hexp[F − iθ]iP w(P, β) ≈ Ω(P ) exp 6Nsite βP + hF iP − hθ2 iP /2 (29)

because lnhexp[F + iθ]iP ≈ hF iP − hθ2 iP /2 in the leading order and w(P, β) can be written as Ω(P ) exp[6βNsite P ] from Eq. (20), where Ω(P ) is the state density in terms of P . Since the behavior of hF iP −hθ2 iP /2 for P < ∼ 0.85 is different, we consider these two regions separately. ∼ 0.85 and P > The configurations below and above P ∼ 0.85 are generated in the low and high temperature phases, respectively. In 2 the region P > ∼ 0.85, the P dependence of hF iP − hθ iP /2 is small. Therefore, the balance of the weight does not < change. On the other hand, for P ∼ 0.85, hF iP − hθ2 iP /2 increases linearly. This has the same effect as when β

11

µ*q/T

4

3 T/Tc=0.76 T/Tc=0.80 T/Tc=0.83 T/Tc=0.87 T/Tc=0.90 T/Tc=0.94 T/Tc=0.98 T/Tc=1.02 T/Tc=1.07 T/Tc=1.11 T/Tc=1.16

2

1 free gas limit 0

0

2

4

8

6

10

12

3

ρ/T

FIG. 10: Derivative of ln ZC as a function of the quark number density.

changes to βeff ≡ β +

dhF iP 1 dhθ2 iP − dP 2 dP

1 , 6Nsite

(30)

and the derivatives of hF iP and −hθ2 iP increase as ρ increases. This means that the peak position of the probability distribution of P , shown in Fig. 3 for µq = 0, moves to the right as ρ increases, according to the change of the effective β. This behavior is consistent with our usual expectation, i.e. a phase transition arises when the density is increased as well as with increasing temperature. For the case that important configurations change, the multiparameter (β) reweighting in Sec. III is effective, since the important configurations are automatically selected among all configurations generated at multi-β, and also this method is useful for the interpolation between the simulation points. We combine all data obtained at 16 points of β using the multi-β reweighting method. In order to observe the change of the important configurations, we calculate the expectation value of the plaquette with the additional weight factor in Fig. 9, hP i (T, ρ) ≈

hP exp [F + iθ]i(T,µq =0) hexp [F + iθ]i(T,µq =0)

.

(31)

The circles in Fig. 9 are hP i at ρ/T 3 = 0 computed without the multi-β reweighting method. The solid lines which connect these circles are the interpolation using the multi-β reweighting method. Of course, hP i for each simulation point and that obtained by the multi-β reweighting method are consistent with each other at the simulation points. However, the line of the interpolation shows waves at some places between the circles. Such a wave appears where the configurations are missing in Fig. 3. For the calculation at finite ρ/T 3 , we calculate the phase factor eiθ from hθ2 iP in Fig. 8. The F dependence of 2 hθ i(P,F ) is neglected because the values of F is approximately given as a function of P as seen in Fig. 5. The results for hP i at ρ/T 3 = 2.0, 4.0 and 6.0 are plotted from below. This quantity indicates the plaquette value of the most important configuration when the weight factor is modified. From this figure, we find that hP i becomes larger for β < βpc = 3.65 and does not change very much for β > 3.65 when ρ/T 3 is increased. This means that, even when we measure an expectation value at small temperature (small β), the configurations generated by simulations at higher temperature (larger P ) are used for the measurement at finite ρ/T 3 . Moreover, for the measurement at sufficiently large ρ/T 3 , the configurations generated in the low temperature phase are completely suppressed by the additional weight factor even when the temperature is small. This property explains why the phase transition happens when the density is increased while keeping the temperature fixed. However, the errors on these results become larger as ρ/T 3 increases. One of the reasons may be the missing configurations between the peaks of the plaquette distributions. As seen in Fig. 3, the configurations are not distributed uniformly in the range of P which is necessary in this analysis, and correct results cannot be obtained if the important configurations are missing. At low temperature, the important value of P changes very much as ρ increases, therefore

12 we calculate µ∗q /T only when the expectation value of P is at the peak positions of the plaquette distributions in Fig. 3. Dashed lines in Fig 9 are the peak positions. Another important point is that h|z0 |iP and hF iP are strongly correlated with each other. Figure 4 is the average of |z0 | as a function of the plaquette value of each configuration. |z0 | increases as P decreases. Because D1 is purely imaginary, −F = Re[V (¯ ρz0 − Nf Nt D1 z0 )] + O(z02 ) becomes large as Re(z0 ) increases, which is seen in Fig. 5, and the contribution from the configurations which have large z0 is suppressed by the additional weight factor. Although the value of |z0 | for each configuration increases monotonically as a function of ρ/T 3 , nontrivial behavior in µ∗q /T is expected due to the suppression factor. We plot the result of µ∗q /T in Fig. 10 as a function of ρ/T 3 for each temperature T /Tc (β). Dashed lines are cubic spline interpolations of these results. The dot-dashed line is the value of the free quark-gluon gas in the continuum theory, µq 1 µq 3 ρ . (32) = N + f T3 T π2 T

In this calculation, we neglected the F dependence of hθ2 i(P,F ) because the values of P and F are strongly correlated. The systematic error due to this approximation is discussed in Appendix C. The error seems to be small. From Fig. 10, we find that a qualitative feature of µ∗q /T changes around T /Tc ∼ 0.8, i.e. µ∗q /T increases monotonically as ρ increases above 0.8, whereas it shows an s-shape below 0.8. This means that there is more than one value of ρ/T 3 for one value of µ∗q /T below T /Tc ∼ 0.8. This is a signature of a first order phase transition. The critical point in the (T, µq ) plane is estimated in [12] by calculating the effective potential in terms of the plaquette value using the same configurations. The estimation is (T /Tc , µq /T ) ≈ (0.76, 2.5). Because the estimation from the effective potential is rather ambiguous, the difference between the new and old results of the critical temperature may be a systematic error. The critical value of µ∗q /T is about 2.4. This is almost consistent with the previous result by the different method. The error from the truncation of the Taylor expansion of the quark determinant is discussed in [12]. The difference between the results at O(µ4q ) and O(µ6q ) is found to be small at µq /T = 2.5 for the data used in this study. Therefore, the error due to the truncation would not affect the qualitative conclusions. Although further studies including justifications of the approximations used in this analysis are necessary for more quantitative investigations, this result suggests the existence of the first order phase transition line in the (T, µq ) plane. VI.

CONCLUSIONS

We studied the canonical partition function as a function of ρ/T 3 performing an inverse Laplace transformation. We analyzed the data obtained with two-flavors of p4-improved staggered quarks in [6] and calculated the derivative of the canonical partition function with respect to ρ. The problems in this calculation were discussed. To avoid the problems, we adopted the following approximations. First, we estimate the quark determinant from the data of a Taylor expansion up to O(µ6q ) because the direct calculation of the quark determinant is still difficult except on a small lattice. Although terms of higher than µ6q are omitted, this analysis is valid in the low density region. Second, we use a saddle point approximation for the inverse transformation, assuming the volume is sufficiently large. Third, we assume that the probability distribution of the complex phase of the operator in the calculation of µ∗q /T 3 can be well approximated by a Gaussian function. Using multiparameter reweighting method, we combined the configurations generated by µq = 0 simulations at 16 simulation points (β) which cover a wide range of the temperature. It is found that the increase of µ∗q /T 3 becomes larger as the temperature decreases in the low density region. However, the contribution from the configurations generated at low temperature gradually decreases in the measurement of µ∗q /T 3 as ρ/T 3 increases even at low temperature. And, µ∗q /T 3 approaches the value of the free quark gas in the high density limit for all temperatures investigated in this study. The most interesting result is that µ∗q /T 3 as a function of ρ/T 3 shows an s-shape at T < ∼ 0.8. This means that the effective potential Veff in terms of the density has two minima. Therefore, this result strongly suggests the existence of the first order phase transition line in the low temperature and high density region. Since the data we used in this study is obtained by a simulation with much heavier quark masses than the physical quark masses, simulations near the physical mass point are very important. It is also necessary to increase the accuracy of the approximations we have used in this study. Acknowledgments

I would like to thank F. Karsch, K. Kanaya, T. Hatsuda, S. Aoki, T. Izubuchi, and Y. Hidaka for discussions and comments. This work has been authored under Contract No. DE-AC02-98CH10886 with the U.S. Department of

13 Energy. I also wish to thank the Sumitomo Foundation for their financial assistance (No. 050408). Appendix A: Multi-parameter (β) reweighting method

We discuss a method to combine configurations obtained by simulations with many β, following the method by Ferrenberg and Swendsen [22]. We define a β-independent distribution function of the plaquette value, Z Ω(P ′ ) = DU δ(P − P ′ )(det M (0))Nf . (33) The relation between Ω(P ) and w(P, β) in Eq. (16) is Ω(P ) exp(6βNsite P ) = w(P, β). Then, the expectation value of an operator O as a function of P is given by the equation Z Z X 1 1 O, ZGC = Ω(P )e6βNsite P dP, (34) O(P )Ω(P )e6βNsite P dP ≈ hOiβ = ZGC Nconf β,{conf.}

P where Nconf. is the number of configurations, and β,{conf.} denotes the sum of O over all configurations generated in a simulation at β. The expectation value hOiβ is also calculated from the data obtained by more than one simulation point, βi (i = 1, 2, · · · , Nβ ). The Eq. (34) is evaluated by Z

O(P )Ω(P )e6βNsite P dP = ≈

Nβ X i=1

Ni ZGC (βi )

Nβ X

X

i=1 βi ,{conf.}

Z

O(P )Ω(P ) PNβ

j=1

O PNβ

j=1

e6Nsite (βi +β)P −1 Nj e6Nsite βj P ZGC (βj )

e6Nsite βP

dP

(35)

−1 Nj e6Nsite βj P ZGC (βj )

where Ni is the number of configurations at simulation points βi . Hence, hOiβ ≈

hOG(β, P )iall , hG(β, P )iall

(36)

Here, the weight factor G(β, P ) is G(β, P ) = PNβ

i=1

e6Nsite βP −1 Ni e6Nsite βi P ZGC (βi )

,

(37)

and h· · · iall means the average over all configurations generated at all βi . The partition function ZGC (βi ) is determined by a consistency condition for each i, ZGC (βi ) =

Nβ X

X

j=1 βj ,{conf.}

G(βi , P ) = hG(βi , P )iall .

(38)

This equation can be solved except for the normalization factor. We should note that G(β, P ) is independent of the simulation points βi at which the operators are measured, and important configurations for each calculation are selected by the weight factor automatically. Appendix B: Effect from non-Gaussian terms in the phase factor

We estimate the phase factor when the distribution is slightly different from Gaussian. We consider a distribution function with small a4 (P, F )/a2 (P, F ), w(P, ¯ F, θ) ≈

r

a2 π

3a4 1− 2 +O 4a2

"

a4 a2

2 #!−1

w ¯′ (P, F )e−(a2 θ

2

+a4 θ 4 )

.

(39)

14 In this case, the phase factor, exp[−1/(4a2 )], changes to Z

" #! r −1 2 2 4 3a 3a a 1 a4 a 4 4 4 2 1 − 2 + ··· . (40) + 3− +O eiθ e−a2 θ −a4 θ dθ ≈ w ¯′ (P, F ) exp − w ¯′ (P, F ) 4 π 4a2 4a2 4a2 16a2 a2

and also the expectation value of θ2 for fixed P and F becomes 2

hθ i(P,F )

1 3a4 = − 3 +O 2a2 2a2

"

a4 a2

2 #

.

(41)

From this equation, the term of 3a4 /(4a32 ) in Eq. (40) is absorbed into hθ2 i/2, hence the leading contribution from a4 in the phase factor is exp[−a4 /(16a42 )]. The value of a4 can be evaluated by the Binder cumulant,

4 " # 2 θ (P,F ) a4 6a4 (P, F ) θ . (42) B4 ≡ + O = 3 − a22 (P, F ) a2 hθ2 i2(P,F ) 2 2 Because a−1 2 ∼ hθ i ∼ O(µq ) for the chemical potential at the saddle point z0 , the effect from a4 becomes larger as 4 6 the density increases. Therefore, for the case of a4 /a2 < ∼ O(1) in the µq = 0 limit, exp[−a4 /(16a2 )] is O(µq ) at most, i.e.

F iθ e e (T,µ =0) = eF exp −hθ2 i(P,F ) /2 + O(µ6q ) (T,µ =0) . (43) q

q

This argument suggests that the approximation by the Gaussian distribution is valid for the investigation of the low density region even if a4 is nonzero, however the estimation of the range of ρ in which the non-Gaussian contribution is small may be important as well as the application range of the Taylor expansion in µq . Appendix C: Suppression factor from complex phase fluctuation

In the calculation of µ∗q /T , we need to calculate the suppression factor from the complex phase fluctuation, exp[−hθ2 i(P, F, ρ/T 3 )]. This factor should be a function of P, F and ρ/T 3, however we neglected the F dependence in the calculation of Fig. 10. In this section, we discuss the error from this approximation. As shown in Fig. 5, the values of P and F are strongly correlated. The solid line is the mean value of F among the configurations having the plaquette value P for each ρ/T 3 . The dashed lines above and below the solid line show the fluctuations. Most of the configurations characterized by P and F are distributed in the narrow region between the two dashed lines. Therefore, we have approximated hθ2 i as a function of only P and ρ/T 3 in Fig. 8. To estimate the importance of the F dependence, we calculate µ∗q /T using another approximation of hθ2 i which changes according to F of each configuration. From the data of hF iP as a function of ρ/T 3 for each P , which is shown in Fig. 5, we find ρ/T 3 which gives F for each configuration and find the value of hθ2 i at this ρ/T 3 using the data of Fig. 8. We then obtain hθ2 i as a function of P and F , but the ρ/T 3 dependence is neglected. We calculate µ∗q /T using this hθ2 i(P,F ) and compare the previous result to estimate the systematic error due to the approximation in hθ2 i. The result is shown in Fig. 11. The dotted lines are the spline interpolation in Fig. 10. The difference between the results of µ∗q /T in Figs. 10 and 11 seems to be small. This suggests that the determination of hθ2 i with three parameter (P, F, ρ/T 3 ) is not very important for the qualitative argument. Note: In the jackknife error estimation of this calculation, we have neglected the dispersion of hθ2 i(P,F ) among the jackknife ensemble. Therefore, statistical errors in Fig. 11 are smaller than those in Fig. 10. The small error does not mean that the analysis in this appendix gives better results with smaller statistical error. In fact, the errors in Fig. 10 become the same size if we neglect the dispersion of hθ2 i in the jackknife analysis.

[1] M. Asakawa and K. Yazaki, Nucl. Phys. A 504, 668 (1989). [2] A. Barducci, R. Casalbuoni, S. De Curtis, R. Gatto and G. Pettini, Phys. Lett. B 231, 463 (1989); Phys. Rev. D 41, 1610 (1990); A. Barducci, R. Casalbuoni, G. Pettini and R. Gatto, Phys. Rev. D 49, 426 (1994). [3] M. Stephanov, K. Rajagopal and E. Shuryak, Phys. Rev. Lett. 81, 4816 (1998). [4] Z. Fodor and S. Katz, JHEP 0203, 014 (2002); JHEP 0404, 050 (2004).

15

µ*q/T

4

3 T/Tc=0.76 T/Tc=0.80 T/Tc=0.83 T/Tc=0.87 T/Tc=0.90 T/Tc=0.94 T/Tc=0.98 T/Tc=1.02 T/Tc=1.07 T/Tc=1.11 T/Tc=1.16

2

1 free gas limit 0

0

2

4

8

6

10

12

3

ρ/T

FIG. 11: Derivative of ln ZC as a function of the quark number density. The estimation of hθ2 i(P, F ) is different from that of Fig. 10.

[5] C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and C. Schmidt, Phys. Rev. D 68, 014507 (2003). [6] C.R. Allton, M. D¨ oring, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and K. Redlich, Phys. Rev. D 71, 054508 (2005). [7] Ch. Schmidt, C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch and E. Laermann, Nucl. Phys. B (Proc. Suppl.) 119, 517 (2003); F. Karsch, C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, E. Laermann and Ch. Schmidt, Nucl. Phys. B (Proc. Suppl.) 129, 614 (2004); S. Ejiri, C.R. Allton, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and Ch. Schmidt, Prog. Theor. Phys. Suppl. 153, 118 (2004). [8] P. de Forcrand and O. Philipsen, Nucl. Phys. B 673, 170 (2003); JHEP 0701, 077 (2007) [9] J.B. Kogut and D.K. Sinclair, arXiv:0712.2625. [10] R.V. Gavai and S. Gupta, Phys. Rev. D 71 114014 (2005). [11] S. Ejiri, Phys. Rev. D 73, 054502 (2006), [12] S. Ejiri, Phys. Rev. D 77, 014508 (2008). [13] D.E. Miller and K. Redlich, Phys. Rev. D 35, 2524 (1987). [14] P.E. Gibbs, Phys. Lett. B 172, 53 (1986). [15] I.M. Barbour, S.E. Morrison, E.G. Klepfish, J.B. Kogut, M.-P. Lombardo, Phys. Rev. D 56, 7063 (1997); Nucl. Phys. B(Proc. Suppl.) 60A, 220 (1998). [16] A. Hasenfratz and D. Toussaint, Nucl. Phys. B 371, 539 (1992). [17] A. Alexandru, M. Faber, I. Horvath and K.-F. Liu, Phys. Rev. D 72, 114513 (2005). [18] S. Kratochvila and P. de Forcrand, PoS (LAT2005) 167 (2005); Nucl. Phys. B (Proc. Suppl.) 153, 62 (2006). [19] C.R. Allton, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann, Ch. Schmidt and L. Scorzato, Phys. Rev. D 66, 074507 (2002). [20] F. Csikor, G. I. Egri, Z. Fodor, S. D. Katz, K. K. Szabo and A. I. Toth, JHEP 0405 (2004) 046. [21] S. Ejiri, Phys. Rev. D 69, 094506 (2004). [22] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 63, 1195 (1989). [23] A. Roberge and N. Weiss, Nucl. Phys. B 275, 734 (1986). [24] A.M. Ferrenberg and R.H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988). [25] A. Gocksch, Phys. Rev. Lett. 61, 2054 (1988). [26] R. Aloisio, V. Azcoiti, G. Di Carlo, A. Galante, A.F. Grillo, Phys. Rev. D 61, 111501(R) (2000). [27] K.N. Anagnostopoulos and J. Nishimura, Phys. Rev. D 66, 106008 (2002); J. Ambjorn, K.N. Anagnostopoulos, J. Nishimura and J.J.M. Verbaarschot, JHEP 0210, 062 (2002). [28] T. Takaishi, Mod. Phys. Lett. A 19,909 (2004) [29] Z. Fodor, S. Katz and C. Schmidt, JHEP 0703, 121 (2007). [30] K. Splittorff, PoS LAT2006, 023 (2006); K. Splittorff and J.J.M. Verbaarschot, Phys. Rev. Lett. 98, 031601 (2007). [31] K. Splittorff and J.J.M. Verbaarschot, Phys. Rev. D 77, 014514 (2008).