Finite-size scaling around the critical point in the heavy quark region of QCD

Finite-size scaling is investigated in detail around the critical point in the heavy-quark region of nonzero temperature QCD. Numerical simulations are performed with large spatial volumes up to the aspect ratio $N_s/N_t=12$ at a fixed lattice spacing with $N_t=4$. We show that the Binder cumulant and the distribution function of the Polyakov loop follow the finite-size scaling in the $Z(2)$ universality class for large spatial volumes with $N_s/N_t \ge 9$, while, for $N_s/N_t \le 8$, the Binder cumulant becomes inconsistent with the $Z(2)$ scaling. To realize the large-volume simulations in the heavy-quark region, we adopt the hopping parameter expansion for the quark determinant: We generate gauge configurations using the leading order action including the Polyakov loop term for $N_t=4$, and incorporate the next-to-leading order effects in the measurements by the multipoint reweighting method. We find that the use of the leading-order configurations is crucially effective in suppressing the overlapping problem in the reweighting and thus reducing the statistical errors.


I. INTRODUCTION
One of the interesting features of the medium described by quantum chromodynamics (QCD) is the existence of phase transitions of various orders. While the finitetemperature QCD transition is an analytic crossover at zero quark chemical potential µ q [1,2], this phase transition is expected to become of first order in dense medium with large µ q [3]. The end point of the first-order transition is called the critical point (CP) at which the transition is of second order. The singularity in thermodynamic observables associated with the second order CP are believed to be useful in detecting the CP in heavyion collision experiments [4,5]. Accordingly, researches called the beam-energy scan are actively performed in experimental facilities all over the world to search for the critical fluctuations around the CP [5][6][7].
The order of the finite temperature QCD transition changes also with variation of quark masses [8]. For the 2 + 1-flavor QCD, it is known that the crossover at the physical quark masses becomes of first order both in the light-and heavy-quark limits; the phase diagram representing this feature is known as the Columbia plot [9,10]. Revealing the nature of the phase transitions with the variation of quark masses is an important subject of QCD at nonzero temperature since it provides us with various insights into the transition at the physical quark masses.
One of the difficulties in these analyses is that observ-ables near the CP are strongly dependent on the spatial volume of the system. The spatial volume dependence is in part described by the finite-size scaling (FSS) [32]. However, the FSS is applicable only for describing the singular part of thermodynamic quantities that dominates over the non-singular part only in the vicinity of the CP for sufficiently large spatial volumes. When the spatial volume is not large enough the FSS of observables is violated due to the contributions of the non-singular part and this makes their analysis based on the FSS problematic. In fact, although the CP of QCD is believed to belong to the three-dimensional Z(2) universality class [8], a clear FSS in this universality class has not been observed in the latest numerical study around the CP in the light quark mass region [23]. In the heavy quark region, on lattices of N t = 6 and 8 with the aspect ratio N s /N t = 6 and 8, the Binder cumulant B 4 is reported to be consistent with the Z(2) FSS using data mostly in the crossover side, while data in the first order side as well as that on N t = 10 lattice show deviation from the Z(2) FSS [31]. These results suggest the necessity to perform numerical analyses with yet larger spatial volumes and with high statistics.
In the present study, we focus on the CP in the heavy quark region and study the behavior of observables by numerical simulations with large spatial volumes corresponding to the aspect ratios up to N s /N t = 12. To carry out analyses on the large spatial volumes with high precision, we fix the temporal lattice extent to be N t = 4 in this study.
We also employ the hopping parameter expansion (HPE) to deal with the quark determinant. In this study, we generate gauge configurations using the leading order (LO) action of the HPE including the Polyakov loop term for N t = 4, and then incorporate the next-toleading order (NLO) effects by a multipoint reweighting method [33]. In our previous study, we generated the configurations in quenched QCD (zero-th order of the HPE) and reweighted them to incorporate LO and NLO effects [27,29,30]. We show that the use of the LO action to generate configurations is quite effective in suppressing the overlapping problem of the reweighting method. This is essential to carry out simulations with large system volumes as studied in the present paper. We also verify the convergence of the HPE by comparing the LO and NLO results.
We perform the Binder cumulant analysis of the Polyakov loop around the CP and find that the numerical results for N s /N t ≥ 9 follow well the FSS in the Z(2) universality class. On the other hand, inclusion of the data at N s /N t ≤ 8 in the analysis gives rise to a statisticallysignificant deviation from the scaling behavior. We further investigate the scaling behavior of the distribution function of the Polyakov loop. We find that the distribution function follows the FSS in the Z(2) universality class for N s /N t ≥ 8. From the deviation pattern of the distribution function for LT = 6 from the Z(2) FSS, we discuss that the violation of the scaling behavior in the Binder cumulant is caused by the deviation in the tails of the distribution for small LT . In this paper, we consider the case of degenerated N f flavors with N f = 1, 2, 3. Generalization of the formalism to non-degenerate cases is straightforward. This paper is organized as follows. In the next section we give a brief review on the FSS. We then explain the setup of our lattice simulation and analyses using the HPE in Sec. III. In Sec. IV, we determine the transition line and perform the Binder cumulant analysis to determine the location of the CP and the critical exponent. In Sec. V, we investigate the FSS of the distribution function of the Polyakov loop. The last section, Sec. VI, is devoted to a summary. In Appendix A, we give definition of cumulants. In Appendix B, we examine the effect of the smearing width used in the calculation of distribution functions. In Appendix C, the HPE of the quark determinant is calculated up to the NLO. In Appendix D, the convergence of the HPE is examined by comparing the Binder cumulants at the LO and the NLO.

II. FINITE-SIZE SCALING
Let us first give a brief review on the FSS and its application to the CP in the heavy-quark region of QCD.
The heavy-quark limit of QCD corresponds to the SU (3) Yang-Mills theory (quenched QCD). This theory has a first-order deconfinement phase transition at nonzero temperature T . When the quark mass m q is finite (throughout this section we assume that the quark masses are degenerate), with increasing 1/m q , this firstorder transition becomes weaker and eventually terminates at the CP, as schematically shown in the phase diagram on the (T, 1/m q ) plane in Fig. 1 (a). This CP, as well as that in the light quark region, is believed to belong to the Z(2) universality class, i.e. the universality class of the three-dimensional Ising model [8].
Near the CP of the three-dimensional Ising model, the relevant scaling parameters are the reduced temperature t and external magnetic field h; extensive variables conjugate to these parameters are the energy and the magnetization, respectively. As shown in Fig. 1 (b), the CP is located at (t, h) = (0, 0) and the first-order transition exists on the t axis for t < 0. The singular part of thermodynamic quantities near the CP is described by the scaling function of t and h. According to the universality, the singular part of thermodynamic quantities near the CP of heavy-quark QCD is described by the same scaling function, where the scaling parameters t and h are encoded into the (T, 1/m q ) plane as schematically shown in Fig. 1 (a); the t axis is parallel to the firstorder line at the CP while the direction of the h axis is not constrained from the universality.
According to the FSS argument [32] the singular part of the dimension-less free energy F (t, h, L −1 ) around the CP obtained at a finite volume V = L 3 has a scaling for arbitrary scale factor b. The values of the exponents y t and y h are specific for each universality class. In the Z(2) universality class these parameters are numerically obtained as [32] y t = 1.588, y h = 2.482.
By setting b = L one has Derivatives of F (t, h, L −1 ) with respect to t and h define the cumulants of the corresponding extensive variables. For example, cumulants of the magnetization M are given by with ∂ h = ∂/∂h. From Eq. (3), L dependence of the cumulants M n c near the CP are written as In Ref. [34], it is suggested that the so-called (fourthorder) Binder cumulant plays a useful role to determine the location of the CP from numerical results obtained at finite L. When the distribution of M obeys the Gauss distribution or the distribution composed of two delta functions with an equal weight, we have respectively. Since the distribution of M approaches these functions in the L → ∞ limit on the crossover and first-order lines at h = 0, respectively, B 4 should approach Eq. (7) on these lines in the L → ∞ limit. Moreover, from Eq. (4), B 4 at h = 0 behaves as a function of t and L as for small t, where b 4 = ∂ 4 hF (0, 0)/(∂ 2 hF (0, 0)) 2 + 3, ν = 1/y t , and c is a constant. Eq. (8) shows that B 4 (t, 0, L −1 ) obtained for various L at h = 0 has a crossing at t = 0. The parameter b 4 is given only fromF (t, h) and thus are specific for each universality class. For the Z(2) universality class, the value is known to be [32] b 4 = 1.604.
When we consider magnetization per unit volume m = M/V , the probability distribution of m is given by At the CP, (t, h) = (0, 0), one finds from Eq. (13) that Equation (13) also suggests that, whenp M (M ; t) has a local extremum at M =M (t), p m (m; t, 0, L −1 ) has corresponding local extremum at This implies that the t and L dependences of the maximum of p m (m; t, 0, L −1 ) are described by a single func-tionM (t).

A. Lattice action and basic observables
In this study we investigate the four-dimensional system described by the lattice action of QCD with S g and S q being the gauge and quark actions. For S g we employ the plaquette action with the gauge coupling parameter β = 6/g 2 and the space-time lattice volume N site = N 3 s ×N t . The plaquette operatorP is given bŷ where U x,µ is the link variable in the µ direction at site x, x+μ is the next site in the µ direction from x, N c = 3, and tr C is the trace over color indices. For S q , we adopt the Wilson quark action with the Wilson quark kernel where x, y represent lattice sites. The color and Diracspinor indices are suppressed for simplicity. κ f is the hopping parameter for the f th flavor. The bare quark mass m f is related to κ f as with the lattice spacing a. The matrix B xy has nonzero values only when lattice sites x and y are located in an adjacent sites. Therefore, this term represents the "hopping" of a quark between adjacent sites. The heavy quark limit m f → ∞ corresponds to κ f → 0.
In the following, we consider degenerated N f flavors with a common hopping parameter κ = κ f corresponding to a common quark mass m q -generalization to non-degenerate cases is straightforward. In this case, the expectation value of a gauge operatorÔ(U ) is calculated as with the partition function Z = DU e −Sg+N f ln detM (κ) .
In the heavy quark limit κ = 0 (m q = ∞), the deconfinement phase transition at nonzero temperature is characterized by the spontaneous symmetry breaking of the global Z(3) center symmetry of the SU (3) gauge symmetry. The most conventional choice for the order parameter of this phase transition is the Polyakov loop where the summation x is over the spatial lattice sites on one time slice. In the heavy quark limit, Ω = 0 below the critical temperature T c , while Ω takes a nonzero value at T > T c . For finite m q , the Z(3) symmetry is explicitly broken by the quark term, and thus Ω becomes non vanishing for all T . Even in this case, when m q is sufficiently large, Ω jumps discontinuously at the first-order transition and thus can be used to detect the first-order transition line and its CP [27,29,30].
In Sec. V, we study the scaling property of the distribution function of the real part of the Polyakov loop defined by In Sec. III E, we also calculate the double distribution function ofP andΩ R defined by In numerical calculation of these distribution functions, because the statistics of the data is finite, we have to replace the delta functions by smeared ones with finite width. A conventional choice for this is the normalized  [27,29]. The width ∆ should be large enough to have statistically meaningful number of data at each point within the width, and simultaneously small enough to resolve the functional shape of the distribution function. Examining the resolution and the statistical error of distribution functions, we adopt ∆ ΩR = 0.002 and ∆ P = 0.0001 in this study. In Appendix B, we confirm that the resulting distribution functions as well as other results of observables discussed in this study are stable under variations of the widths around these values.

B. Hopping parameter expansion
To calculate Eq. (23) around the heavy-quark limit, in the present study we adopt the hopping parameter expansion (HPE) for ln detM (κ): where the matrix B xy is defined by Eq. (21) and Tr is the trace over all indices. In Eq. (28), the contribution at κ = 0 is subtracted for convenience. Since B xy takes nonzero values only for adjacent lattice sites x and y, B n is graphically represented by trajectories of n links [35]. Because of the trace in Eq. (28), non-vanishing contributions are given by closed trajectories. The lowest-order contributions of Eq. (28) start from n = 4, and all contributions for odd n vanishes when N t is even. By writing one finds for N t = 4 that the LO contribution is given by the plaquette and the Polyakov loop as The NLO term S NLO consists of the six-step Wilson loops and bent Polyakov loops as Here,Ŵ rec ,Ŵ chair , andŴ crown represent the six-step Wilson loops of the rectangular, chair, and crown types, respectively, as illustrated in Fig. 2.Ω n are the bent Polyakov loops illustrated in Fig. 3, which run one step in a spatial direction, n steps in the temporal direction and return to the original line. All the Wilson loops and Polyakov-loop-type loops are normalized such that W rec =Ŵ chair =Ŵ crown = 1 andΩ n = 1 in the weak coupling limit, U x,µ = 1. Explicit definitions of these operators as well as the derivation of Eqs. (30) and (31) are given in Appendix C.

C. Numerical implementation with HPE
In this study, we generate the gauge configurations with respect to the action at the LO in the HPE, i.e. with We then perform the measurements at the NLO by incorporating the effect of S NLO by the multipoint reweighting method [29,33,36]. In this subsection we discuss the numerical implementation of these analyses.
In the Monte Carlo simulations of pure gauge theory, thanks to the locality of the action S g , it is possible to adopt the pseudo-heat-bath (PHB) and over-relaxation (OR) algorithms for updating gauge configurations. Focusing on a link variable U µ (x), the dependence of S g on U µ (x) is given by where the staple with U x,−µ = U † x−μ,µ for µ > 0, is graphically shown in Fig. 4 (a). In the PHB and OR algorithms, the link variable U µ (x) is updated according to the probability determined by Eq. (34). The fact that Eq. (34) is represented only by local variables near U µ (x) makes this procedure efficient especially on the memory-distributed parallel computing.
When the LO term, Eq. (30), is included into the action, the contribution of a temporal link U 0 (x) to S g+LO is modified as which is schematically shown in Fig. 4 is given by the product of N t −1 link variables along the temporal direction. The contribution of a spatial link to S g+LO is unchanged from Eq. (34). These results on ∆S g+LO (U µ (x)) suggest that the Monte Carlo updates of U µ (x) can be performed by the PHB and OR efficiently even for S g+LO , provided that the temporal direction is not separated into different parallel nodes and Y 0 (x) can be calculated efficiently. Satisfying this condition is not difficult to attain for largevolume simulations. By taking this advantage, in this study we perform update of gauge fields at the LO by combining PHB and OR 1 . Compared with the puregauge simulation, the increase of the numerical cost to deal with S g+LO in this method is small since the additional multiplications of SU(3) matrices required for an update are only N t − 2 times for the temporal links and the cost to update the spatial links is unchanged.
In the measurement of observables, we include all the contribution of S NLO using the multipoint reweighting method 2 . The expectation value of a gauge observablê O(U ) at the NLO at the parameter set (β, κ) is given from the LO simulations at (β,κ) as with and · LO β,κ is the expectation value taken with the action S g+LO (β, κ). In short, we generate gauge configurations for the LO action S g+LO at several values of (β,κ), and evaluate the expectation values to the NLO at various (β, κ) by the multipoint reweighting method.

D. Simulation parameters
In this study, we perform Monte Carlo simulations with fixed temporal lattice size N t = 4, while the spatial extent N s is changed from 24 to 48. This allows us to perform simulations with large aspect ratio N s /N t = LT up to 12, where L is the lattice size along the spatial direction in physical units. For each N s , the gauge configurations are generated for 3 to 5 sets of (β * , λ) shown in Table. I, which are chosen so that β * is close to the transition line at the LO. All numerical results shown in Secs. IV and V are generated by the multipoint reweighting method from these configurations.
The gauge configurations are updated with the LO action S g+LO using the PHB and OR algorithms as discussed in Sec. III C. Gauge configurations are updated by five OR steps after each PHB step. We measure observables every two sets of the PHB+OR updates, i.e. totally ten OR steps and two PHB steps. For all parameters we have performed 6 × 10 5 measurements in this way. In the following, we set N f = 2 to show the numerical results unless otherwise stated.
In Monte Carlo simulations near a first-order transition, because the transition between the coexisting two phases becomes rare when the spatial volume of the system is large, observables averaged over the two phases tend to have quite long autocorrelations. In Fig. 5 we show the Monte Carlo time history ofΩ R for N s /N t = LT = 10 and 12 at the smallest λ, at which the autocorrelation is the longest. The horizontal axis represents the Monte Carlo time in the unit of measurements. To as a function of τ max for the two Monte Carlo trajectories shown in Fig. 5, where τ represents the Monte Carlo time. We estimate the autocorrelation time from the value of τ int at τ max 4τ int , i.e. the crossing point between τ int and the dashed line in Fig. 6, at which the τ max dependence is well saturated. The values of τ int thus determined on each lattice are listed in Table I.
Throughout this study, we estimate the statistical errors of observables by the jackknife method unless otherwise stated, adopting the binsize of 10, 000 measurements which is sufficiently larger than the estimated autocorrelation lengths. We checked that the statistical errors thus estimated are roughly stable within a variation of the binsize from 5, 000 to 30, 000 measurements.

E. Overlapping problem
Our main objective to use the LO action for configuration generation is to avoid the overlapping problem in the reweighting. In Fig. 7 we show the contour plot for the probability distribution function p(P, Ω R ) defined by Eq. (27), obtained at λ = 0.003 (red) and 0.007 (blue) on 40 3 × 4 lattices adopting ∆ ΩR = 0.002 and ∆ P = 0.0001 for the smearing widths. We checked that p(P, Ω R ) hardly changes under variation of ∆ ΩR and ∆ P around this choice. The solid lines are the contours for the distribution measured on the LO configurations generated with S g+LO . Each contour curve is drawn such that the probability inside the contour is 0.9, 0.7, · · · , and 0.1. In Fig. 7, we also show by the dotted lines the contours of p(P, Ω R ) at the NLO, which is obtained by reweighting the LO data at the same (β * , λ). The meaning of the contours is the same as the solid lines. We find that the deviation of the NLO distribution from the original one at the LO is not significant, suggesting that the effects of the NLO contribution are not large. The large overlap of the LO and NLO distributions ensures that, for observables constructed fromP andΩ R , the NLO results obtained by reweighting the LO data at the same β * and λ are statistically reliable. In the analysis of the CP, the overlapping of the distributions is even more improved after adjusting the parameters to the transition line.
On the other hand, from this figure, we find that the overlapping of the distributions at λ = 0.003 and 0.007 is quite poor -the regions with probability larger than 0.7 are not overlapping at all with each other. This means that, if we were to calculate observables at λ = 0.007 by reweighting data obtained at λ = 0.003, or vice versa, the statistical quality of the results would be quite low. In Refs. [27,29,30], the CP in heavy quark region was investigated on lattices with N s /N t = 4-6, by reweighting from pure gauge configurations, i.e. those obtained at λ = 0. Because the overlapping problem becomes quickly severe as the system volume becomes large, the same strategy is not applicable to the present study in which much larger system volumes up to N s /N t = LT = 12 are simulated. Figure 7 shows that, to incorporate the NLO effects by reweighting, the use of the LO action for configuration generation is sufficiently effective in suppressing the overlapping problem. The smallness of the NLO effects in Fig. 7 further suggests that the effects of dynamical quarks are dominated by the LO term for these values of λ.

A. Transition line
We first determine the location of the transition line that corresponds to h = 0 in terms of the Ising parameters; see Fig. 1. In the coupling parameter space (β * , λ), we denote the transition line as β * tr (λ). In this study, we determine β * tr at each λ adopting the following three conventional choices: Fig. 8, we show the LT dependence of β * tr determined by these definitions for several values of λ. The figure shows that the maximum of Ω 2 R c has a visible LT dependence. On the other hand, the zero point of Ω 3 R c and the minimum of B Ω 4 do not have statisticallysignificant LT dependence for LT ≥ 8. This result shows that the zero point of Ω 3 R c and minimum of B Ω 4 are sufficiently close to the β * tr in the L → ∞ limit in this range of LT . In the following, we employ the minimum of B Ω 4 for the definition of the transition line β * tr for each N t . In Fig. 9, we show the transition line on the (β * , λ) and (β, λ) planes obtained at LT = 10.
In Fig. 10, we show the distribution function of Ω R , Eq. (26), on the transition line for several values of λ, where the delta function in Eq. (26) is smeared by the Gauss function as before with the width ∆ ΩR = 0.002. As discussed in Appendix B, dependence of these results on ∆ ΩR is well suppressed around this ∆ ΩR . The shaded bands represent the statistical errors. At λ = 0.003, we see a clear two-peak structure in p(Ω R ) and find that the peaks become sharper as LT becomes larger. This behavior suggests the first-order phase transition at λ = 0.003. At λ = 0.007, on the other hand, while two peaks are observed for LT ≤ 9, the two peaks cease to exist as LT becomes large. This suggests the crossover transition in the L → ∞ limit at this λ.

B. Binder cumulant
Next, let us determine the position of the CP on the (β, κ) plane. As discussed in Sec. II, it is convenient to employ the Binder cumulant B 4 of Ω R This quantity approaches the known values given in Eq. (7) in the L → ∞ limit depending on the order of the transition. Furthermore, provided that Ω R corresponds In the upper panel of Fig. 11, we show B Ω 4 along the transition line as a function of λ for five values of LT . λ is varied continuously by the multipoint reweighting method. The lower panel is an enlargement of the upper panel around the crossing point. The figure shows that B Ω 4 has a crossing at λ = λ c 0.005 and is an increasing (decreasing) function of LT for λ > λ c (λ < λ c ). The existence of the CP at λ 0.005 is suggested from this result.
To determine λ c and the critical exponent ν quantitatively, we fit the numerical results of B Ω 4 by a fitting function motivated by Eq. (8): where b 4 , λ c , ν, c are the fit parameters. In this study, we can vary λ continuously by the multipoint reweighting method. However, because data at different λ on the same volume are correlated, it is not meaningful to use too many λ values. Using the data at three largest volumes, LT = 12, 10, 9, two or three λ values (6 or 9 data points, respectively) should be sufficient for the four parameter fit of Eq. (42). We thus repeat the fit for several choices of λ values, taking the covariance between data at different λ into account in the calculation of χ 2 .
In Table II, we summarize the results of the fit using the data at three largest volumes, LT = 12, 10, 9, and at λ values listed in the left column of the table. The statistical error in the table is estimated by the standard chi-square analysis. The table shows that the value of χ 2 /dof are smaller than unity in the fits with two λ values, but χ 2 /dof is unacceptably large with three λ values, while all the results for the fitting parameters are well consistent within errors. We choose the fit result for λ = (0.0048, 0.0053) depicted by bold characters in Table II as the central value and include the uncertainty in the fits with two λ values as the systematic error.
We repeat the analysis also with other sets of system volumes. The results of the fits with four and five largest volumes, together with the fit with three largest volumes, are summarized in Table III. For the fit with four and five largest volumes, LT = 12-8 and 12-6, we now have 8 and 10 data points for the four parameter fit with two  Table II. The results with the three largest volumes are shown by black triangles, while those with four and five largest volumes are shown by blue squares and green pentagons, respectively. In the figure, b 4 expected from the Z(2) universality, Eq. (9), is shown by the dashed horizontal line.
From Table III and Fig. 11, we find that, when we adopt the fit with the three largest volumes, LT ≥ 9, the fit result of b 4 is consistent with the Z(2) value within about 1σ. On the other hand, when we include smaller volumes, LT ≥ 8 or LT ≥ 6, b 4 from the fits show statistically significant deviation from the Z(2) value. In Table III, we also summarize the results for the critical exponent ν. From the Z(2) universality class, we expect ν = 1/y t ≈ 0.630. We find that the result of the fit with LT ≥ 9 is consistent with the Z(2) value within the error, while the result of the fit with LT ≥ 6 has a significant deviation from the Z(2) value, though the values of χ 2 /dof are all smaller than unity.
We thus conclude that the FSS in the Z(2) universality class is confirmed when the system volume is large enough, LT ≥ 9 -lattices with LT ≤ 8 are not large enough to apply an FSS analysis for B Ω 4 . The value of λ c thus determined is also shown in Fig. 9.
In Appendix D, we perform the analysis of B Ω 4 at the LO of the HPE and compare the results with those at the NLO discussed in this Section. We find that the LO result for λ c is about 2.6% larger than the NLO value. This small difference suggests that the truncation error of the HPE is well under control at the NLO around λ c .

C. Mixing with energy-like observable
So far, we have performed the analyses of B Ω 4 assuming that Ω R corresponds to the magnetization m = M/V of the Ising model. Although our numerical results thus far are in good agreement with this assumption, a possible mixing with the energy-like observable [14,22] in Ω R is not excluded in general. In this case, the behavior of B Ω   (14)(2) near the CP is modified from Eq. (8) as [22] B Ω 4 (λ, To investigate the effect of this possible mixing, we try fits of B Ω 4 based on Eq. (43). We use the values of B Ω 4 at three λ for the fits to increase the number of data points. We find that the six parameter fits with Eq. (43) with the fitting parameters b 4 , λ c , ν, c, d, y t − y h are quite unstable, suggesting that χ 2 has many local minima. The model space of Eq. (43) would be too large against the data. As a next trial, we perform five parameter fits with Eq. (43) by fixing y t − y h = −0.894. In this case, we find that χ 2 still has many local minima, and χ 2 /dof becomes larger compared with the four parameter fit. It is also found that the value of d is consistent with zero within the error for all trials with the variation of λ values. This suggests that the mixing of the energy-like observable in Ω R is negligible around the CP in the heavy-quark region.
In Table IV, we summarize our final results for the location of the CP, (β c , κ c ). In the table, we also show the results for N f = 1 and 3. We note that the N f dependence of the HPE is trivial at the LO in the sense that N f enters the action Eq. (32) at this order only through the combination λ = 64N c N f κ 4 after the replacement β → β * . Therefore, λ c does not depend on N f . At the LO, this allows us to obtain the value of κ c for various N f from the value of κ c at N f = 2 [29]. Because such a simple scaling is no longer applicable at the NLO, we made individual numerical analyses at N f = 1 and 3. From Table IV, we find that the results of λ c are almost insensitive to N f . This means that the NLO effects on λ c are small.

V. DISTRIBUTION FUNCTION OF ΩR
In this section, we study the scaling behavior of the distribution function p(Ω R ) to further investigate the consistency with the Z(2) universality class around the CP.

A. Scaling of distribution function
Let us first focus on the LT dependence of p(Ω R ) at the CP. In the following, instead of p(Ω R ) itself, we study the effective potential defined from p(Ω R ): as this quantity is more convenient in comparing the results at different LT [27,29]. From Eq. (14), the LT dependence of V (Ω R , λ, LT ) at the CP will be described by a single functionṼ (x) as up to an additive constant, where Ω R is subtracted from Ω R to adjust the center of the distribution. To see if the scaling behavior of Eq. (45) is satisfied, we show in Fig. 12 the effective potential V (Ω R ; λ c , LT ) at the CP obtained at five values of LT , as a function of Ω R − Ω R (LT ) 3−y h , where we set 3 − y h = 0.518. For the figure, we adjust the arbitrary constant term of V (Ω R ; λ c , LT ) such that V (Ω (1) ; λ c , LT ) + V (Ω (2) ; λ c , LT ) = 0, where Ω (1) and Ω (2) (> Ω (1) ) are the values of Ω R at the two local minima of V (Ω R ; λ c , LT ). No further adjustments are made in the figure. The error bands do not include the uncertainty of the additive constant. The lower panel is an enlargement of the region indicated by the dotted rectangle in the upper panel. From  Fig. 12, we find that the numerical results for LT = 8-12 agree almost completely within the errors with the scaling relation Eq. (45). This result nicely supports the FSS in the Z(2) universality class at the CP. From the upper panel of Fig. 12, we note that the effective potential for LT = 6 shows a clear deviation from the results for larger volumes at Ω R Ω (1) and Ω R Ω (2) , while it agrees well with them in the range Ω (1) Ω R Ω (2) . This suggests that the deviation from the Z(2) FSS by lattices with small LT , discussed in Sec. IV B, is due to that in the tails of the distribution p(Ω R ) for small LT .

B. Gap between the two minima
Using Eq. (13), the argument of Sec. V A on the effective potential can be extended away from the CP along the transition line. In this subsection, we study the gap between the two local minima of V (Ω R ; λ c , LT ), According to Eq. (15), this quantity should behave around the CP as provided that p(Ω R ) obeys the FSS. In Fig. 13, we show the λ dependence of Ω (1) and Ω (2) . As seen from Fig. 10, a clear two peak structure of p(Ω R ) disappears when λ exceeds some value depending on LT . Even before the disappearance of the two peaks, identification of local maxima of p(Ω R ) becomes unstable by statistical fluctuations. In Fig. 13, we thus truncate the plots for Ω (1) and Ω (2) at finite λ. The shaded areas in the figure represent statistical errors estimated by the jackknife method, for which we repeat the analysis of Ω (1,2) for p(Ω R ) obtained in each jackknife sample with the smearing width of ∆ ΩR = 0.002. 3 As shown in Appendix B, ∆ ΩR dependence of these results is well sup- 3 We see that the errors in Fig. 13 become occasionally large. We find that this is due to statistical oscillations in the shape of p(Ω R ) around the peak: Though the oscillations are within the statistical errors, the peak position in each jackknife sample can jump discontinuously when oscillation appears just at the peak position as we vary λ. This makes the resulting jackknife error large there. From this observation, we think that these large errors are overestimated. pressed at this ∆ ΩR .
From Fig. 13 we extract ∆Ω as a function of λ. In It is interesting to note that the scaling behavior of ∆Ω is observed even at LT = 6, although the FSS of B Ω 4 is violated already at LT = 8. As discussed in the previous subsection, we may understand this when the violation of the FSS for B Ω 4 is due to the violation in the tails of the distribution function p(Ω R ). As seen in Fig. 12, V (Ω R ) at various volumes agrees well for Ω (1) Ω Ω (2) even for small values of LT . As the higher-order cumulants are sensitive to the whole structure of the distribution, B Ω 4 will be more sensitive to the violation of FSS at the tails of the distribution. On the other hand, ∆Ω is insensitive to them by definition. We also note that the statistical error of ∆Ω is naturally small because it is defined by the peaks of the distribution. Therefore, ∆Ω is useful in seeing the FSS around the CP.
From Eq. (47), ∆Ω should behave linearly as a function of (LT ) y h −3 at the CP λ = λ c . In Fig. 15, we show ∆Ω on the transition line at various values of λ, as a function of (LT ) y h −3 . In the same figure, the dashed lines show linear functions ∆Ω = k(LT ) y h −3 for various values of k. Figure 15 suggests that the linear behavior is realized at λ 0.005, which is consistent with our estimation λ c = 0.00503 (14)(2) from the analysis of B Ω 4 .

C. Discussions
Let us comment on the relation of the present results with those given in Refs. [27,29,30]. In these studies, the CP is defined as the point at which the two peak structure of p(Ω R ) disappears. On lattice with finite LT , this leads to λ which is larger than our value of λ c in the L → ∞ limit. In fact, in Ref. [27] the location of the CP is estimated as κ c = 0.0658(3)( +4 −11 ) for N f = 2 and N t = 4, which is about 10% larger than that given in Table IV. Values of λ c (κ c ) which are smaller than Ref. [30] for each N t are also reported by a recent study of the CP in the heavy quark region on fine lattices (N t = 6-10) using the Binder cumulant method [31]. Though the difference may be removed in the L → ∞ limit, a careful extrapolation will be required. Because the FSS is clearly identified in this study, we think that the extrapolation to the large volume limit is stably performed with the present analysis.
We also note that the latent heat at the deconfinement transition in the SU (3) Yang-Mills theory (κ = 0) has been measured in Ref. [38] recently. It was found that the latent heat becomes larger with increasing the spatial volume. This may be attributed to a remnant of the FSS around the Z(2) CP, like ∆Ω studied in the present study.

VI. CONCLUSIONS
In this paper, we studied the distribution function of the Polyakov loop and its cumulants around the CP in the heavy quark region of QCD. Large volume simulations up to the aspect ratio N s /N t = LT = 12 have been carried out to see the finite-size scaling, while the lattice spacing is fixed to N t = 4. We have performed the measurement of observables using the hopping parameter expansion for the quark determinant; the measurement has been performed at the next-to-leading order of the hopping parameter expansion by the multipoint reweighting method evaluated on the gauge configurations generated for the action at the leading order. We found that this analysis is quite effective in reducing statistical errors by avoiding the overlapping problem of the reweighting method, while the numerical cost hardly changes from the pure Yang-Mills simulations. The convergence of the hopping parameter expansion at the next-to-leading order at the critical point is also verified by the comparison with the leading order result.
Using the data on p(Ω R ) thus obtained, we have performed the Binder cumulant analysis for determining the location of the critical point and evaluating the critical exponent. We found that the critical exponent ν and the value of the Polyakov-loop Binder cumulant B Ω 4 at the critical point is consistent with the Z(2) universality class when LT ≥ 9 data are used for the analysis. On the other hand, statistically-significant deviation from the Z(2) scaling is observed when the data at LT = 8 is included, which suggests that this spatial volume is not large enough to apply the finite-size scaling.
The scaling behavior near the critical point is further studied using the structure of the distribution function of the real part of the Polyakov loop, p(Ω R ). We found that the structure of p(Ω R ) for various LT obeys the finite-size scaling well especially near the peaks of p(Ω R ). We have also proposed the use of the gap ∆Ω between the peaks of p(Ω R ) for the finite-size scaling analysis. We showed that the λ and LT dependence of ∆Ω is in good agreement with the Z(2) scaling over a wide range of λ and LT . On the other hand, the deviation of p(Ω R ) around the tails of the distribution is observed on small lattices, which would give rise to the violation of the finite-size scaling of B Ω 4 in small volumes. as From Eqs. (A3) and (A5), one easily finds that the cumulants are related to the moments as and etc. with δx = x − x . The cumulants are useful in representing properties of p(x) than the moments for various purposes. In particular, in statistical mechanics the cumulants of an extensive variable are extensive variables; see, for example, Ref. [5].
Appendix B: Effect of smearing width for distribution function In Secs. IV and V, we calculate the distribution function p(Ω R ) with smearing the delta function in Eq. (26) by the normalized Gauss function with the width ∆ ΩR . With the statistics of this study, we choose ∆ ΩR = 0.002 from an examination of the statistical error and resolution of p(Ω R ).
In this Appendix, we examine the ∆ ΩR dependence of the numerical results, picking up the middle panel of Fig. 10, Fig. 13 and Fig. 14 as representative results. In Fig. 16, we compare these results with those obtained with ∆ ΩR = 0.001 and 0.004. In Fig. 16, the results with ∆ ΩR = 0.001 and 0.004 are shown by the dashed and dotted lines, while the solid lines show the results with ∆ ΩR = 0.002. From this figure, we find that ∆ ΩR dependence is well small around ∆ ΩR = 0.002 and does not affect our conclusions. In this appendix we derive Eqs. (30) and (31). Throughout this Appendix we assume general value for N t .
As in Eq. (28), the HPE of ln[det M (κ)] is given by Tr[B n ]. Since the matrix B has nonzero contributions only between neighboring lattice sites, Tr[B n ] are graphically represented by the closed trajectories with n links [35]. However, trajectories including "appendices" shown in Fig. 17 do not contribute to the HPE because for such a path the product of the matrix in Dirac space vanishes at the tip of the appendix as (1 − γ µ )(1 + γ µ ) = 0 [35]. With this exception, all possible closed trajectories composed of n links contribute to the HPE at the order κ n . In the following, we calculate their contributions by classifying the trajectories by the shape.
Let us start from the plaquette, i.e. 1 × 1 rectangle, which gives a lowest-order contribution to Eq. (28) at the order κ 4 . The plaquette operatorP in Eq. (18) is defined in such a way that P = 1 in the weak coupling limit with U µ (x) = 1. To satisfy this condition Eq. (18) has a coefficient 1/(N c M plaq ) = 1/18, where M plaq is the number of different plaquettes per lattice site; M plaq = 4 C 2 = 6, which is the number of combinations of axes (µ, ν) at which the plaquettes are located.
The contribution from all plaquettes to Eq. (28) is calculated to be Here, D plaq is the coefficient from the trace in the Dirac space where tr D means the trace over the Dirac indices. The factor 2 in Eq. (C1) comes from two directions for each trajectory, which have to be distinguished in the HPE. The factor 1/n in Eq. (28) is canceled by the number of four starting points of a trajectory; this cancellation occurs for all trajectories. Next, let us consider the Wilson loops of length 6 without windings along the temporal direction. At this order there are three types of trajectories; 1 × 2 rectangle, chair, crown, which are shown in Fig. 2. We define the operators,Ŵ rect ,Ŵ chair ,Ŵ crown , corresponding to these trajectories as W chair = 1 N c M chair N site x µ,ν<ρ,ν =µ =ρ s,t=±1 Re tr C U x,sν U x+sν,µ U † x+μ,sν U x+μ,tρ U † x+tρ,µ U † x,tρ , (C4) Re tr C U x,µ U x+μ,sν U x+μ+sν,tρ U † x+sν+tρ,µ U † x+tρ,sν U † x,tρ , with U x,−µ = U † x−μ,µ and U † x,−µ = U x−μ,µ for µ > 0. Eqs. (C3)-(C5) are defined so that Ŵ rect = Ŵ chair = Ŵ crown = 1 in the weak coupling limit again; to satisfy this condition, the operators are divided by M rect = 12, M chair = 48, and M crown = 16, respectively, corresponding to the number of trajectories per lattice site.
The contribution of these trajectories to the HPE of ln det M (κ) is given by with s =rect, chair, crown. D s is the coefficient from the trace in the Dirac space. For the 1 × 2 rectangle we have and similar manipulations lead to D chair = D crown = −16.
Next we consider the Polyakov-loop type operators having a winding along the temporal direction. The lowest-order contribution among them is the Polyakov loop, i.e. the straight lines of length N t . To calculate its contribution, one needs to pay a special attention to the fact that there is only one independent Polyakov loop for each spatial coordinate on one time slice, not for each lattice site. Therefore, their contributions to the HPE is given by where the factor −1 is to be applied because of the antiperiodic boundary condition of the quark determinant. The real part ofΩ has to be taken after multiplying the factor 2 because two directions of a trajectory are independently taken into account. The factor from the Dirac trace for the Polyakov loop is calculated to be Finally, we consider the contribution of the bent Polyakov loops shown in Fig. 3, whose explicit definitions are given by tr C U x,si U x+sî,4 U † x+4,si U x+4,4 U x+2·4,4 · · · U x+(Nt−1)·4,4 , (C10) tr C U x,si U x+sî,4 U x+sî+4,4 U † x+2·4,si U x+2·4,4 · · · U x+(Nt−1)·4,4 , and so on. The factor M bent = 6 is needed to make Ω n = 1 in the weak coupling limit. From the definition ofΩ n we haveΩ n =Ω Nt−n . Also, when N t is even Ω Nt/2 counts each trajectory twice, and thus its contribution to the HPE should be divided by 2. Bearing these facts in mind, the contribution fromΩ n to the HPE of ln det M (κ) is given by where the last term ReΩ Nt/2 should be omitted for odd N t . The contribution of the Dirac trace is calculated to be D bent = 2 Nt+1 . Accumulating these results gives Eqs. (30) and (31).

Appendix D: Comparison of LO and NLO results
In this appendix, to see the convergence of the HPE at the NLO, we repeat the analyses in Sec. IV B at the LO and compare its results with those at the NLO. In Fig. 18, we show the Binder cumulant B Ω 4 obtained at the LO along the transition line. In this figure, the NLO results of Fig. 11 are also shown by the thin dashed lines. We find that, though the difference between the LO and the NLO results grows as λ becomes larger, the deviation is within a few percent level around the crossing point λ c .
In Fig. 18, we also show the result of the four parameter fit with Eq. (42) using data at the LO on LT = 12, 10 and 9 lattices by the black square. The same result at the NLO is shown by gray circle for comparison. We find that the LO result of the fit is consistent with the NLO result within statistical errors: The values of b 4 = 1.630 (24) and ν = 0.620(47) at the LO are hardly changed from the corresponding NLO values given in Table III. Though the central value of λ c = 0.00516 (15) at the LO is about 2.6% larger than the NLO value, it is consistent with the NLO result within errors, suggesting that the truncation error of the HPE at the NLO is well suppressed in these quantities.
The success of the B Ω 4 fit together with the consistency of the fit results with the Z(2) values suggests that the Z(2) scaling is realized also with the LO action when the system volume is sufficiently large. This is reasonable since the scaling properties near the CP are insensitive to detailed structures of the theory.