Endpoint of the first-order phase transition of QCD in the heavy quark region by reweighting from quenched QCD

1Department of Physics, Niigata University, Niigata 950-2181, Japan 2Graduate School of Science and Technology, Niigata University, Niigata 950-2181, Japan 3Track Maintenance of Shinkansen, Rail Maintenance 1st Department, East Japan Railway Company Niigata Branch, Niigata, Niigata 950-0086, Japan 4Tomonaga Center for the History of the Universe, University of Tsukuba, Tsukuba 305-8571, Japan 5Faculty of Pure and Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan 6Department of Physics, Osaka University, Toyonaka, Osaka 560-0043, Japan 7J-PARC Branch, KEK Theory Center, Institute of Particle and Nuclear Studies, KEK, 203-1, Shirakata, Tokai, Ibaraki 319-1106, Japan 8Graduate School of Education, Hiroshima University, Higashihiroshima, Hiroshima 739-8524, Japan (Dated: December 22, 2019))


I. INTRODUCTION
Quantum chromodynamics (QCD) is known to have a rich phase structure as a function of temperature T , quark chemical potential µ, and quark masses m q . The confined phase at low T and small µ turns into the deconfined phase when T exceeds the transition temperature T c . The nature of this deconfinement phase transition varies as a function of µ and m q . From lattice QCD simulations, the transition is an analytic crossover at small µ around the physical quark masses [1,2] but is expected to change to a first-order transition when we increase µ or vary m q . Determination of the critical point where the crossover changes to the first-order transition is important in understanding the nature of quark matter in heavy-ion collisions [3][4][5][6][7][8][9][10].
In 2+1 flavor QCD, which contains dynamical up, down and strange quarks, there exist two first-order phase transition regions in the quark mass parameter space at µ = 0. One region locates around the light quark limit (chiral limit): When all the three quarks are massless, the spontaneously broken chiral symmetry at low T is recovered by a first-order chiral transition at which confinement is also resolved [11]. Many studies have been made to determine the critical mass where the first-order transition terminates in the quark mass parameter space [6,7,[12][13][14][15][16][17][18][19]. The critical mass was found to be close to the physical point and thus a quantitative determination of it is phenomenologically important. However, so far, the continuum limit of the critical mass has not been obtained conclusively. Another first-order transition region locates around the heavy quark limit [20][21][22][23][24][25][26]: When all the quarks are infinitely heavy, QCD is just the pure gauge SU(3) Yang-Mills theory (quenched QCD), which is known to have a weakly first-order deconfinement transition [27]. This first-order transition changes to crossover when the quark mass becomes smaller than a critical value. Simulations on coarse lattices suggest that the critical quark mass is large. However, its continuum limit is also not well understood.
In this paper, we study the endpoint of the first-order deconfinement transition in the heavy quark region of two and 2+1 flavor QCD, and evaluate the critical quark mass at the endpoint. In our previous papers [23,24], we studied the issue at µ = 0 and finite µ on a 24 3 × 4 lattice. From quenched QCD simulations combined with the hopping parameter expansion, we found that the first-order deconfinement transition in the heavy quark limit becomes weaker and eventually disappears as the quark mass decreases. We determined the location of the critical surface separating the first-order and crossover regions around the heavy quark limit. In this paper, we extend the studies using finer lattices with the temporal lattice size N t = 6 and 8, and compare the results with those at N t = 4 to discuss the lattice spacing dependence of the critical point. We also compute the pseudo-scalar meson mass at the critical point.
This paper is organized as follows. In the next section, we define our model and introduce the histogram method [28,29] to find the endpoint of the first-order phase transition. In Sec. IV, we determine the critical point in two-flavor QCD adopting a leading order (LO) approximation of the hopping parameter expansion. The overlap problem in the reweighting method is also discussed. In Sec. V, influence of the next-to-leading order (NLO) terms of the hopping parameter expansion is evaluated to examine applicability of the hopping parameter expansion in the study of the critical point. The pseudo-scalar meson mass is computed at the endpoint of the first-order phase transition in Sec. VI. The lattice spacing dependence of the critical quark mass is discussed for two flavor QCD. We also determine the boundary of the first-order region in the mass parameter plane of 2 + 1 flavor QCD in Sec. VII. We summarize our results in Sec. VIII. In Appendix A, we introduce the multi-point histogram method used in this study.

II. FORMULATION
The action we study consists of the gauge action and the quark action, where the gauge action is the standard plaquette action given by with β = 6/g 2 the gauge coupling parameter, N site = N 3 s × N t the space-time lattice volume, and P the plaquette operator: Here, U x,µ is the link variable in the µ direction at site x and x +μ the next site in the µ direction from x. For quarks, we adopt the standard Wilson quark action given by where M xy is the Wilson quark kernel with K f the hopping parameter for the f -th flavor.

A. Histogram and effective potential
The order parameter of the deconfinement transition of QCD in the heavy quark limit is given by the Polyakov loop defined aŝ where x is a summation over one time slice. The Polyakov loop provides us with an order parameter for the Z(3) center symmetry, which spontaneously breaks down when T exceeds a critical value. In our previous studies [23,24,29], we found that the histogram of the Polyakov loop is useful in determining the phase structure around the heavy quark limit. 1 We focus on the absolute value of the Polyakov loop |Ω|. Denoting K = (K 1 , · · · , K N f ), we define the histogram of |Ω| by The probability distribution function of |Ω| is given by Z −1 (β, K) w(|Ω|; β, K) with Z the partition function defined by Z(β, K) = w(|Ω|; β, K) d|Ω|, (8) and the naive histogram, obtained by just counting the number of configurations with |Ω| in a simulation, is given by N Z −1 (β, K) w(|Ω|; β, K) with N the total number of configurations. When the transition is of first-order, we expect multiple peaks in the histogram. We can thus detect the endpoint of a first-order transition through a change in the shape of the histogram. For convenience, we define the effective potential of |Ω| as where the constant term is fixed later. The effective potential has two minima when the histogram has two peaks. Using the reweighting method, the parameter β and K can be changed from the simulation point (β 0 , K 0 ) with K 0 = (K 0;1 , · · · , K 0;N f ). The probability distribution of |Ω| at (β, K) is given by where is the expectation value at (β 0 , K 0 ). In principle, Eq. (10) should enable us to predict the shape of w(|Ω|) at any (β, K). However, a statistically reliable data of w(|Ω|) is available only around |Ω| ≈ |Ω| (β 0 ,K 0 ) . When |Ω| (β,K) at the target point (β, K) shifts a lot from |Ω| (β 0 ,K 0 ) , it is not easy to obtain a reliable prediction about the nature of the vacuum at (β, K). This is the overlap problem, which is severe around a first-order transition point on large lattices. The overlap problem can be avoided, to some extent, by combining information obtained at several different simulation points. For the histogram, the multi-point histogram method is developed based on the reweighting method [23,24]. The method is introduced in Appendix A.

B. QCD in the heavy quark region
To investigate the quark mass dependence of the effective potential in the heavy quark region, we evaluate the quark determinant by the hopping parameter expansion in the vicinity of the heavy quark limit K = 0. For each flavor, we have with where (∂M/∂K) xy is the term following K f in the right hand side of Eq. (5). Therefore, the non-vanishing contributions to D n are given by Wilson loops and Polyakov-loop-type loops.
In this study, we compute the leading order (LO) and the next-to-leading order (NLO) terms from the Wilson loops and Polyakov-loop-type loops. The LO contribution consists of the smallest Wilson loop, plaquetteP , defined by Eq. (3), and the shortest Polyakov-loop-type loop, the Polyakov loop Ω, defined by Eq. (6). Because of the anti-periodic boundary condition and gamma matrices in the hopping terms, up to the NLO contributions, Eq. (12) reads where a constant term does not appear because the left hand side vanishes at K = 0. We assume that N t is an even number. The NLO contribution from Wilson loop is given by 6-step Wilson loops: rectangleŴ rec , chair-typeŴ chair , and crown-typeŴ crown . They are illustrated in Fig. 1. The Polyakov-loop-type loops for the NLO contribution are (N t + 2)-step bended Polyakov loopsΩ n with n = 1, · · · , N t /2, which run one step in a space direction, n steps in the time direction and return to the original line, e.g., They are illustrated in Fig. 2. All the Wilson loops and Polyakov-loop-type loops are normalized such thatŴ rec =Ŵ chair =Ŵ crown = 1 andΩ n = 1 in the weak coupling limit, U x,µ = 1. The first term of Eq. (14) is proportional toP and thus can be absorbed by a shift β → β * in , Ω 2 (middle) and Ω 3 (right) for N t = 6. The vertical upward direction is the temporal direction.
the gauge action. Collecting the contribution from all flavors, we find We note that the 6-step Wilson loops are typical operators in improved gauge actions. Thus the contributions of these operators can also be absorbed by a shift of improvement parameters of the gauge action. Because a shift in improvement parameters only affects the amount of lattice discretization errors within the same universality class, the 6-step Wilson loop terms will not affect characteristic physical properties of the system in the continuum limit, such as the critical exponents of the phase transition. In contrast, the terms proportional to the Polyakov-loop-type loops act like external magnetic fields in spin models, and thus may change the nature of the phase transition. Therefore, in our study of NLO contributions, we concentrate on the effects of the Polyakov-loop-type terms on the nature of the phase transition, disregarding the effects of the 6-step Wilson loop terms in Eq. (14). In Sec. IV C, we first study the phase structure only with the LO terms of O(K Nt ) in ln det M . We then study the influence of the NLO terms of O(K Nt+2 ) in Sec. V. On the other hand, the 6-step Wilson loops will do affect detailed properties of the system, such as the values of critical hadron masses. In Sec. VI, we examine their effects on the pseudo-scalar meson mass by comparing with the results of direct full QCD simulations at K = 0.

III. SIMULATION PARAMETERS
We calculate the effective potential of |Ω| in the heavy quark mass region performing simulations of quenched QCD on 24 3 × 6, 32 3 × 6, 36 3 × 6, and 24 3 × 8 lattices. We generate the quenched configurations using the pseudo heat bath algorithm of SU(3) gauge theory with over relaxation. Because the effective potential must be investigated in a wide range of |Ω| to study the nature of the phase transition, we perform simulations at several points of β and combine these data by the multi-point histogram method. Details of the multi-point histogram method are given in Appendix A. We mainly investigate on N t = 6 lattices. The spatial volume dependence is also studied in this case. Adjusting the simulation parameters to the critical point, the lattice spacing a is given by a = (N t T c ) −1 in terms of N t and the transition temperature T c . To study lattice spacing dependence, i.e., N t -dependence, we perform an additional simulation on a lattice with N t = 8, and also use the results of N t = 4 lattices obtained in Ref. [23]. The values of β and the number of configurations are summarized in Table I.
To check validity of the reweighting method with the hopping parameter expansion, we also perform full QCD simulations with two-favors of dynamical Wilson quarks on 16 3 × 32 lattices, In numerical evaluation of the histogram defined by Eq. (7), we adopt a Gaussian approximation for the delta function: . We choose the width parameter ∆ by considering a balance between the resolution of the histogram and its statistical error. The values of ∆ we adopt are given in Table II. The statistical errors are estimated by the jackknife method.

IV. CRITICAL POINT OF TWO-FLAVOR QCD IN THE LEADING ORDER HOPPING PARAMETER EXPANSION
We investigate how the quark mass affects the shape of the effective potential around the firstorder transition point using the reweighting method explained in Sec. II. We first study the case of two-flavor QCD with degenerate u and d quarks. Extension to the case of 2 + 1 flavor QCD is discussed in Sec. VII. We calculate the effective potential V eff at small K by reweighting from K = 0. In the first-order transition region, V eff has two minima. At each K, we first adjust β to the transition point by adjusting β so that two minimum values of V eff are equal. We then measure the difference ∆V eff between the peak height in the middle of V eff and the minimum value as illustrated in Fig. 3. 2 We define the critical point K ct , where the first-order transition line terminates, as K where ∆V eff vanishes. The logarithm of the quark determinant ln det M , Eq. (14), in the reweighting factor is calculated by the LO terms of O(K Nt ), i.e., Polyakov loop, in this section, and the Polykov-looptype loops up to the NLO terms of O(K Nt+2 ) in Sec. V. By comparing the results of K ct using the LO contribution only and that using both LO and NLO contributions, we estimate the truncation error of the hopping parameter expansion.  flat at the minimum at K = 0.135. We plot ∆V eff as a function of K in Fig. 5. Fitting the smallest four data by a linear function (the dashed line), we obtain the critical point K ct = 0.1359 (30).
To study finite-size effect in this result at N t = 6, we perform simulations also on 32 3 × 6 and 36 3 × 6 lattices. The results of V eff (|Ω|) on the 32 3 × 6 lattice is shown in Fig. 6. We find the double-well shape becomes milder as we increase K. The height of the potential barrier ∆V eff is shown in Fig. 7. ∆V eff decreases towards zero at K ∼ 0.12. Fitting the last seven data points by a linear function, we obtain K ct = 0.1286(40), which is roughly the same as K ct obtained on the 24 3 × 6 lattice within the statistical error. The effective potential V eff on the 36 3 × 6 lattice is shown in Fig. 8. The figure shows that the double-well shape of V eff remains up to our largest K. To check reliability of these results, we study the overlap problem in the calculation of V eff . As discussed in Appendix A, when the number of original configurations is not large enough around the peaks of the target histogram, the results by the reweighting method is not reliable. In general, the overlap problem is more likely to occur as the volume increases.
In Figs. 9 and 10, we plot the distribution of (P, |Ω|) from original configurations at K = 0 obtained on the 24 3 × 6 and 32 3 × 6 lattices, respectively. We change color depending on β of the original configurations used in the multi-point histogram method. The bigger circle symbols denote the two peak positions of the reweighted histogram at the transition point at various K. We see that the peak positions remain within the region of the original distribution up to the largest  K we study. We thus conclude that the overlap problem does not contaminate the results on the 24 3 × 6 and 32 3 × 6 lattices. The distribution of (P, |Ω|) on the largest 36 3 × 6 lattice is shown in Fig. 11. The two peak positions of the histogram for K = 0.0, 0.107 and 0.119 are plotted by bigger circle symbols. We note that the peak positions for K = 0.119 are clearly out of the distribution, meaning that the overlap problem occurs around there. We conclude that the result on the 36 3 × 6 lattice is not reliable at K ≈ 0.119. This explains the strange behavior of the histogram in Fig. 6 at K = 0.119. To obtain a reliable V eff on this lattice, we would need full QCD simulations or quenched simulations with an external source term of the Polyakov loop [32]. The result of V eff (|Ω|) obtained on the 24 3 × 8 lattice is shown in Fig. 12. We find that the double-well shape of V eff is quite stable up to K = 0.180. At the same time, from the (P, |Ω|) distribution shown in Fig. 13, we find that these peak positions of the histogram for K = 0.0, 0.15 and 0.18 are in the area where the number of configurations is enough. Therefore, the result of the  Because this value of K ct is quite large, we examine the location of the chiral limit. We perform zero-temperature simulations of two-flavor full QCD at β = 5.992, which corresponds approximately to the transition point at these K's by the reweighting study using quenched configurations. Details of the study on the zero-temperature lattice are given in Sec. VI. In Fig. 14, we plot the results of the pseudo-scalar meson mass squared (m PS a) 2 by cross symbols as function of 1/K. Because they align on a straight line we perform a linear fit to obtain the chiral limit around K ≈ 0. 15. This implies that the hopping parameter expansion is not applicable at K > ∼ 0.15. The large value of K ct for N t = 8 is obtained by an analysis with the LO hopping parameter expansion. The problem must be due to the truncation error of the hopping parameter expansion. We discuss the issue in Sec. V. TABLE III: Summary of K ct and K ct,eff for two-flavor QCD using the leading order (LO) or up to the next-to-leading order (NLO) terms of ln det M . K ct on the 24 3 × 4 lattice was obtained in Ref. [23].
lattice K ct up to LO K ct up to NLO K ct,eff with effective NLO 24 3 × 4 0.0658(3)( +4 −11 ) -0.0640(10) 24 3 × 6 0.1359 (30) 0.1202 (19) 0.1205 (23)  Our results of K ct in two-flavor QCD using the LO terms of ln det M are summarized in Table III. Results with terms up to the NLO are discussed in the next section. If N t = 4 and 6 lattices are both in the asymptotic scaling region, we would expect that K ct for N t = 6 is 3/2 times larger than that for N t = 4 because K ∼ (m q a) −1 = N t T c /m q in the limit m q a 1. However, K ct for N t = 6 turned out to be about twice as large as that for N t = 4. Though the central values of K ct for N t = 6 may be indicating decrease of K ct with increasing the spatial volume, our data are yet not sufficient to take a large volume extrapolation. The spatial volume dependence will be studied in more detail in a separate paper [32].

V. INFLUENCE OF NEXT-TO-LEADING ORDER TERMS
We now study the influence of the NLO terms of the hopping parameter expansion in the reweighting factor ln det M given by Eq. (14). As discussed in Sec. II B, we disregard the effects of the 6-step Wilson loops in this study. We thus concentrate on how the O(K Nt+2 ) bended Polyakov loops Ω n with N t + 2 steps, as shown in Fig. 2 for the case of N t = 6, affect the phase structure computed with the LO Polyakov loop Ω only.
The NLO results of V eff (|Ω|) and ∆V eff obtained on the 24 3 × 6 lattice are shown in Figs. 15 and 16. The critical point extracted by a linear fit using the smallest four data points of ∆V eff turns out to be K ct = 0.1202 (19).

A. Effective NLO method
Before proceeding, let us discuss a method, introduced in Ref. [24], to effectively incorporate NLO terms in the LO calculation of V eff and K ct . The method is based on the strong linear correlation between Ω and Ω n observed on N t = 4 lattices [24]. In Fig. 17, we show the distribution of (ReΩ, ReΩ n ) measured on each configuration of the 24 3 × 6 lattice near the transition point. Red, green, and blue dots are for ReΩ 1 , ReΩ 2 , and ReΩ 3 , respectively. We find that ReΩ and ReΩ n have strong linear correlation with each other also on the N t = 6 lattice.
With the strong linear correlation, we may introduce an approximation that ReΩ n ≈ c n ReΩ, n = 1, 2, · · · , N t /2, where c n = ReΩ n /ReΩ . With this approximation, the O(K Nt ) and O(K Nt+2 ) terms in Eq. (14) read with This means that the NLO effects can be reproduced by a shift of the hopping parameter in front of ReΩ in the LO calculation. Our results for ReΩ n /ReΩ are summarized in Table IV for 24 3 × 4 and 24 3 × 6 lattices. In particular, the critical point K ct,LO obtained by the LO calculation can be effectively translated to the critical point K ct,eff to the NLO accuracy by Solving this relation with the K ct,LO obtained on the 24 3 × 6 lattice, we find K ct,eff = 0.1205 (23). This is well consistent with K ct = 0.1202 (19) obtained with the direct NLO calculation. We thus find that the effective NLO method works well. The method is useful to avoid repeating similar analyses, e.g., for various number of flavors. We also improve the calculation of K ct,eff for N t = 4 in Ref. [24], by taking the slight ndependence of c n into account. We obtain K ct,eff = 0.0640(10) on the 24 3 × 4 lattice.

B. Critical point with NLO contributions
Our results of K ct are summarized in Table III. On the 24 3 × 4 lattice, we obtain K ct,eff = 0.0640(10) by the effective NLO method. This is 3% smaller than K ct,LO = 0.0658(3)( +4 −11 ) computed with the LO contribution only, but the difference, i.e., the truncation error of the hopping parameter expansion of ln det M , is small in comparison with the statistical errors.
On the other hand, the truncation error turned out to be not negligible for N t = 6. With NLO contributions, we find K ct = 0.1202 (19), which is about 10% smaller than the LO value 0.1359 (30). As expected, the NLO terms affect to reduce K ct . However, the reduction is not sufficient to achieve K ct ≈ 0.1 expected from the naive scaling with the N t = 4 result.
We reserve a study of NLO effects on N t = 8 lattices for future. The effective method, with which calculation of bended Polyakov loops can be in part reduced, may be useful on a lattice with larger N t in which wider variety of bended Polyakov loops have to be calculated.

VI. MESON MASS AT THE CRITICAL POINT
To make clearer the physical implication of the values for the critical hopping parameter K ct calculated in previous sections, we calculate the pseudo-scalar meson mass m PS corresponding to K ct by performing additional zero-temperature simulations. In this study, we perform the following two simulations: One is the direct two-flavor full QCD simulation adopting the same combination of gauge and quark actions and adjusting the simulation parameters (β, K) to the critical point obtained in the finite-temperature study. Another is the quenched QCD simulation combined with the reweighting method, as adopted in the determination of K ct at finite temperatures. In the latter approach, the LO hopping parameter expansion is adopted, though the influence of Ω is quite small because Ω ≈ 0 on zero-temperature configurations. As discussed in Sec. II B, the effect of the plaquette term can be absorbed by the shift β → β * with β * given by Eq. (17).
Our simulation parameters are summarized in Table V. The combination (β, K) is used in the two-flavor full QCD simulations, while β * is used in the quenched QCD simulations. The hopping parameter K is adjusted to the critical point K ct obtained by finite-temperature simulations on 24 3 ×4 and 24 3 ×6 lattices using the LO and the NLO reweighting factor. In both full and quenched simulations, we generate configurations by the hybrid Monte Carlo algorithm on 16 3 × 32 lattices. The number of configurations is 52 for each simulation point. The meson correlation functions are computed every 10 trajectories after thermalization using a point quark source. We also perform additional simulations at a few points of K around the simulation points given in Table V to evaluate d(m PS a)/dK. The information of the derivative is used for the error estimation of m PS due to the ambiguity of K ct , and also for the calculation of m PS at K ct,eff which is slightly different from the simulation point K ct .
To evaluate m PS a, we fit the pseudo-scalar meson correlation function π(t)π(0) by the following fit function, where m PS a and C are the fit parameters. The fit range is decided by consulting the effective mass defined by This effective mass is constant with m PS a when the correlation function is well approximated by Eq. (22). We choose the fit range where the effective mass shows a plateau.
Our results for the effective mass are shown in Fig. 18 as function of t. The left panels are obtained by the quenched simulation with reweighting, and the right panels are by two-flavor full QCD simulations. The top two plots are obtained at K ct of the LO for the 24 3 × 4 lattice. The two plots in the second row are obtained at K ct of the NLO for the 24 3 × 4 lattice, those in the third row are at K ct of the LO for the 24 3 × 6 lattice, and those in the bottom row are at K ct of the NLO for the 24 3 × 6 lattice. The horizontal solid lines represent the fit range and simultaneously the results of m PS a with their statistical errors shown by the dashed lines.
The results of the pseudo-scalar meson mass are summarized in Table VI. The errors for m PS a in Table VI represent the statistical errors from the hadron mass fit. On the other hand, in the errors for m PS /T c = N t × (m PS a), we include those propagated from the error of K ct using d(m PS a)/dK obtained by additional simulations near K ct .
Comparing the results of two-flavor full QCD simulation with those of the quenched simulation with reweighting, we find that the results for N t = 4 are well consistent with each other. Though the results of m PS a show slight deviation for the case of N t = 6, the deviation is much smaller than the error propagated from the error of K ct , as shown by m PS /T c . Recall that the 6-step Wilson loops are neglected in the reweighting calculation, while the quark determinant is fully included in the full QCD simulation. Hence, the consistency of reweighting and full QCD results means that the truncation error of the hopping parameter expansion is small in the calculation of m PS at the critical points on N t = 4 and 6 lattices. Because the computational cost of the quenched simulation is much smaller than that of the full QCD, we may compute other hadron masses at K ct by the quenched simulation with reweighting.
On the other hand, the truncation error in the calculation of K ct is large for the N t = 6 lattice as discussed in Sec. V. In accordance with this, the value of m PS /T c at the LO K ct is about 1.4 times larger than that at the NLO K ct for N t = 6. If we assume that the systematic error from the truncation of the hopping parameter expansion is roughly the same as the difference between the LO and NLO results, we obtain m PS /T c = 11.15(42)(372) for N t = 6. This is smaller than the result m PS /T c = 15.73 (14) (26) for N t = 4. Hence, our results suggest that the critical pseudo-scalar meson mass measured on the N t = 6 lattice is smaller than that on the N t = 4 lattice.

VII. CRITICAL LINE IN 2+1 FLAVOR QCD
We extend the analysis of the critical point to 2 + 1 flavor QCD. Up to the LO of the hopping parameter expansion, the quark determinant term is given by where K ud and K s are the hopping parameters of the degenerate up and down quarks and of the strange quark, respectively. The case of two-flavor QCD is reproduced by setting K s = 0 in Eq. (24). As discussed in Sec. II B, the effects of the plaquette term can be absorbed by a shift of the gauge coupling β. We then find that the reweighting factor for 2 + 1 flavor QCD is obtained by replacing 2K Nt in two-flavor QCD with 2K Nt ud + K Nt s . In particular, the critical line in the (K ud , K s ) plane of N f = 2 + 1 QCD and the LO critical point K ct,LO of N f = 2 QCD as [23] 2K Nt ct,ud + K Nt ct,s = 2K Nt ct,LO .
We can easily extend this relation in the LO to that in the NLO applying the effective NLO method discussed in Sec. V A. By substituting Eq. (21) into Eq. (25), we obtain an effective NLO relation for the critical line, 2K Nt ct,ud 1 + C Ω N t K 2 ct,ud + K Nt ct,s 1 + C Ω N t K 2 ct,s = 2K Nt ct,LO , with C Ω given by Eq. (20). We solve Eqs. (25) and (26)  Finally, we calculate the pseudo-scalar meson mass on the critical line adopting the reweighting method discussed in Sec. VI. In the reweighting method, we shift β to β * = β + 48(2K 4 ud + K 4 s ) for 2 + 1 flavor QCD. Here, β * should be fine-tuned to the transition point depending on the values of K ct,ud and K ct,s . However, because m PS /T c is not so sensitive to β, we adopt β * = 5.680 and 5.870 for all (K ct,ud , K ct,s ) on the N t = 4 and 6 lattices, respectively. The results of m PS a and m PS /T c as functions of (K ct,ud , K ct,s ) are summarized in Table VII. We have checked at several simulation points that m PS /T c does not change within the error even if β * is fine-tuned to the transition point. Our results of m PS /T c are plotted in Fig. 21 (N t = 4) and Fig. 22 (N t = 6) as function of K s . As in the case of two-flavor QCD, we find that m PS /T c on the N t = 6 lattice is smaller than that on the N t = 4 lattice. We also note that the truncation error of the hopping parameter expansion is large for N t = 6.

VIII. CONCLUSIONS
We studied the endpoint of the first-order phase transition in the heavy quark region of QCD, by computing the effective potential V eff (|Ω|) defined as the logarithm of the probability distribution function of the absolute value of the Polyakov loop Ω. We evaluated V eff (|Ω|) at various values of the hopping parameter K using the reweighting method from the heavy quark limit. We evaluated the reweighting factor at small K up to the leading O(K Nt ) as well as the next-to-  leading O(K Nt+2 ) of the hopping parameter expansion for the Polyakov-loop-type operators. To investigate the effective potential in a wide range of |Ω| covering the hot and cold phases at the first-order phase transition, we combine the data obtained at several simulation points using the multi-point reweighting method. We found that the clear double-well shape of V eff (|Ω|) at K = 0 changes gradually as we increase K, and eventually turns into a single-well shape. We define the critical point K ct as the point where the double-well shape disappears.
Our results of K ct for two-flavor QCD are summarized in Table III. We found that, though the volume dependence of K ct is small, K ct seems to decrease as the volume increases. On our largest lattice, 36 3 × 6, however, we could not determine K ct due to the overlap problem. We reserve a study of the volume dependence for a future work [32]. Comparing the results of K ct calculated up to the leading and next-to-leading orders of the hopping parameter expansion, we found that the truncation error of the hopping parameter expansion is not small at N t = 6, while it was negligible at N t = 4. A reason that the convergence of the hopping parameter expansion is worse at N t = 6 will be that K ct at N t = 6 is larger than that at N t = 4, as expected from the fact that 1/K is approximately proportional to the quark mass times the lattice spacing. However, by comparing the results for N t = 4 and 6, we note that K ct at N t = 6 is much larger than that expected from naive scaling with the data at N t = 4. To make implication of large values of K ct clearer, we also calculated the pseudo-scalar meson mass m PS at zero-temperature just on the K ct for N t = 4 and N t = 6. Although the truncation error of the hopping parameter expansion is large for N t = 6, our results suggest that the critical pseudo-scalar meson mass measured on N t = 6 lattice is smaller than that on N t = 4 lattice. This means that the critical quark mass decreases as the lattice spacing decreases.
We also extend the study of two-flavor QCD to 2 + 1 flavor QCD and determine the critical line at the boundary of the first-order transition region in the (K ud , K s ) plane. Our results for the critical line in 2 + 1 flavor QCD are summarized in Figs. 19 and 20 for N t = 4 and 6, respectively. The pseudo-scalar meson mass on the critical line are shown in Figs. 21 and 22. General tendencies are the same as the two-flavor case.
To extract quantitative physical predictions, we need to reduce the lattice spacing by increasing N t . Our test study on an N t = 8 lattice however suggests that the reweighting from quenched QCD using the low order hopping parameter expansion is not applicable at N t > 6. Careful treatments including higher order terms of the hopping parameter expansion are required. One possible approach is to take into account the effects of further higher order terms of the hopping parameter expansion (12). Another approach is to combine with full QCD simulations. We leave them as future studies.
technique as w(P, |Ω|; β i , K) = e 6N site (β i −β)P w(P, |Ω|; β, K). (A2) Let us first consider the case to combine information obtained at several simulation points, (β i , K) with i = 1, 2, . . . , N sim , at a fixed K. The naive histogram using all the configurations, disregarding the differences in the simulation parameters, is given by where N i and Z i = Z(β i , K) are the number of configurations and the partition function at the i-th simulation point (β i , K), respectively. From this relation, we obtain a useful expression for the histogram at (β, K) as with β = (β 1 , · · · , β N sim ) and For notational simplicity, we suppress the dependence on K in G. . (A6) The right-hand side is just the naive summation of G(P ; β, β) observed on all the configurations disregarding the differences in the simulation parameters. The partition function at (β i , K) can be determined, up to an overall factor, by the consistency relations, for i = 1, · · · , N sim . Denoting f i = − ln Z i , these equations can be rewritten as , i = 1, · · · , N sim .
Starting from appropriate initial values of f i , we solve these equations numerically by an iterative method. Note that, in these calculations, one of the f i 's must be fixed to remove the ambiguity corresponding to the undetermined overall factor. The histogram of |Ω| is obtained by the integral of w(P, |Ω|; β, K) in terms of P , The expectation value of an operatorÔ[P , |Ω|] at (β, K) can be evaluated as whereÔ[P , |Ω|] is an operator which is written by (P , |Ω|), and Z(β, K) is given by Eq. (A6). When one want to change also the hopping parameters, from K 0 to K, the probability distribution function can be evaluated as where G's in the right-hand side are defined at K 0 . Again, N sim i=1 N i · · · G (β i ,K 0 ) in the right-hand side is just the naive sum of · · · G over all the configurations disregarding the differences in the simulation parameters. The expression for Ô [P , |Ω|] (β,K) is given similarly. Figure 23 shows an example of the histogram w(P, |Ω|; β) in the heavy quark limit K = 0 just at the transition point of the pure gauge theory, β = 5.69121, obtained on a 24 3 × 4 lattice. Using configurations of Ref. [23], data at five different simulation points are combined by the multi-point histogram method. We see a double-peak distribution at K = 0. The values of P and |Ω| of each configuration are plotted in Fig. 24 -red, green blue purple and cyan dots are results at β = 5.68, 5.685, 5.69, 5.6925 and 5.70, respectively. The yellow circles are the expectation values calculated by the multi-point histogram method varying β from 5.66 to 5.76, using the data at β = 5.68-5.70. To test the reliable range of the multi-point histogram method, we made additional simulations at β = 5.665 and 5.76. The black circles are the expectation values of P and |Ω| directly calculated without the reweighting method at β = 5.665 and 5.76 besides the above-mentioned five β's. We find that the results of direct simulations at β = 5.665 and 5.76 deviate from the results of reweighting using the data at β = 5.68-5.70. These data points are in the regions where the distribution of P and |Ω| used in the reweighting is not overlapping with the target expectation values. Because the reweighting method only changes the weight of each configuration in the average over the configurations, the method fails when the important region of relevant operators is out of the original distribution.