Auxiliary field Monte-Carlo simulation of strong coupling lattice QCD for QCD phase diagram

We study the QCD phase diagram in the strong coupling limit with fluctuation effects by using the auxiliary field Monte-Carlo method. We apply the chiral angle fixing technique in order to obtain finite chiral condensate in the chiral limit in finite volume. The behavior of order parameters suggests that chiral phase transition is the second order or crossover at low chemical potential and the first order at high chemical potential. Compared with the mean field results, the hadronic phase is suppressed at low chemical potential, and is extended at high chemical potential as already suggested in the monomer-dimer-polymer simulations. We find that the sign problem originating from the bosonization procedure is weakened by the phase cancellation mechanism; a complex phase from one site tends to be canceled by the nearest neighbor site phase as long as low momentum auxiliary field contributions dominate.


I. INTRODUCTION
Quantum Chromodynamics (QCD) phase diagram is attracting much attention in recent years. At high temperature (T ), there is a transition from quark-gluon plasma (QGP) to hadronic matter via the crossover transition, which was realized in the early universe and is now extensively studied in high-energy heavy-ion collision experiments at RHIC and LHC. At high quark chemical potential (µ), we also expect the transition from baryonic to quark matter, which may be realized in cold dense matter such as the neutron star core. Provided that the high density transition is the first order, the QCD critical point (CP) should exist as the end point of the first order phase boundary. Large fluctuations of the order parameters around CP may be observed in the beam energy scan program at RHIC.
The Monte-Carlo simulation of the lattice QCD (MC-LQCD) is one of the first principle non-perturbative methods to investigate the phase transition. We can obtain various properties of QCD: hadron masses and interactions, color confinement, chiral and deconfinement transitions, equation of state, and so on. We can apply MC-LQCD to the low µ region, but not to the high µ region because of the notorious sign problem. The fermion determinant becomes complex at finite µ, then the statistical weight is reduced by the average phase factor e iθ , where θ is the complex phase of the fermion determinant. There are many attempts to avoid the sign problem such as the reweighting method [1], the Taylor expansion method [2], the analytic continuation from imaginary chemical potential [3], the canonical ensemble method [4], the fugacity expansion [5], the histogram method [6], and the complex Langevin method [7]. Many * Electronic address: t-ichi@ruby.scphys.kyoto-u.ac.jp † Electronic address: ohnishi@yukawa.kyoto-u.ac.jp ‡ Electronic address: takashi-nakano@kke.co.jp of these methods are useful for µ/T < 1, while it is difficult to perform the Monte-Carlo simulation in the larger µ region.
Recent studies suggest that CP may not be reachable in phase quenched simulations [8]: In the phase quenched simulation for N f = 2, the sampling weight at finite µ is given as | det D(µ)| 2 = det D(µ)(det D(µ)) * = det D(µ) det D(−µ * ), where D represents the fermion matrix for a single flavor. The phase quenched fermion determinant for real quark chemical potential µ d = µ u = µ ∈ R is the same as that at finite isospin and vanishing quark chemical potentials, µ d = −µ u = µ. Thus the phase quenched phase diagram in the temperaturequark chemical potential (T, µ) plane would be the same as that in the temperature-isospin chemical potential (T, δµ) plane, as long as we can ignore the mixing of u and d condensates. We do not see any critical behavior in the finite δµ simulations outside of the pion condensed phase [9]. By comparison, pion condensed phase appears at large δµ, where the above correspondence does not apply. We may have CP inside the pion condensed phase. Gauge configurations in the pion condensed phase, however, would be very different from those of compressed baryonic matter which we aim to investigate. Therefore, we need to find methods other than the phase quenched simulation in order to directly sample appropriate configurations in cold dense matter for the discussion of CP and the first order transition.
The strong coupling lattice QCD (SC-LQCD) is one of the methods to study finite µ region based on the strong coupling expansion (1/g 2 expansion) of the lattice QCD. There are some merits to investigate QCD phase diagram using SC-LQCD, while the strong coupling limit (SCL) is the opposite limit of the continuum limit. First, the effective action is given in terms of color singlet components, then we expect suppressed complex phases of the fermion determinant and a milder sign problem. We obtain the effective action by integrating out the spatial link variables before the fermion field integral. This point is different from the standard treatment of MC-LQCD, in which we integrate out the fermion field before the link integral. Second, we can obtain insight into QCD phase diagram from the mean-field studies at strong coupling. The chiral transition has systematically and analytically been studied in the strong coupling expansion (1/g 2 expansion) under the mean-field approximation: the strong coupling limit (leading order, O(1/g 0 )) [10][11][12][13][14][15][16][17], the nextto-leading order (NLO, O(1/g 2 )) [12][13][14][15][16][17], and the nextto-next-to-leading order (NNLO, O(1/g 4 )) [15,17].
It is necessary to go beyond the mean-field treatment and to include the fluctuation effects of the order parameters for quantitative studies of the finite density QCD. Monomer-dimer-polymer (MDP) simulation is one of the methods beyond the mean-field approximation. We obtain the effective action of quarks after the link integral, and evaluate the fermion integral by summing up monomer-dimer-polymer configurations [18]. The phase diagram shape is modified to some extent, compared with the mean-field results on an isotropic lattice: the chiral transition temperature is reduced by 10-20 % at µ = 0, and the hadronic phase expands to higher µ direction by 20-30 % [19]. Until now, we can perform MDP simulations only in the strong coupling limit, 1/g 2 = 0, and the finite coupling corrections are evaluated in the reweighting method [20]. Since both finite coupling and fluctuation effects are important to discuss the QCD phase diagram, we need to develop a theoretical framework which includes both of these effects.
In this work, we study the QCD phase diagram by using an auxiliary field Monte-Carlo (AFMC) method as a tool to take account of the fluctuation effects of the auxiliary fields. AFMC is widely used in nuclear manybody problems [21,22] and in condensed matter physics such as ultra cold atom systems [22,23]. In AFMC, we introduce the auxiliary fields to decompose the fermion interaction terms and carry out the Monte-Carlo integral of auxiliary fields, which is assumed to be static and constant in the mean-field approximation. We can thus include the fluctuation effects of the auxiliary fields in AFMC beyond the mean-field approximation.
Another important aspect of this paper is how to fix the chiral angle, the angle between the scalar and pseudoscalar modes. In finite volume, symmetry of the theory is not broken spontaneously and an order parameter, in principle, vanishes. In spin systems, a root mean square order parameter is applied to obtain the appropriate order parameter [24]. We here use a similar method, chiral angle fixing (CAF). In CAF, we rotate all fields by the chiral angle, and obtain quantities by using rotated new fields.
This paper is organized as follows. In Sec. II, we explain the formulation of AFMC in SC-LQCD. In Sec. III, we show the numerical results on the order parameters, phase diagram, and the average phase factor. In Sec. IV, we numerically confirm a source of the sign problem in AFMC, and discuss the order of the phase transition based on the volume dependence of the chiral suscep-tibility. In Sec. V, we devote ourselves to a summary and discussion.

A. Lattice action
We here consider the lattice QCD with one species of unrooted staggered fermion for color SU (N c ) in the anisotropic Euclidean spacetime. Throughout this paper, we work in the lattice unit a = 1, where a is the spatial lattice spacing, and the case of color SU(N c = 3) in 3+1 dimension (d = 3) spacetime. Temporal and spatial lattice sizes are denoted as N τ and L, respectively.
The partition function and action are given as, where χ x , U ν,x , U Pτ and U Ps represent the quark field, the link variable, and the temporal and spatial plaquettes, respectively. η j,x = (−1) x0+···+xj−1 is the staggered sign factor, and V ± x and M x are mesonic composites. Chemical potential µ is introduced in the form of the temporal component of vector potential. The physical lattice spacing ratio is introduced as f (γ) = a phys s /a phys τ . The lattice anisotropy parameters, γ and ξ, are introduced as modification factors of the temporal hopping term of quarks and the temporal and spatial plaquette action terms. Temporal and spatial plaquette couplings should satisfy the hypercube symmetry condition in the isotropic limit (ξ → 1), g τ (g 0 , 1) = g s (g 0 , 1) = g 0 . In the continuum limit (a → 0 and g 0 → 0), two anisotropy parameters should correspond to the physical lattice spacing ratio, f (γ) = γ = ξ, when we construct lattice QCD action requiring a phys s /a phys τ = γ in the continuum region, then we can define temperature as T = f (γ)/N τ = γ/N τ . By comparison, it seems to be more appropriate to define temperature as T = γ 2 /N τ due to quantum corrections in the strong coupling limit (SCL) [14]. For example, the critical temperature is predicted to be proportional to γ 2 rather than γ in the mean field treatment in SCL [14]. We follow this argument and adopt f (γ) = γ 2 .
In SCL, we can ignore the plaquette action terms S G , which are proportional to 1/g 2 . The above lattice QCD action in the chiral limit m 0 → 0 has chiral symmetry

B. Effective action
In the present formulation, we have four main steps to obtain physical observables. First, we integrate out the lattice partition function over spatial link variables in the strong-coupling limit. Second, we introduce the auxiliary fields for the mesonic composites and convert the four-Fermi interaction terms to the fermion bilinear form. Third, we perform the integral over the fermion fields and temporal link variables analytically, and obtain the effective action of the auxiliary fields. Finally, we carry out the Monte-Carlo integral over the auxiliary fields.
In the second step, we transform the four-Fermi interactions, the second terms in Eq. (9), to the fermionbilinear form. By using the Fourier transformation in spatial directions M x=(x,τ ) = k e ik·x M k,τ , the interaction terms read where f (k) = j cos k j andk = k + (π, π, π). For later use, we divide the momentum region into the positive (f (k) > 0) and negative (f (k) < 0) modes. In last line of Eq. (10), we use the relation f (k) = −f (k). We introduce the auxiliary fields via the extended Hubbard-Stratonovich (EHS) transformation [16]. We can bosonize any kind of composite product by introduc-ing two auxiliary fields simultaneously, where ψ = ϕ + iφ and dψ dψ * = dReψ dImψ = dϕdφ. When the two composites are the same, A = B, Eq. (11) corresponds to the bosonization of attractive interaction terms. For the bosonization of interaction terms which lead to repulsive potential in the mean-field approximation, we need to introduce complex number coefficients, The bosonization of the interaction terms in Eq. (10) is carried out as where k,τ , and α = L 3 /4N c . We introduce σ k,τ and π k,τ as the auxiliary fields of M k,τ and iM −k,τ , respectively. σ k,τ (π k,τ ) includes the scalar (pseudoscalar) and some parts of higher spin modes. By construction, σ k,τ and π k,τ satisfy the relation σ −k,τ = σ * k,τ and π −k,τ = π * k,τ , which means that σ x , π x ∈ R.
In the third step, we carry out the Grassmann and temporal link (U 0 ) integrals analytically [12][13][14]. We find the partition function and the effective action as, where is a known function of m x and can be obtained by using a recursion formula [12][13][14], as summarized in Appendix B. When m x=(x,τ ) is independent of τ (static), we obtain X Nτ = 2 cosh(N τ arcsinh (m x /γ)).
In the last step, we carry out AFMC integral [25,26]. We numerically integrate out the auxiliary fields (σ k,τ , π k,τ ) based on the auxiliary field effective action Eq. (19) by using the Monte-Carlo method, then we could take auxiliary field fluctuation effects into account.
When we perform integration, we have a sign problem in AFMC [25,26]. The effective action S AF eff in Eq. (19) contains the complex terms X Nτ via the spatial diagonal parts of the fermion matrix I x = m x /γ. Auxiliary fields are real in the spacetime representation, σ x , π x ∈ R, but the negative auxiliary field modes appear with imaginary coefficients as iε x π x , which come from the EHS transformation. The imaginary part of the effective action gives rise to a complex phase in the statistical weight exp(−S AF eff ), and leads to the statistical weight cancellation.
It should be noted that the weight cancellation is weakened in part by the phase cancellation mechanism in low momentum auxiliary field modes. In AFMC, the fermion determinant is decomposed into the one at each spatial site. Since negative modes π k,τ involve iε x , the phase on one site from low momentum π k,τ modes tend to be canceled by the phase on the nearest neighbor site. Thus we could expect that the statistical weight cancellation is not severe when low momentum modes mainly contribute. By comparison, strong weight cancellation might arise from high momentum modes. We discuss the contributions from high momentum modes in Sec. IV B.
While we have the sign problem in AFMC, we anticipate that we could study the QCD phase diagram since the long wave modes are more relevant to phase transition phenomena. We show the results of the QCD phase transition phenomena based on AFMC in the next section, Sec. III.

III. QCD PHASE DIAGRAM IN AFMC
We show numerical results in the chiral limit (m 0 = 0) on 4 3 × 4, 6 3 × 4, 6 3 × 6 and 8 3 × 8 lattices. We have generated the auxiliary field configurations at several temperatures on fixed fugacity (fixed µ/T ) lines. We here assume that temperature is given as T = γ 2 /N τ [14]. Statistical errors are evaluated in the jack-knife method; we consider an error to be the saturated value after the autocorrelation disappears as shown later in Fig. 2.

A. Chiral Angle Fixing
It is a non-trivial problem how to describe the spontaneous symmetry breaking in Monte-Carlo calculations on a finite size lattice: the expectation value of the order parameter generally vanishes since the distribution is symmetric under the transformation. Rigorously, we need to take the thermodynamic limit with explicit symmetry breaking term, and to take the limit of the vanishing explicit breaking term, as schematically shown in Fig. 1 in the case of chiral symmetry. This procedure is time consuming and is not easy to carry out when we have the sign problem.
We here propose a chiral angle fixing (CAF) method as a prescription to calculate the chiral condensate on a finite size lattice. The effective action Eq. (9) is invariant under the chiral transformation, The chiral symmetry is kept in the bosonized effective action by introducing the chiral U(1) transformation for auxiliary fields as, where (σ k , π k ) are the temporal Fourier transform of In order to obtain the chiral condensate rigorously, we need to put a finite mass, first take thermodynamic limit and finally take the chiral (massless) limit as shown in the upper panels. In CAF, we take chiral rotation to make the π0 field vanish, and we get the finite chiral condensate (center bottom panel), which would be close to the correct value.
(σ k,τ , π k,τ ), Because of the chiral symmetry, the chiral condensate σ 0 vanishes as long as the auxiliary field configurations are taken to be chiral symmetric, as explicitly shown in Appendix A.
In order to avoid the vanishing chiral condensate, we here utilize CAF. We rotate σ 0 and π 0 modes toward the positive σ 0 direction as schematically shown in Fig. 1. All the other fields are rotated with the same angle, −α = − arctan(π 0 /σ 0 ), in each Monte-Carlo configuration. We use these new fields to obtain order parameters, susceptibilities, and other quantities, and eventually obtain finite chiral condensate. Chiral condensate obtained in CAF should mimic the spontaneously broken chiral condensate in the thermodynamic limit. Similar prescriptions are adopted in other field of physics. For example, we take a root mean square order parameter to obtain the appropriate value in spin systems [24].

B. Sampling and Errors
We generate auxiliary field configurations by using the Metropolis sampling method. We generate Markov chains starting from two types of initial conditions: the Wigner phase (σ x = 0.01, π x = 0) and the Nambu-Goldstone (NG) phase (σ x = 2, π x = 0) initial conditions.
For each τ , we generate a candidate auxiliary field configuration (σ ′ k,τ , π ′ k,τ ) by adding random numbers to the current configuration (σ k,τ , π k,τ ) for all spatial momenta k at a time, and judge whether the new configuration is accepted or not. Since it is time consuming to update each auxiliary field mode separately, we update all spatial momentum modes in one step at the cost of an acceptance probability. It should be noted that the acceptance probability is larger in the the present (σ k,τ , π k,τ ) sampling procedure in each τ compared with updating (σ k , π k ) in the whole momentum space at a time.
We evaluate errors of calculated quantities in the jackknife method. The evaluated errors of the chiral condensate φ are shown as a function of bin size in the right middle panel of Fig. 2. Since the Metropolis samples are generated sequentially in the Markov chain, subsequent events are correlated. This autocorrelation disappears when the Metropolis time difference is large enough. In the jack-knife method, we group the data into bins and regard the set of configurations except for those in a specified bin as a jack-knife sample. We find that the autocorrelation disappears for the bin size larger than 30 in this case. The jack-knife error increases with increasing bin size, and eventually saturates. We adopt the saturated value of the jack-knife error after the autocorrelation disappears as the error of the calculated quantity as in the standard jack-knife treatment. The errors are found to be small enough, for example ∆φ 0.01, compared with its mean value shown in Fig. 3 and to discuss the phase transition.

C. Order Parameters
In Fig. 3, we show the chiral condensate, φ = σ 0 , and the quark number density ρ q after CAF, as a function of temperature (T ) on a 8 3 × 8 lattice. Necessary formulae to obtain these quantities are summarized in Appendix B. We also show the distribution of φ in Fig. 4. The order parameters, φ and ρ q , clearly show the phase transition behavior. With increasing T for fixed µ/T , the chiral condensate φ slowly decreases at low T , shows rapid or discontinuous decrease at around the transition temperature, and stays to be small at higher T . The quark number density ρ q also shows the existence of phase transition at finite µ. The order of the phase transition can be deduced from the behavior of φ, ρ q and the φ distribution on a small lattice [25,26]. The chiral condensate φ and the quark number density ρ q smoothly change around the (pseudo-)critical temperature (T c ) at small µ/T . Additionally, the φ distribution has a single peak as shown in the top panel of Fig. 4. These observations suggest that the phase transition is crossover or the second order at small µ/T on a large size lattice. We refer to this µ/T region as the would-be second order region.
By comparison, the order parameters show hysteresis behavior in the large µ/T region. As shown by dashed lines in Fig. 3, two distinct results of φ and ρ q depend on the initial conditions, the Wigner phase and the NG phase initial conditions. The temperature of sudden φ change for the NG initial condition is larger than that for the Wigner initial condition. The distribution of φ shows a double peak as shown in the bottom panel of Fig. 4. In terms of the effective potential, the dependence of initial conditions indicates that there exist two local minima, which are separated by a barrier. In the hysteresis region, the transition between the two local minima is suppressed by the barrier and Metropolis samples stay around the local minimum close to the initial condition. At the temperature of sudden φ change, the barrier height becomes small enough for the Metropolis samples to overcome the barrier. These results suggest that the phase transition is the first order at large µ/T . We refer to this µ/T region as the would-be first order region.

D. Phase Diagram
We shall now discuss the QCD phase diagram in AFMC. In Fig. 6, we show the QCD phase diagram for various lattice sizes. We define the (pseudo-)critical temperature T c as a peak position of the chiral susceptibility χ σ shown in Fig. 5 in the would-be second order region. We determine the peak position by fitting the susceptibility with a quadratic function. The errors are comprised of both statistical and systematic errors. We fit χ σ as a function of T with statistical errors obtained in the jack- knife method. In order to evaluate the systematic error, we change the fitting range as long as the fitted quadratic function describes an appropriate peak position. We take notice that we do not fit χ σ as a function of T in each jack-knife sample.
In the would-be first order region of µ/T , we determine the phase boundary by comparing the expectation values of effective action S eff in the configurations sampled from the Wigner and NG phase initial conditions. We define T c as the temperature where S eff with the Wigner initial condition becomes lower than that with the NG initial condition as shown in Fig. 5. We have adopted this prescription, since it is not easy to obtain equilibrium configurations over the two phases when the thermodynamic potential barrier is high. At large µ/T , Metropolis samples in one sequence stay in the local minimum around the initial condition, and we need very large sampling steps to overcome the barrier.
In Fig. 6, we compare the AFMC phase boundary with that in the mean field approximation [11,16,17] and in the MDP simulation [11,19] in the strong coupling limit. Compared with the MF results, T c at low µ is found to be smaller, and NG phase is found to be extended in the finite µ region in both MDP [19] and AFMC. As found in previous works [25,26], the phase boundary is approximately independent of the lattice size in the would-be second order region. The would-be first order phase boundary is insensitive to the spatial lattice size but is found to depend on the temporal lattice size. With increasing temporal lattice size, the transition chemical potential µ c becomes larger, which is consistent with MDP [19]. Phase boundary extrapolated to N τ → ∞ is shown by the shaded area, and is found to be consistent with the continuous time MDP results with the same limit, N τ → ∞ with keeping γ 2 /N τ finite.
Spatial lattice size independence of the phase boundary may be understood as a consequence of almost decoupled pions. The zero momentum pion can be absorbed into the chiral condensate via the chiral rotation and has no effects on the transition. Finite momentum pion modes have finite excitation energy, then we do not have soft modes in the would-be first order region on a small size lattice. For a more serious estimate of the size dependence, we need larger lattice calculations.
We find that the would-be first order phase boundary has a positive slope, dµ/dT > 0, at low T . The Clausius-Clapeyron relation reads dµ/dT | 1st = −(s W − s NG )/(ρ W q − ρ NG q ), where s W,NG and ρ W,NG q are the entropy density and quark number density in the Winger and NG phases, respectively. Since ρ q is higher in the Winger phase as shown in Fig. 3, the entropy density should be smaller in the Winger phase. This is because ρ q is close to the saturated value, ρ q ∼ 3 = N c , in the Wigner phase, then the entropy is carried by the hole from the fully saturated state. Similar behavior is found in the mean-field treatment in the strong coupling limit [11]. In order to avoid the quark number density saturation, which is a lattice artifact, we may need to adopt a larger N τ [19] or to take account of finite coupling effects [16,17].

E. Average Phase Factor
In Fig. 7, we show the average phase factor e iθ as a function of T on 8 3 × 8 and 4 3 × 4 lattices, where θ is a complex phase of the fermion determinant in each Monte-Carlo configuration. The average phase factor shows the severity of the statistical weight cancellation; we have almost no weight cancellation when e iθ ≃ 1, and the weight cancellation is severe in the cases where e iθ ≃ 0. The average phase factor has a tendency to increase at large µ except for the transition region. This trend can be understood from the effective action in Eq. (19). The complex phase appears from X Nτ terms containing auxiliary fields, and their contribution generally becomes smaller compared with the chemical potential term, 2 cosh(3N τ µ/γ 2 ), at large µ. In the phase transition region, fluctuation effects of the auxiliary fields are decisive and finite momentum auxiliary fields might contribute significantly, which leads to a small average phase factor. The average phase factor on a 4 3 ×4 lattice, e iθ 0.9, is practically large enough to keep statistical precision. By comparison, the smallest average phase factor on a 8 3 ×8 lattice is around 0.1 at low temperature on a µ/T = 2.4 line. Even with this average phase factor, uncertainty of the phase boundary shown in Fig. 6 is found to be small enough to discuss the fluctuation effects.
We show the severity of the sign problem in AFMC in Fig. 8. The severity is characterized by the difference of the free energy density in full and phase quenched (p.q.) MC simulations, ∆f = f full − f p.q. which is related to the average phase factor, e −Ω∆f = e iθ p.q. , where Ω = N τ L 3 is the spacetime volume. While ∆f takes smaller values on a 4 3 × 4 lattice, it takes similar values on lattices with larger spatial size L ≥ 6. We expect that ∆f in AFMC for larger lattices would take values similar to those on a 8 3 × 8 lattice.
We find that ∆f in AFMC is about twice as large as that in MDP when we compare the results at similar (µ, T ) [19]. It means that the sign problem in AFMC is more severe than that in MDP. It is desired to develop a scheme to reduce ∆f in AFMC on larger lattices. In Sec. IV B, we search for a possible way to weaken the weight cancellation by cutting off high momentum auxiliary fields.

A. Volume Dependence of Chiral susceptibility
We investigate the volume dependence of the chiral susceptibility to discuss the phase transition order in the low chemical potential region. We expect that the phase transition is the second order at small µ/T according to the mean-field results and O(2) symmetry arguments. The latter states that the fluctuation induced first order phase transition is not realized as for O(2) symmetry [27].
In Fig. 9, we show the chiral susceptibility for fixed µ/T = 0.2 on various size lattices. From this comparison, we find that χ σ has a peak at the same T for different lattice sizes, and that the peak height on 6 3 × 4 and 6 3 × 6 lattices are almost the same. These two findings suggest that it is reasonable to define the temperature as T = γ 2 /N τ in the strong coupling limit. We also find that the peak height of the susceptibility increases with increasing spatial lattice size. The divergence of the susceptibility in the thermodynamic limit signals the first or second order phase transition.
In order to find finite size scaling of chiral susceptibility, we plot 1/χ σ as a function of inverse spatial lattice volume in Fig. 10. The chiral susceptibility is proportional to spatial volume V = L 3 in the first order phase transition region and to V (2−η)/3 in the second order phase transition region for a d = 3 O(2) spin systems, where the O(2) critical exponent is η = 0.0380(4) [28]. By comparison, χ σ does not diverge when the transition is crossover. It seems to suggest that the chiral phase transition at low µ is not the first order, and we cannot exclude the possibility of the crossover transition with the present precision in comparison with the above three scaling functions shown in Fig. 10 in AFMC. The current analysis implies that the phase transition is the second order or crossover phase transition. In order to conclude the order of the phase transition firmly, we need higherprecision and larger volume calculations.

B. High momentum mode contributions
We quantitatively examine the influence of high momentum auxiliary field modes on the average phase factor and the order parameters. For this purpose, we compare the results by cutting off high momentum auxiliary field modes having j sin 2 k j > Λ, where Λ is a cutoff parameter. The parameter Λ is varied in the range 0 ≤ Λ ≤ d = 3 to examine their cutoff effects; we include all Monte-Carlo configurations when Λ = 3, while we only take account of the lowest momentum modes when Λ = 0.
The average phase factor might become large in the cases where high momentum mode contributions are negligible as discussed in Sec. II B, so we anticipate that the weight cancellation becomes weaker for smaller Λ. In the left top panel of Fig. 11, we show the Λ dependence of the average phase factor on a 8 3 ×8 lattice for µ/T = 0.6. The average phase factor has a large value when Λ → 0, where we improve the statistical weight cancellation. These results are consistent with our expectation for the statistical weight cancellation with high momentum modes. We could here conclude that high momentum modes are closely related to severe weight cancellation.
In the right bottom panel of Fig. 11, we show the chiral condensate φ on a 8 3 × 8 lattice for µ/T = 0.6. We here utilize φ = τ σ k=0,τ /N τ . This expression is equivalent to Eq. (24) for full configurations. The chiral condensate does not depend on the parameter Λ since the lowest modes of the integration variables (σ k,τ , π k,τ ) in AFMC consist of the scalar and pseudoscalar modes. In Fig. 11, we also plot the cutoff dependence of other quantities: quark number density (ρ q ), chiral susceptibility (χ σ ) and quark number susceptibility (χ µ,µ ). We find that these quantities do not strongly depend on the cutoff as long as Λ ≥ 2. By contrast, the quantities are affected by the cutoff parameter for Λ < 2. We have already known that the average phase factor becomes large if we set Λ ≤ 2.5. Thus, this analysis implies a probable presence of an optimal cutoff Λ o , with which the order parameter values are almost the same as those of the full ensemble results and the reliability of numerical simulation is improved. We conclude that there is a possible way to study the QCD phase diagram for larger lattice by cutting off or approxi-mating the high momentum modes without changing the behavior of the order parameters.

V. SUMMARY
We have investigated the QCD phase diagram and the sign problem in the auxiliary field Monte-Carlo (AFMC) method with chiral angle fixing (CAF) technique. In order to obtain the auxiliary field effective action, we first integrate out spatial link variables and obtain an effective action as a function of quark fields and temporal link variables in the leading order of the 1/g 2 and 1/d expansion with one species of unrooted staggered fermion. By using the extended Hubbard-Stratonovich (EHS) transformation, we convert the four-Fermi interactions into the bilinear form of quarks. The auxiliary field effective action is obtained after analytic integration over the quark and temporal link variables. We have performed the auxiliary field integral using the Monte-Carlo technique.
We have obtained auxiliary field configurations in AFMC and the order parameters: the chiral condensate and quark number density. Both of order parameters show phase transition behavior. In the low chemical potential region, the chiral condensate decreases smoothly with increasing temperature, while the quark number density increases gently. This behavior suggests that the order of the phase transition is the second or crossover, which is consistent with the analysis of the distribution of the chiral condensate. We call the low chemical potential region the would-be second order region. In order to deduce the phase boundary, we here define (pseudo-)critical temperature as a peak position of the chiral susceptibility. One finds that the critical temperature is suppressed compared with the mean-field results on a isotropic lattice and almost independent of lattice size as shown in the monomer-dimer-polymer simulations (MDP) at the would-be second order phase transition [19]. We also give some results of finite size scaling to guess the phase transition order. While one could expect the second order phase transition from the mean-field and O(2) symmetry arguments in the low chemical potential region, it is not yet conclusive to decide whether the transition is the second order or crossover at the present precision.
At high chemical potential, the order parameters show sudden jump and hysteresis, and depend on initial conditions: the Wigner and Nambu-Goldstone initial conditions. The distribution of the chiral condensate has a double peak around the phase transition region. These results imply that the order of the phase transition is the first order owing to the existence of the two local minima with a relatively high barrier compared to the Metropolis jumping width. We call this phase transition the would-be first phase transition in the present paper. We here regard transition temperature as a crossing point of the expectation value of the effective action with two initial conditions. According to our analysis, the Nambu-Goldstone phase is enlarged toward the high chemical potential region compared with the mean-field results. The phase boundary depends very weakly on spatial lattice size and more strongly on temporal lattice size. This behavior is also found in MDP [19].
We find that we have a sign problem in AFMC. The origin of the weight cancellation is the bosonization of the negative modes in the extended Hubbard-Stratonovich (EHS) transformation; an imaginary number must be introduced in the fermion matrix. The fermion determinant becomes complex, and the statistical weight cancellation arises when we numerically integrate auxiliary fields. In our framework, we have a phase cancellation mechanism for low momentum auxiliary fields; a phase on one site is canceled out by the nearest neighbor site phase. We quantitatively show that the high momentum modes contribute to the statistical weight cancellation by cutting off these modes. We also confirm the cutoff dependence on order parameters and susceptibilities. We find that there is a cutoff parameter region where the behavior of the quantities are not altered from the full configurations and the statistical weight cancellation is weakened. Therefore, there is a possibility to investigate phase transition phenomena using cutoff or approximation scheme for high momentum modes.
While we have a sign problem in AFMC, statistical weight cancellation is not serious on small lattices adopted here (∼ 8 3 × 8 size) because of the phase cancellation mechanism for the low momentum modes. The phase boundary in AFMC is found to be consistent with that in MDP [19].
In this paper, we utilize CAF in order to obtain the order parameters and susceptibilities in the chiral limit on a finite size lattice. The chiral condensate in finite volume should vanish in a rigorous sense due to the chiral symmetry between the scalar and pseudoscalar modes. In order to simulate the non-vanishing chiral condensate to be obtained in the rigorous procedure of the thermodynamic limit followed by the chiral limit, the chiral transformation of auxiliary fields are carried out in each configuration so as to fix the chiral angle to be in the real positive direction (positive scalar mode direction). We could evaluate the adequate chiral condensate and chiral susceptibility by using CAF.
The AFMC method could be straightforwardly applied to include finite coupling effects since bosonization technique is applied in the mean-field analysis [16,17]. Both fluctuations and finite coupling effects are important to elucidate features of the phase transition phenomena, so the AFMC would be a possible way to include these two effects at a time. The sign problem might be severer than that in the strong coupling limit when we include finite coupling effects. One of methods to avoid lower numerical reliability is to invoke shifted contour formulation [29]. We hope that we may apply the formulation with finite coupling effects or on a larger lattice. We also obtain appropriate order parameters in a relatively hassle-free CAF method compared to a rigorous way. We might use this CAF method with higher-order corrections in the strong coupling expansion to investigate the phase diagram.
φ k,ω and α k,ω are chiral radius and chiral angle respecting each chiral partner. We find that the chiral condensate ideally vanishes according to Eq. (A1). In CAF, we rotate the negative chiral angle (−α) with respect to all fields and set π 0 = 0. We obtain the finite chiral condensate in the Nambu-Goldstone (NG) phase as The resultant chiral condensate in CAF should simulate the spontaneously broken chiral condensate in the thermodynamic limit.
We have some advantages in CAF. One is that the chiral condensate is finite in the NG phase and the chiral susceptibility may have a peak. In the cases where the chiral condensate vanishes ( σ 0 = 0) because of the chiral symmetry, the chiral susceptibility, which is proportional to ∂ 2 ln Z/∂m 2 0 = −∂ χχ /∂m 0 = ∂ σ 0 /∂m 0 , is expressed as ∂ 2 ln Z/∂m 2 0 = σ 2 0 . We could expect that the chiral susceptibility increases with lower temperature. After we utilize CAF, we obtain the chiral susceptibility with a peak, ∂ 2 ln Z/∂m 2 0 = σ 2 0 − σ 0 2 because of the non-vanishing chiral condensate as shown in Fig. 9. Another merit of CAF is that when we calculate the chiral condensate and the chiral susceptibility, we could take into account the information on the pseudoscalar mode which is mixed with scalar mode in the chiral limit.