Nature of the phase transition for finite temperature $N_{\rm f}=3$ QCD with nonperturbatively O($a$) improved Wilson fermions at $N_{\rm t}=12$

We study the nature of the finite temperature phase transition for three-flavor QCD. In particular we investigate the location of the critical endpoint along the three flavor symmetric line in the light quark mass region of the Columbia plot. In the study, the Iwasaki gauge action and the nonperturvatively O($a$) improved Wilson-Clover fermion action are employed. We newly generate data at $N_{\rm t}=12$ and set an upper bound of the critical pseudoscalar meson mass in the continuum limit $m_{\rm PS,E}\lesssim 110$MeV.


I. INTRODUCTION
The finite temperature transition in QCD is an important subject in elementary particle physics and cosmology. The nature of the finite temperature phase transition has been studied over a number of years. So far there are some analytic attempts to investigate the nature of the phase transition; effective theories based on the universality argument [1][2][3][4] were systematically studied and recently an anomaly matching argument [5,6] has been developed. Although these approaches can capture qualitative aspects, it is hard to provide its quantitative information on the nature of the phase transition without fully taking the nonperturbative effects of QCD. Lattice QCD simulations play central roles in revealing the quantitative aspects, and in fact many efforts have been devoted for this aim. See reviews [7][8][9][10] for a current status of the QCD phase structure with the finite temperature and quark number density. In such studies, the so-called Columbia plot [11] is often used to express the nature of the phase transition in various parameter space. A whole structure of the plot is basically dictated by a critical point, line or surface, which separates the first order and crossover region, depending on the dimensionality of the parameter space. It is therefore crucial to figure out the shape of such critical boundaries. The standard Columbia plot in the case of zero density has two axes: the up-down and strange quark masses (See Fig. 1). In the heavy mass region of the plot, especially the static limit is well established as the first order phase transition [12,13] and the heavy region apart from the static limit is also studied [14,15]. On the other hand, the light quark mass region is still under debate and we will closely investigate such region in the following. Although there are interesting issues for two-flavor QCD, for example restoration of the U A (1) symmetry and so on (see Refs. [7][8][9][10] for recent progress), in this paper we restrict ourselves to the three-flavor symmetric case where all three quark masses are degenerated. In particular, our goal is to locate the critical endpoint along the flavor-symmetric line on the standard Columbia plot. Let us look back on a historical background of the location of the critical endpoint for the three-flavor QCD. A rough but first estimate of the critical endpoint was provided by Iwasaki et al., [16] using the Wilson-type fermions and their critical quark mass is relatively heavy m q,E 140 MeV, equivalently m PS,E O(1) GeV in terms of the pseudoscalar meson mass. Subsequently a study with the standard staggered fermion action was carried out by JLQCD collaboration [17] and they estimated am q,E ≈ 0.03 using screening mass analysis. A similar study was done by Liao [18] and similar conclusion was drawn. Then Karsch et al., [19] reported m PS,E ≈ 290 MeV using the Binder intersection method with the combination of the standard staggered fermions and the Wilson plaquette gauge action, and in addition they also estimated m PS,E ≈ 190 MeV using an improved staggered-type fermion action, i.e., p4-action. In Ref. [20], they updated m PS,E = 67 (18) MeV for the p4-action with replacing the gauge action to the Symanzik-improved one, then a large cutoff effect on the critical endpoint was indicated. Although the Ralgorithm [21] was used in the staggered fermion studies mentioned above, de Forcrand and Philipsen [22] performed rational hybird Monte Carlo (RHMC) simulation [23,24] and found am q,E = 0.0260(5) which is significantly smaller than the previous value with the R-algorithm. Smith and Schmidt [25] examined the RHMC results using larger spatial volumes and it was confirmed that the critical point belongs to the three-dimensional Z 2 universality class. In the above staggered studies, a single lattice spacing was exclusively used (the temporal lattice size was fixed to be N t = 4), however, de Forcrand et al., [26] extended their study to see the lattice cutoff dependence and found that m PS,E /T E , T E is temperature at the critical endpoint, decreases from 1.680(4) to 0.954 (12) as increasing N t from 4 to 6. This explicitly shows that it is important to control the cutoff effects on the critical point and also suggests that the critical mass in the continuum limit may be quite small. Further studies were continued using the improved staggered fermions with smearing techniques [27][28][29][30], but they could not even detect a critical point. Instead, for example, Ding et al., [30] quoted an upper bound of the critical mass m PS,E 50 MeV.
In the above situation, we embarked a study of the nature of the phase transition using the O(a) improved Wilson fermions instead of the staggered-type fermions. Such a study is important to check the universality when taking the continuum limit and our formulation is completely free of the rooting issue [31,32]. In the early stage of our study [33] with coarse lattice spacings N t = 4, 6 and a part of N t = 8, we observed a quite large scaling violation in the continuum extrapolation of the critical endpoint. Therefore we extended our study to N t = 8 and 10 [34] together with the multiensemble reweighting technique. Then we confirmed the universality class of the critical endpoint to be Z 2 universality class for N t = 4 and 6, while it is assumed for N t = 8 and 10 and we used a modified fitting form of the Binder (kurtosis) intersection analysis. And then we set an upper bound m PS,E < 170 MeV. In the current paper, we further extend our study and generate the new data set of N t = 12 in order to take the continuum limit smoothly and make sharpe the prediction of the critical point if exists.
The rest of the paper is organized as follows. We describe the simulation setup and the analysis methods in Sec. II. In Sec. III, we locate the critical point by applying two analysis methods for a cross check. Then we discuss the continuum limit of the critical pseudoscalar mass and the critical temperature. Our conclusions are summarized in Sec. IV. Results of zero temperature simulations for scale setting are summarized in Appendix A.

II. SETUP AND METHODS
Our finite temperature N f = 3 QCD simulations are performed with the Iwasaki gauge action [35] and nonperturvatively O(a) improved Wilson-Clover fermion action [36]. In this paper we report our newly generated data with the temporal lattice size of N t = 12. To carry out the finite size scaling analysis 5 different spatial lattice sizes with N s = 16, 20, 24, 28, and 32 are used. As we will see soon, the smaller spatial lattices 16 and 20 are used only when estimating the transition point in the thermodynamic limit and the critical points will be determined using the larger volumes 24, 28, and 32, which satisfy m PS L 4. Gauge configurations are generated with the RHMC algorithm [24] implemented with the Berlin QCD code [37], where the acceptance rate is tuned to be around 70-80%. Observables are measured at every 10th molecular dynamics trajectory whose length is set to unity. There is a single hopping parameter κ for three degenerate dynamical flavors in our simulations, which is adjusted to search for a transition point at each β, where β values are chosen in a range between 1.80 and 1.82. See Table I for the parameter sets and their statistics.
We follow the same analysis methods as our previous studies [33,34,38], which are summarized in the following: 1. The chiral condensate and its higher order moments up to the fourth are measured. The definition of the moments is given in our previous paper [38].
2. The multiensemble reweighting [39] in only κ but not β is used, which enables us to smoothly interpolate the moments. The reweighting factor, which is given by the ratio of fermion determinants at different κ values, is calculated with an expansion of the logarithm of the determinant [38]. Adopting an expansion form for the moments in the reweighting method, we can evaluate the moments at continuously many points at a relatively low cost.
3. From these moments the susceptibility, the skewness, and the kurtosis, which is equivalent to the Binder cumulant up to an additional constant, are calculated.
4. The κ value at the transition point is estimated from the peak position of the susceptibility at each β.

5.
After repeating the procedure 1-4 for a few spatial lattice sizes, the location of the critical point is estimated by the kurtosis intersection analysis [19], where we search for a point at which the kurtosis value for the phase transition is independent of the volume as schematically illustrated in Fig. 2. In the determination of the critical point we use a fit ansatz with the inclusion of the energy-like observable contribution [34]. kurt. kurt.

FIG. 2. An illustration of the kurtosis intersection analysis.
Kt denotes the kurtosis value for the phase transition and E indicates the critical endpoint, where Kt is independent of the volume.

A. Moments and location of the transition point
As an illustration of the data, we show the susceptibility and the kurtosis of the chiral condensate for β = 1.80 and 1.81 in Fig. 3 together with the κ-reweighting results. From the peak position of the susceptibility, we extract the value of κ at transition points denoted as κ t (β, N s ), whose values are summarized in Table II for various N s and β. The peak height of the susceptibility and the minimum of the kurtosis are also shown in the table. For each value of β, the infinite volume limit of the transition point κ t (β, N s = ∞) is carried out by using a fitting form with an inverse spatial volume correction term 1 , where κ t (β, N s = ∞) and c(β) are fitting parameters. For β = 1.80, three smaller volumes N s = 16, 20, and 24 are used while all five volumes are used for β = 1.81 and 1.82. The quality of the fitting is reasonable with for all cases. The resulting phase diagram in the bare parameter space is summarized in Fig. 4. The phase transition line in the thermodynamic limit (we denote κ t (β) = κ t (β, ∞)) at N t = 12 is determined by the linear interpolation,

B. Kurtosis intersection analysis
The minimum of kurtosis is plotted in Fig. 5 to perform kurtosis intersection analysis at N t = 12. The left panel of Fig. 5 includes all N s data points. At β = 1.80, a typical behavior of the first order phase transition is clearly seen; the kurtosis tends to be smaller with increasing the volume. On the other hands the results at β = 1.81 show volume independent behavior. For β = 1.82, apart from N s = 32 data point which has a relatively larger error bar, the crossover behavior is seen in the volume dependence; for larger volume the kurtosis tends to be larger. Therefore it is likely that there is a crossing point between β = 1.80 and 1.82. To keep away from the finite size effects we use N s ≥ 24 data points (m PS L ≥ 4) for the kurtosis intersection fitting. We employ the following fitting form [34] which incorporates the correction term associated with the contribution of the energy-like observable, where K E , β E , A, ν, B, and y t − y h are basically fitting parameters. Following Ref. [34], we assume the threedimensional Z 2 universality class for K E = −1.369, ν = 0.630, and y t − y h = −0.894, namely the actual fitting parameters are now β E , A and B. The fitting results are shown in Table III together with the previous smaller N t ≤ 10 results. The quality of the fitting for N t = 12 is reasonable. In the table, κ E is estimated by an interpolated transition line in Eq.
(2) together with the corresponding β E as an input.   [34]. The open symbols represent a transition point (TP) while the filled symbols are the critical endpoints (CEP) determined by the kurtosis intersection with a correction term. On the transition line, the left (right) hand side of an critical endpoint is the first order phase transition (crossover) side. The phase transition line for each Nt is a polynomial interpolation. For Nt = 12 the interpolation formula is given in Eq. (2). In the plot, κc is the pseudoscalar massless point with N f = 3 which is determined by zero temperature simulations as shown in Table VII in Appendix A.   As seen above, the kurtosis intersection analysis is not fully satisfactory since we have heavily relied on the assumption of the universality class of the critical point. Therefore we should cross check the location and the universality class of the critical point. For that purpose, we investigate the scaling of the susceptibility peak height for the chiral condensate, At a critical point, the exponent should be b = γ/ν, where γ and ν are critical exponents. As discussed in Ref. [34], in general there are correction terms in the above formula but here we neglect them just for a simplicity. The data in Table II is fitted with the above functional form as seen in Fig. 6. The resulting exponent b is plotted in Fig. 7 along the transition line projected on β. Assuming the Z 2 universality class (γ = 1.237 and ν = 0.630) provides an estimation of the critical point of β, we confirm that it is consistent with that of the kurtosis intersection for N t = 12 as well. This cross check assures that our analysis is working well.

D. Estimate of critical mass in continuum limit
For scale setting, we perform the zero temperature simulations, which roughly cover the parameter range of the critical endpoints, and their results are summarized in Appendix A. For example, in Fig. 8, the pseudoscalar meson mass m PS and the Wilson flow scale √ t 0 [40] in lattice units are plotted as a function of κ at β = 1.80 and 1.81. The blue vertical line represents the location of the transition point κ t (β, ∞) in Table II    the interpolation range, the monotonic behavior of data points suggests that such a short extrapolation should be harmless.
The dimensionless combination of the hadronic quantities √ t 0 T , m PS /T , and √ t 0 m PS along the transition line projected on β are plotted in Fig. 9 for N t = 10 and 12. The vertical red line represents the location of the critical point determined by the kurtosis intersection method, and the plot allows us to obtain the hadronic quantities at the critical point. From an interpolation, one can obtain the critical value of the dimensionless quantities for each temporal size N t . The actual numbers are summarized in Table IV.    Figure 10 shows the continuum extrapolation of √ t 0 m PS,E , m PS,E /T E , and √ t 0 T E . As for √ t 0 T E (lower right panel of Fig. 10), even though the new data point at N t = 12 is included, a stable continuum extrapolation is observed and we obtain √ t 0 T E = 0.09943(34) which has no significant difference compared with the previous one √ t 0 T E = 0.09970(37) in Ref. [34]. In terms of the physical units the critical temperature is given by T E = 134(3) MeV, where we have used the Wilson flow scale 1/ √ t 0 = 1.347(30) GeV in Ref. [41] as an input. On the other hand, in Fig. 10, √ t 0 m PS,E and m PS,E /T E show significantly large scaling violation. In the extrapolation procedure, we try some functional forms including up to quadratic correction term and examine the fitting range dependence. As a result, their dependence turns out to be quite large as shown in Fig. 10   program does not work well in our parameter range. Thus to avoid such an extrapolation, in the future we should do large β > 1.90 simulations, that is, very large N t simulations where the O(a) improvement works well. In such a simulation, O(a 2 ) scaling may be seen and an extrapolation to the negative value could be avoided. In any case, here we conservatively quote an upper bound of the critical value √ t 0 m PS,E 0.08, which is taken from the maximum continuum value among all the fits except for (A). In physical units, this bound is m PS,E 110 MeV, which is smaller than our previous estimate (m PS,E 170 MeV) [34]. A Columbia-like plot, whose axes are given by hadron masses, is shown in Fig. 11 to display the current situation of our study. For future references, we mention the continuum extrapolation of m PS,E /T E in Fig. 10 (lower left) where large cutoff dependence is seen as well. Thus we quote an upper bound m PS,E /T E 0.93 obtained with the same criteria as that of √ t 0 m PS,E .

IV. SUMMARY AND OUTLOOK
In this study, we performed the large scale simulations for N t = 12 by using the Wilson-type fermions. This is an extension of our previous works at the smaller temporal size simulations N t ≤ 10 [33,34]. By using the modified formula of the kurtosis intersection analysis, the critical endpoint is determined with assuming 3D Z 2 universality class. The continuum limit for the critical temperature is smoothly taken and we obtain T E = 134(3) MeV which is essentially the same as before. On the other hand, for the critical mass, the continuum extrapolation significantly dominates the systematic error, thus here we conservatively quote upper bounds where we have made the upper bound about 40% smaller than before.
In fact, the studies using the staggered-type fermions suggested much lower bound m PS,E 50 MeV in Ref. [30]. Thus it is likely that the critical mass is so small that modern computers cannot access it directly, or it could be zero. Moreover an insightful result for N f = 4 QCD was reported by de Forcrand and D'Elia [42] where the standard staggered fermions are used to study the critical point. They found large cutoff effects compared with N f = 3 case and the critical mass tends to be zero with decreasing the lattice spacing. A similar tendency is observed even in the Wilson-type fermions by our group [43]. Since there is no rooting issue when the number of flavor is a multiple of 4, the feature that the critical mass is extremely small for multiple-flavor QCD seems to be robust.
Of course, in order to make a quantitative conclusion, one has to carry out large N t simulations or use the improved lattice actions. Another possibility is to invent a new analysis method which is useful to study such a near-zero critical mass.
very hard, since one needs extremely small lattice size Ns 6 for such a low β case in order to keep the physical length scale constant. Of course one can change the constant physics condition such that a low β simulation is feasible with a reasonable lattice size say Ns = 6, but in that case the physical lattice size is larger and then a high β simulation requires quite large lattice size Ns 6 and moreover infrared cutoff thanks to the boundary condition gets weaker.   (27)