Nature of the phase transition for finite temperature N f =3 QCD with nonperturbatively O( a ) improved Wilson fermions at N t =12
Yoshinobu Kuramashi, Yoshifumi Nakamura, Hiroshi Ohno, Shinji Takeda
UUTHEP-744, UTCCS-P-129, KANAZAWA-20-01
Nature of the phase transition for finite temperature N f = 3 QCDwith nonperturbatively O( a ) improved Wilson fermions at N t = 12 Yoshinobu Kuramashi, Yoshifumi Nakamura, Hiroshi Ohno, and Shinji Takeda ∗ Center for Computational Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8577, Japan RIKEN Center for Computational Science, Kobe 650-0047, Japan Institute of Physics, Kanazawa University, Kanazawa 920-1192, Japan (Dated: March 17, 2020)We study the nature of the finite temperature phase transition for three-flavor QCD. In particularwe investigate the location of the critical endpoint along the three flavor symmetric line in thelight quark mass region of the Columbia plot. In the study, the Iwasaki gauge action and thenonperturvatively O( a ) improved Wilson-Clover fermion action are employed. We newly generatedata at N t = 12 and set an upper bound of the critical pseudoscalar meson mass in the continuumlimit m PS , E (cid:46)
110 MeV.
PACS numbers: 11.15.Ha,12.38.Gc ∗ [email protected] a r X i v : . [ h e p - l a t ] M a r 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 aresome analytic attempts to investigate the nature of the phase transition; effective theories based on the universalityargument [1–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 thenature of the phase transition without fully taking the nonperturbative effects of QCD. Lattice QCD simulationsplay central roles in revealing the quantitative aspects, and in fact many efforts have been devoted for this aim. Seereviews [7–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 variousparameter space. A whole structure of the plot is basically dictated by a critical point, line or surface, which separatesthe first order and crossover region, depending on the dimensionality of the parameter space. It is therefore crucial tofigure 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 limitis well established as the first order phase transition [12, 13] and the heavy region apart from the static limit is alsostudied [14, 15]. On the other hand, the light quark mass region is still under debate and we will closely investigatesuch region in the following. Although there are interesting issues for two-flavor QCD, for example restoration of theU A (1) symmetry and so on (see Refs. [7–10] for recent progress), in this paper we restrict ourselves to the three-flavorsymmetric case where all three quark masses are degenerated. In particular, our goal is to locate the critical endpointalong the flavor-symmetric line on the standard Columbia plot. FIG. 1. Columbia plot for N f = 2 + 1 QCD at zero density. Let us look back on a historical background of the location of the critical endpoint for the three-flavor QCD. Arough but first estimate of the critical endpoint was provided by Iwasaki et al. , [16] using the Wilson-type fermionsand their critical quark mass is relatively heavy m q , E (cid:38)
140 MeV, equivalently m PS , E (cid:38) O (1) GeV in terms ofthe pseudoscalar meson mass. Subsequently a study with the standard staggered fermion action was carried outby JLQCD collaboration [17] and they estimated am q , E ≈ .
03 using screening mass analysis. A similar study wasdone by Liao [18] and similar conclusion was drawn. Then Karsch et al. , [19] reported m PS , E ≈
290 MeV using theBinder intersection method with the combination of the standard staggered fermions and the Wilson plaquette gaugeaction, 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 actionto the Symanzik-improved one, then a large cutoff effect on the critical endpoint was indicated. Although the R-algorithm [21] was used in the staggered fermion studies mentioned above, de Forcrand and Philipsen [22] performedrational hybird Monte Carlo (RHMC) simulation [23, 24] and found am q , E = 0 . 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 . . N t from 4to 6. This explicitly shows that it is important to control the cutoff effects on the critical point and also suggeststhat the critical mass in the continuum limit may be quite small. Further studies were continued using the improvedstaggered fermions with smearing techniques [27–30], but they could not even detect a critical point. Instead, forexample, Ding et al. , [30] quoted an upper bound of the critical mass m PS , E (cid:46)
50 MeV.In the above situation, we embarked a study of the nature of the phase transition using the O ( a ) improved Wilsonfermions instead of the staggered-type fermions. Such a study is important to check the universality when takingthe continuum limit and our formulation is completely free of the rooting issue [31, 32]. In the early stage of ourstudy [33] with coarse lattice spacings N t = 4, 6 and a part of N t = 8, we observed a quite large scaling violation inthe continuum extrapolation of the critical endpoint. Therefore we extended our study to N t = 8 and 10 [34] togetherwith the multiensemble reweighting technique. Then we confirmed the universality class of the critical endpoint tobe Z universality class for N t = 4 and 6, while it is assumed for N t = 8 and 10 and we used a modified fitting formof the Binder (kurtosis) intersection analysis. And then we set an upper bound m PS , E <
170 MeV. In the currentpaper, we further extend our study and generate the new data set of N t = 12 in order to take the continuum limitsmoothly 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 thecontinuum limit of the critical pseudoscalar mass and the critical temperature. Our conclusions are summarized inSec. 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 nonpertur-vatively O( a ) improved Wilson-Clover fermion action [36]. In this paper we report our newly generated data withthe 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 whenestimating the transition point in the thermodynamic limit and the critical points will be determined using the largervolumes 24, 28, and 32, which satisfy m PS L (cid:38)
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%. Observablesare measured at every 10th molecular dynamics trajectory whose length is set to unity. There is a single hoppingparameter κ for three degenerate dynamical flavors in our simulations, which is adjusted to search for a transitionpoint at each β , where β values are chosen in a range between 1 .
80 and 1 .
82. See Table I for the parameter sets andtheir 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 themoments is given in our previous paper [38].2. The multiensemble reweighting [39] in only κ but not β is used, which enables us to smoothly interpolate themoments. 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 themoments in the reweighting method, we can evaluate the moments at continuously many points at a relativelylow cost.3. From these moments the susceptibility, the skewness, and the kurtosis, which is equivalent to the Binder cumulantup 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 estimatedby the kurtosis intersection analysis [19], where we search for a point at which the kurtosis value for the phasetransition is independent of the volume as schematically illustrated in Fig. 2. In the determination of the criticalpoint we use a fit ansatz with the inclusion of the energy-like observable contribution [34]. TABLE I. Simulation parameters and the number of configurations. β κ N s = 16 N s = 20 N s = 24 N s = 28 N s = 321 .
80 0 . − − − − . − − − . − − . − − − − . − − − − .
81 0 . − − − − . − − − − . − − − − . − − . − .
82 0 . − − − − . − − − − . − − −
960 14000 . −
600 4210 2510 15500 . − − − − . − − − − E1st order crossover K t light heavy V>V>V suscept. kurt.-2 suscept.kurt. FIG. 2. An illustration of the kurtosis intersection analysis. K t denotes the kurtosis value for the phase transition and E indicates the critical endpoint, where K t is independent of the volume. III. RESULTSA. 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 and1 .
81 in Fig. 3 together with the κ -reweighting results. From the peak position of the susceptibility, we extract thevalue 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 s u sc ep t i b ili t y k u r t o s i s κβ =1.80 N s =16 N s =20 N s =24 s u sc ep t i b ili t y k u r t o s i s κβ =1.81 N s =16 N s =20 N s =24 N s =28 N s =32 FIG. 3. The susceptibility (upper half) and kurtosis (lower half) of the chiral condensate as a function κ for several spatial sizes, N s = 16–32. The left (right) panel is for β = 1 .
80 ( β = 1 . σ band) are plotted. spatial volume correction term , κ t ( β, N s ) = κ t ( β, N s = ∞ ) + c ( β ) /N , (1)where κ t ( β, N s = ∞ ) and c ( β ) are fitting parameters. For β = 1 .
80, three smaller volumes N s = 16, 20, and 24 areused while all five volumes are used for β = 1 .
81 and 1 .
82. The quality of the fitting is reasonable with χ / d . o . f . < . κ t ( β ) = κ t ( β, ∞ )) at N t = 12 is determined by the linear interpolation, κ t ( β ) = 0 . − ( β − . × . . (2) B. Kurtosis intersection analysis
The minimum of kurtosis is plotted in Fig. 5 to perform kurtosis intersection analysis at N t = 12. The left panelof Fig. 5 includes all N s data points. At β = 1 .
80, a typical behavior of the first order phase transition is clearlyseen; the kurtosis tends to be smaller with increasing the volume. On the other hands the results at β = 1 .
81 showvolume 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. Thereforeit 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] whichincorporates the correction term associated with the contribution of the energy-like observable, K = (cid:104) K E + AN /ν s ( β − β E ) (cid:105) (1 + BN y t − y h s ) , (3)where K E , β E , A , ν , B , and y t − y h are basically fitting parameters. Following Ref. [34], we assume the three-dimensional Z universality class for K E = − . ν = 0 . y t − y h = − . β 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 interpolatedtransition line in Eq. (2) together with the corresponding β E as an input. There is no specific meaning in using 1 /N correction form. As seen in Table II, we do not observe a significant N s -dependence on κ t ,thus a choice of extrapolation form seems irrelevant. In fact, we performed a fit with exp( − m PS L ) correction term with fixed m PS = 0 . β and κ value. As a result, the thermodynamic value of κ t is consistentwith that obtained with 1 /N correction form. TABLE II. Summary of the value of κ t : the value of κ at the transition point, χ max : the maximum of susceptibility and K min :the minimum of kurtosis for each N s and β . The thermodynamics limit of κ t is taken with the fitting form in Eq. (1). Theerrors of χ max and K min are estimated by the jackknife analysis. β N s κ t ( β, N s ) χ max K min .
80 16 0 . . . − . . . . − . . . . − . ∞ . .
81 16 0 . . − . . . . − . . . . − . . . . − . . . . − . ∞ . .
82 16 0 . . − . . . − . . . − . . . . − . . . . − . ∞ . κ β κ c N t =6 TPN t =8 TPN t =10 TPN t =12 TPN t =6 CEPN t =8 CEPN t =10 CEPN t =12 CEP FIG. 4. Phase diagram for the bare parameter space ( β, κ ) at N t = 12 together with N t = 6, 8, and 10 results in Ref. [34].The open symbols represent a transition point (TP) while the filled symbols are the critical endpoints (CEP) determined bythe kurtosis intersection with a correction term. On the transition line, the left (right) hand side of an critical endpoint is thefirst order phase transition (crossover) side. The phase transition line for each N t is a polynomial interpolation. For N t = 12the interpolation formula is given in Eq. (2). In the plot, κ c is the pseudoscalar massless point with N f = 3 which is determinedby zero temperature simulations as shown in Table VII in Appendix A. -2-1.5-1-0.5 0 1.79 1.795 1.8 1.805 1.81 1.815 1.82 1.825 1.83 k u r t o s i s β all dataN s =16N s =20N s =24N s =28N s =323D Z -2-1.5-1-0.5 0 1.79 1.795 1.8 1.805 1.81 1.815 1.82 1.825 1.83 k u r t o s i s β FittingN s =24N s =28N s =32CEP3D Z FIG. 5. Kurtosis intersection for the chiral condensate at N t = 12. The left panel contains all data points of N s = 16 −
32. Theright panel includes larger volumes N s ≥
24 together with the fitting function in Eq. (3) which assumes the 3D Z universalityclass and contains a correction term. The black pentagon represents the resulting critical value of β .TABLE III. Fit results for kurtosis intersection with the fitting form in Eq. (3) for N t = 4 −
12. In the fitting we assume the3D Z universality class, namely K E = − . ν = 0 .
630 and y t − y h = − .
894 are fixed in the fitting procedure. Using thecentral value of β E as an input, κ E is obtained from the interpolation formula of the transition line in Eq. (2). The error of κ E contains that from only the interpolation procedure but not the error of β E . N t β E κ E A B χ / d . o . f . . . . . .
776 1 . . . − . .
708 1 . . . − . . . . . − . . . . . − . .
3) 1 . C. Analysis for exponent of susceptibility peak height
As seen above, the kurtosis intersection analysis is not fully satisfactory since we have heavily relied on the assump-tion of the universality class of the critical point. Therefore we should cross check the location and the universalityclass of the critical point. For that purpose, we investigate the scaling of the susceptibility peak height for the chiralcondensate, χ max ∝ ( N s ) b . (4)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 datain Table II is fitted with the above functional form as seen in Fig. 6. The resulting exponent b is plotted in Fig. 7along the transition line projected on β . Assuming the Z universality class ( γ = 1 .
237 and ν = 0 . β , we confirm that it is consistent with that of the kurtosis intersection for N t = 12as 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 thecritical endpoints, and their results are summarized in Appendix A. For example, in Fig. 8, the pseudoscalar mesonmass m PS and the Wilson flow scale √ t [40] in lattice units are plotted as a function of κ at β = 1 .
80 and 1 . κ t ( β, ∞ ) in Table II for corresponding β . Fromthis figure, we obtain the hadronic quantity at the transition point. Although the transition point is slightly out of
10 100 16 24 32 χ m a x N s β =1.80 β =1.81 β =1.82 FIG. 6. The volume scaling of susceptibility peak height for N t = 12. Both axes are scaled logarithmically. The filled symbolsare included in the fit but open ones are not b β N t =4N t =6N t =8N t =10N t =123D Z FIG. 7. Exponent of the susceptibility peak height along the transition line projected on β value for N t = 4, 6, 8, 10, and12. The line connecting the data points is to guide readers’ eyes. The point where the line for each N t intersects the (green)horizontal line is an estimate of the critical point assuming the Z universality class. On the other hand, the shaded areasrepresent the critical β determined by the kurtosis intersection analysis. the interpolation range, the monotonic behavior of data points suggests that such a short extrapolation should beharmless.The dimensionless combination of the hadronic quantities √ t T , m PS /T , and √ t m PS along the transition lineprojected on β are plotted in Fig. 9 for N t = 10 and 12. The vertical red line represents the location of the criticalpoint determined by the kurtosis intersection method, and the plot allows us to obtain the hadronic quantities atthe critical point. From an interpolation, one can obtain the critical value of the dimensionless quantities for eachtemporal size N t . The actual numbers are summarized in Table IV. ( 𝑎 𝑚 PS ) 𝑡 ⁄ 𝑎 PS ) 𝑡 ⁄𝑎𝜅 𝑡 ( 𝑎 𝑚 PS ) 𝑡 ⁄ 𝑎 PS ) 𝑡 ⁄𝑎𝜅 𝑡 FIG. 8. The hadronic quantities in lattice units ( am PS ) , √ t /a as a function of κ at β = 1 .
80 (left) and 1 .
81 (right). Thevertical blue line shows the location of the transition point for κ at the corresponding β with N t = 12.TABLE IV. The hadronic dimensionless quantities at the critical endpoint for N t = 4, 6, 8, 10, and 12. Note that N t = 10results are updated compared with the previous work [34] since the hadronic quantities at β = 1 .
78 are updated as shown inTable VII. N t √ t m PS , E √ t T E m PS , E /T E √ t m PS , E χ / d . o . f . √ t T E χ / d . o . f . m PS , E /T E χ / d . o . f . A (fit) a + a /N N t = 8-12 0.1243(52) 12 .
07 0.09943(34) 0 .
63 1.491(50) 13 .
64B (fit) a + a /N t + a /N N t = 6-12 − .
18 – – − .
72C (solve) a + a /N N t = 10-12 0.061(19) – – – 0.71(22) –D (solve) a + a /N t + a /N N t = 8-12 − − Figure 10 shows the continuum extrapolation of √ t m PS , E , m PS , E /T E , and √ t T E . As for √ t T E (lower right panelof Fig. 10), even though the new data point at N t = 12 is included, a stable continuum extrapolation is observed and weobtain √ t T E = 0 . √ t T E = 0 . T E = 134(3) MeV, where we have usedthe Wilson flow scale 1 / √ t = 1 . √ t m PS , E and m PS , E /T E show significantly large scaling violation. In the extrapolation procedure, we try some functional formsincluding up to quadratic correction term and examine the fitting range dependence. As a result, their dependenceturns out to be quite large as shown in Fig. 10 (upper left for √ t m PS , E and lower left for m PS , E /T E ) and Table V. Thefitting named as (A) is not acceptable since the χ / d . o . f . is very large. For other cases, χ / d . o . f . is reasonable but insome cases, the extrapolated mass is negative. Furthermore we plot √ t m PS , E as a function of 1 /N t in Fig. 10 (upperright), which shows a linear scaling behavior thus the leading scaling violation seems O( a ) for this quantity. Sinceour value of c sw around β ∼ .
81 is out of the interpolation range ( β ≥ . a ) improvement In ref [36], the constant physics condition is used to determine c sw . Actually this condition makes the determination of c sw at low β 𝑡 𝑇 t 𝑚 PS , t ⁄ 𝑇 t 𝑡 𝑚 PS , t 𝛽 𝑡 𝑇 t 𝑚 PS , t ⁄ 𝑇 t 𝑡 𝑚 PS , t 𝛽 FIG. 9. √ t T , m PS /T , and √ t m PS along the transition line projected on β for N t = 10 (left) and 12 (right). The verticalred line shows the location of the critical value of β determined by the kurtosis intersection analysis. program does not work well in our parameter range. Thus to avoid such an extrapolation, in the future we shoulddo large β > .
90 simulations, that is, very large N t simulations where the O( a ) improvement works well. In such asimulation, O( a ) scaling may be seen and an extrapolation to the negative value could be avoided. In any case, herewe conservatively quote an upper bound of the critical value √ t m PS , E (cid:46) .
08, which is taken from the maximumcontinuum value among all the fits except for (A). In physical units, this bound is m PS , E (cid:46)
110 MeV, which is smallerthan our previous estimate ( m PS , E (cid:46)
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 continuumextrapolation of m PS , E /T E in Fig. 10 (lower left) where large cutoff dependence is seen as well. Thus we quote anupper bound m PS , E /T E (cid:46) .
93 obtained with the same criteria as that of √ t 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 anextension of our previous works at the smaller temporal size simulations N t ≤
10 [33, 34]. By using the modifiedformula of the kurtosis intersection analysis, the critical endpoint is determined with assuming 3D Z universalityclass. The continuum limit for the critical temperature is smoothly taken and we obtain T E = 134(3) MeV which isessentially the same as before. On the other hand, for the critical mass, the continuum extrapolation significantlydominates the systematic error, thus here we conservatively quote upper bounds m PS , E (cid:46)
110 MeV , m PS , E /T E (cid:46) . , (5)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 (cid:46)
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 bezero. Moreover an insightful result for N f = 4 QCD was reported by de Forcrand and D’Elia [42] where the standardstaggered fermions are used to study the critical point. They found large cutoff effects compared with N f = 3 caseand the critical mass tends to be zero with decreasing the lattice spacing. A similar tendency is observed even in theWilson-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 improvedlattice actions. Another possibility is to invent a new analysis method which is useful to study such a near-zero criticalmass. very hard, since one needs extremely small lattice size N s (cid:28) β 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 N s = 6,but in that case the physical lattice size is larger and then a high β simulation requires quite large lattice size N s (cid:29) crossover 1st order 𝑡 𝑚 PS , E SU(3) symmetric pointfit: 𝑎 + 𝑎 ⁄𝑁 fit: 𝑎 + 𝑎 ⁄𝑁 t + 𝑎 ⁄𝑁 solve: 𝑎 + 𝑎 ⁄𝑁 solve: 𝑎 + 𝑎 ⁄𝑁 t + 𝑎 ⁄𝑁 crossover 1st order 𝑡 𝑚 PS , E t SU(3) symmetric point crossover 1st order 𝑚 PS , E ⁄ 𝑇 E fit: 𝑎 + 𝑎 ⁄𝑁 fit: 𝑎 + 𝑎 ⁄𝑁 t + 𝑎 ⁄𝑁 solve: 𝑎 + 𝑎 ⁄𝑁 solve: 𝑎 + 𝑎 ⁄𝑁 t + 𝑎 ⁄𝑁 𝑡 𝑇 E FIG. 10. Continuum extrapolation of the critical endpoint. In upper panels, the left one is for √ t m PS , E vs 1 /N while theright one is for √ t m PS , E vs 1 /N t . In the lower panel m PS , E /T E (left) and √ t T E (right) are plotted as a function of 1 /N . ACKNOWLEDGEMENTS
This research was supported by Multidisciplinary Cooperative Research Program in CCS, University of Tsukubaand projects of the RIKEN Supercomputer System.
Appendix A: Wilson flow scale and pseudoscalar meson mass at zero temperature
Simulation parameters, results for the pseudoscalar meson mass am PS , and Wilson flow scale parameter √ t /a aresummarized in Table VI. Result of following combined fit is given in Table VII,( am PS ) = a (cid:18) κ − κ c (cid:19) + a (cid:18) κ − κ c (cid:19) , (A1) √ t a = b + b (cid:18) κ − κ c (cid:19) + b (cid:18) κ − κ c (cid:19) . (A2) [1] R. D. Pisarski and F. Wilczek, Phys. Rev. D , 338 (1984).[2] S. Gavin, A. Gocksch and R. D. Pisarski, Phys. Rev. D , R3079 (1994) [hep-ph/9311350].[3] A. Butti, A. Pelissetto and E. Vicari, JHEP , 029 (2003) [hep-ph/0307036]. 𝑁 𝑓 = 𝑁 t ≤ 10𝑁 t ≤ 12 physical point 𝑚 𝜂 𝑠 [ G e V ] 𝑚 [GeV ]Upper bound of critical endpoint FIG. 11. Columbia-like plot with axes m π and m η s in physical units. The blue symbol denotes the upper bound obtained inthis work N t ≤
12 while the red one is given in our previous study N t ≤
10 [34]. The physical point is also shown just for areference.[4] P. Calabrese and P. Parruccini, JHEP , 018 (2004) [hep-ph/0403140].[5] H. Shimizu and K. Yonekura, Phys. Rev. D , no. 10, 105011 (2018) [arXiv:1706.06104 [hep-th]].[6] K. Yonekura, JHEP , 062 (2019) [arXiv:1901.08188 [hep-th]].[7] C. Schmidt and S. Sharma, arXiv:1701.04707 [hep-lat].[8] H. T. Ding, PoS LATTICE , 022 (2017) [arXiv:1702.00151 [hep-lat]].[9] S. Sharma, PoS LATTICE , 009 (2019) [arXiv:1901.07190 [hep-lat]].[10] O. Philipsen, arXiv:1912.04827 [hep-lat].[11] F. R. Brown, F. P. Butler, H. Chen, N. H. Christ, Z. h. Dong, W. Schaffer, L. I. Unger and A. Vaccarino, Phys. Rev. Lett. , 2491 (1990).[12] F. R. Brown, N. H. Christ, Y. F. Deng, M. S. Gao and T. J. Woch, Phys. Rev. Lett. , 2058 (1988).[13] M. Fukugita, M. Okawa and A. Ukawa, Phys. Rev. Lett. , 1768 (1989).[14] H. Saito et al. [WHOT-QCD Collaboration], Phys. Rev. D , 054502 (2011) Erratum: [Phys. Rev. D , 079902 (2012)][arXiv:1106.0974 [hep-lat]].[15] C. Czaban and O. Philipsen, PoS LATTICE , 056 (2016) [arXiv:1609.05745 [hep-lat]].[16] Y. Iwasaki, K. Kanaya, S. Kaya, S. Sakai and T. Yoshie, Phys. Rev. D , 7010 (1996) [hep-lat/9605030].[17] S. Aoki et al. [JLQCD Collaboration], Nucl. Phys. Proc. Suppl. , 459 (1999) [hep-lat/9809102].[18] X. Liao, Nucl. Phys. Proc. Suppl. , 426 (2002) [hep-lat/0111013].[19] F. Karsch, E. Laermann and C. Schmidt, Phys. Lett. B , 41 (2001) [hep-lat/0107020].[20] F. Karsch, C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, E. Laermann and C. Schmidt, Nucl. Phys. Proc. Suppl. ,614 (2004) [hep-lat/0309116].[21] S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken and R. L. Sugar, Phys. Rev. D , 2531 (1987).[22] P. de Forcrand and O. Philipsen, JHEP , 077 (2007) [hep-lat/0607017].[23] M. A. Clark, A. D. Kennedy and Z. Sroczynski, Nucl. Phys. Proc. Suppl. , 835 (2005) [hep-lat/0409133].[24] M. A. Clark and A. D. Kennedy, Phys. Rev. Lett. , 051601 (2007) [hep-lat/0608015].[25] D. Smith and C. Schmidt, PoS LATTICE (2011) 216 [arXiv:1109.6729 [hep-lat]].[26] P. de Forcrand, S. Kim and O. Philipsen, PoS LAT , 178 (2007) [arXiv:0711.0262 [hep-lat]].[27] G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, PoS LAT , 182 (2007) [arXiv:0710.0998 [hep-lat]].[28] H.-T. Ding, A. Bazavov, P. Hegde, F. Karsch, S. Mukherjee and P. Petreczky, PoS LATTICE , 191 (2011)[arXiv:1111.0185 [hep-lat]].[29] L. Varnhorst, PoS LATTICE , 193 (2015). TABLE VI. Simulation parameters κ , N s , and N t and measured √ t /a and am PS at β = 1 .
77, 1 .
78, 1 .
80, and 1 .
81. Note thatthe data at β = 1 .
78 is updated compared with the previous work [34]. β κ N s N t √ t /a am PS κ c and coefficients for the pseudoscalarmeson mass am PS and the Wilson flow parameter √ t /a for β = 1 .
77, 1 .
78, 1 .
80, and 1 .
81. Note that the data at β = 1 .
78 isupdated compared with the previous work [34]. β κ c a a b b b χ / d . o . f . fit range κ > − − − − − − , no. 7, 074505 (2017) [arXiv:1701.03548 [hep-lat]].[31] M. Creutz, PoS LATTICE , 007 (2007) [arXiv:0708.1295 [hep-lat]].[32] C. Bernard, M. Golterman, Y. Shamir and S. R. Sharpe, Phys. Rev. D , 114504 (2008) [arXiv:0711.0696 [hep-lat]].[33] X. Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda and A. Ukawa, Phys. Rev. D , no. 1, 014508 (2015) [arXiv:1411.7461[hep-lat]].[34] X. Y. Jin, Y. Kuramashi, Y. Nakamura, S. Takeda and A. Ukawa, Phys. Rev. D , no. 3, 034523 (2017) [arXiv:1706.01178[hep-lat]].[35] Y. Iwasaki, Report No. UTHEP-118 (1983), arXiv:1111.7054 [hep-lat].[36] S. Aoki et al. [CP-PACS and JLQCD Collaborations], Phys. Rev. D , 034501 (2006) [hep-lat/0508031].[37] Y. Nakamura and H. Stuben, PoS LATTICE , 040 (2010) [arXiv:1011.0199 [hep-lat]].[38] Y. Kuramashi, Y. Nakamura, S. Takeda and A. Ukawa, Phys. Rev. D , no. 11, 114507 (2016) [arXiv:1605.04659 [hep-lat]]. [39] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. , 2635 (1988).[40] M. L¨uscher, JHEP , 071 (2010) Erratum: [JHEP , 092 (2014)] [arXiv:1006.4518 [hep-lat]].[41] S. Borsanyi et al. , JHEP , 010 (2012) [arXiv:1203.4469 [hep-lat]].[42] P. de Forcrand and M. D’Elia, PoS LATTICE , 081 (2017) [arXiv:1702.00330 [hep-lat]].[43] H. Ohno, Y. Kuramashi, Y. Nakamura and S. Takeda, PoS LATTICE2018