The Chiral Phase Transition in three-flavor QCD from Lattice QCD

We analyze the pseudo-critical behavior of three-flavor QCD using highly improved staggered quarks (HISQ) on lattices with temporal extent $N_\tau =8$ and for quark masses corresponding to a pseudoscalar Goldstone mass in the range $80 ~ {\rm MeV} ~ \lesssim ~ m_\pi ~ \lesssim ~ 140 ~ {\rm MeV}$. Our findings are consistent with the occurrence of a second order chiral phase transition at vanishing values of the quark masses. The chiral phase transition temperature at this finite value of the lattice spacing is determined to be $T_c = 98_{-6}^{+3}~{\rm MeV}$. A comparison with a corresponding analysis performed in (2+1)-flavor QCD suggests that the continuum limit extrapolated chiral phase transition temperature in 3-flavor QCD will turn out to be below $90 ~ {\rm MeV}$.

We analyze the pseudo-critical behavior of three-flavor QCD using highly improved staggered quarks (HISQ) on lattices with temporal extent Nτ = 8 and for quark masses corresponding to a pseudoscalar Goldstone mass in the range 80 MeV < ∼ mπ < ∼ 140 MeV. Our findings are consistent with the occurrence of a second order chiral phase transition at vanishing values of the quark masses. The chiral phase transition temperature at this finite value of the lattice spacing is determined to be Tc = 98 +3 −6 MeV. A comparison with a corresponding analysis performed in (2+1)-flavor QCD suggests that the continuum limit extrapolated chiral phase transition temperature in three-flavor QCD will turn out to be below 90 MeV.

I. INTRODUCTION
The nature of the QCD chiral phase transition as a function of the number of light flavors and the value of the quark masses, has been the subject of intense and ongoing study ever since the first work of Pisarski and Wilczek [1]. It is now well-established that the transition is a crossover for the physical values of the light and strange quark masses [2]. The transition is, however, expected to become second order and first order in the two flavor and three flavor chiral limits, respectively. Consequently, a critical value for the strange quark mass is expected to exist where a tri-critical point separates the region of first and second order chiral phase transitions. This is sketched in the well-known Columbia plot for the quark mass dependence of phase transitions in QCD [3].
Lattice studies have provided strong evidence for the expected second order nature of the chiral phase transition in the 2-flavor chiral limit also for physical value of the strange quark mass. The transition is expected to belong to the 3-d, O(4) universality class if the anomalous U (1) A symmetry remains sufficiently broken; otherwise, a second order phase transition belonging to the U (2) × U (2) universality class is also possible [4]. However, as critical exponents in both universality classes are similar a differentiation between both possibilities will be difficult in practice [4].
While the axial anomaly plays a decisive role for the symmetry breaking pattern in 2-flavor QCD, it is of less relevance in three-flavor QCD [4]. Irrespective of whether or not the U (1) A symmetry is effectively restored at the chiral phase transition temperature, a renormalization group (RG) analysis [5] suggests that the transition is first order, as originally predicted in [1]. The functional RG (FRG) analysis of the three-flavor chiral transition, however, suggests that a first order transition, only occurs for rather small values of the pion mass, m π < ∼ 25 MeV. This is in accordance with lattice QCD calculations using staggered fermions. Although such calculations find first order transitions in calculations with unimproved gauge and fermion actions [6], the bounds on the critical mass show a strong cut-off and discretization scheme dependence [7][8][9]. In calculations using improved staggered fermions no direct evidence has been found for a first order transition on lattices with temporal extent N τ = 6 for m π > ∼ 80 MeV and a bound on the pseudoscalar Goldstone mass, above which no first order transition exists, has been estimated to be m c π 50 MeV [10]. In calculations with O(a) improved Wilson fermions [11] first order transitions have been found at non-zero values of the quark masses and the bound on the critical mass is weaker, m c π < ∼ 110 MeV. This bound, however, also is consistent with a continuous transition in the continuum limit [11,12]. In fact, a recent analysis of the order of the chiral transition as function of the number of flavors, performed with staggered fermions and extrapolated to the continuum limit [12], suggests that the chiral phase transition in three-flavor QCD is second order, which is in contrast to RG analyses. However, as has been pointed out, such analyses are based on a Landau-Ginsburg effective action for the order parameter, which is arrived at by integrating out all gauge degrees of freedom ending up with a φ 4 effective Lagrangian for the order parameter field. The role of gauge fluctuations, however, is subtle and may also influence the order of the chiral phase transition [13,14]. It also has been argued that a φ 6 contribution to the effective Lagrangian for the order parameter may be of relevance and may allow for a second order chiral phase transition to occur in QCD with number of massless flavors being larger than two [12]. A continuous transition in the chiral limit of three-flavor QCD thus may not be ruled out entirely.
In addition to the exploration of the flavor dependence of the chiral phase transition the determination of the chiral phase transition temperature is of central interest. The chiral phase transition temperature is expected to drop with increasing number of flavors, n f , [15] and will eventually vanish at a critical value of the number of flavors, corresponding to the conformal limit of QCD [16], n * f ∼ 10. In recent years, different lattice estimates for both the chiral crossover T pc as well the 2+1flavor chiral transition temperature T n f =(2+1) c , obtained using different actions and with different choices of observables for setting the scale, have converged [17][18][19][20].
The crossover temperature, occurring for physical values of the two (degenerate) light quark masses (m u = m d ) and a strange quark mass, m phys s 27m u , has been determined to within 1% to be T pc = 156.5 (1.5) MeV [17] and T pc = 158.0 (0.6) MeV [18]. In the limit of vanishing two light quark masses, keeping the strange quark mass fixed to its physical value, the phase transition temperature in the continuum limit has been found to be T [20]. A value consistent with this number has recently been obtained also in calculations with twisted mass Wilson fermions in a (2 + 1 + 1)flavor calculation [21].
In three-flavor QCD calculations the boundary for a first order transition has been found at T n f =3 c = 134(3) MeV and a pseudoscalar mass, m c π 110 MeV [11]. Irrespective of whether below this mass indeed a first order transition exists or whether this only corresponds to pseudo-critical temperature, this value for T n f =3 c will drop further when approaching the chiral limit, leading to a chiral transition temperature below that of (2+1)-flavor QCD. In general it is expected that the phase transition temperature drops with increasing n f . In fact, the FRG analysis presented in [15] suggests that the chiral phase transition temperature in QCD with n f = 3 massless quarks is smaller by about 25 MeV compared to the 2-flavor case.
In the study presented here we want to further explore the nature of the chiral phase transition in 3-flavor QCD and, in particular, provide a first estimate of the chiral phase transition temperature based on calculations with the highly improved staggered quark (HISQ) action. These calculations extend earlier studies performed with the HISQ action [10] on coarser lattices. We work here at N τ = 8 compared to N τ = 6 used in [10].
The paper is organized as follows. In the next section we introduce the chiral observables we will study to determine the chiral phase transition temperature in three-flavor QCD. Section III summarizes basic relations needed for our discussion of finite size scaling (FSS) of chiral observables. In Sec. IV we summarize some results on the chiral transition in three-flavor QCD obtained previously on coarser lattices and present a first comparison with data on the disconnected part of the chiral susceptibility obtained in our new study. In Sec. V we finally present our results for the FSS analysis of the chiral order parameter and its susceptibility, from which we deduce the chiral phase transition temperature on lattices with temporal extent N τ = 8. We give our conclusions in Sec. VI.

A. Simulation parameters
The results that we present here were obtained from lattice QCD simulations with three degenerate flavors.
In our calculations with staggered fermions we use the HISQ action and a tree-level improved Symanzik gauge action. This framework is identical to that used previously in finite temperature calculations for (2+1)-flavor QCD [22][23][24] as well as the three-flavor QCD calculations [10] performed on lattices with temporal extent N τ = 6.
The temporal extent of our lattices was fixed at N τ = 8 throughout while the spatial extent was chosen to be one of N σ = 24, 32 or 40. In order to control finitevolume effects in our calculations we typically performed calculations at a given quark mass value for two different values of N σ . The larger one is chosen such that m π L ≥ 3 in the region of the pseudo-critical temperatures.
We used the Bielefeld GPU code [25] to generate around 10,000-50,000 hybrid Monte Carlo trajectories separated by 0.5 Time Units (TU) for N σ = 40 and 1 TU for the smaller volumes. These data sets have been generated in 10 independent streams that have been decorrelated initially using about 200 trajectories. The rational hybrid Monte Carlo algorithm [26,27] was used to generate the configurations and the step sizes were tuned so as to achieve acceptance rates of 60-80%. We saved gauge field configurations after every 5th time unit and performed calculations of various chiral observables on these configurations.
We generated data sets for different values of the quark masses at up to 17 values of the temperature in the range 110 MeV < ∼ T < ∼ 170 MeV. As a guidance for our choice of bare quark masses we used the line of constant physics (LCP) determined in Ref. [24] for the case of (2+1)flavor QCD. This LCP defines the value of the strange quark mass, m phys s (β), as a function of the gauge coupling β. It is tuned to its physical value by demanding the mass of the η ss meson to stay constant on this LCP. We choose different sets of three degenerate light quark masses, m q , corresponding to m q = H m phys s (β), with H = 1/27, 1/40, 1/60 and 1/80. Further details of our scale setting are discussed in the next subsection.

B. Scale setting
In the continuum limit the relation between lattice cutoff and gauge coupling β = 10/g 2 is controlled by the universal, asymptotic 2-loop β-function of three-flavor QCD, with b 0 = 9/(16π 2 ) and b 1 = 1/(4π 4 ) and Λ L denoting the QCD Λ-parameter for the three-flavor lattice discretization scheme. Eq. (1) receives corrections at nonzero values of the lattice spacing that depend on the observable used to set the scale. We use here the nonperturbative β-functions f K a(β) and a/r 1 (β) determined in (2+1)-flavor QCD for the kaon decay constant, f K a, as well as the parameter r 1 /a deduced from the slope of the heavy quark potential [24].
In all our figures we use the f K scaling function to introduce a temperature on lattice with the temporal extent N τ = 8, where we set the scale for the temperature by using the value f K = 156.1/ √ 2 MeV for the kaon decay constant as has been done also in the determination of the chiral phase transition temperature in (2+1)-flavor QCD [20]. We note that this value agrees within errors with the Flavor Lattice Averaging Group (FLAG) average [28]. This scale setting also is used in all our fits to data. In some cases we use fits based on the r 1 /a parametrization to obtain some idea about systematic errors in our results that are not yet continuum extrapolated. As a physical value for r 1 we use r 1 = 0.3106 fm.
As mentioned in the previous subsection we choose the three degenerate quark masses in our simulations as fractions of the strange quark mass used in the (2+1)-flavor QCD calculations to define a LCP. For this we use the parametrization of m phys s (β) given in Appendix C of Ref. [24]. Our choice of three-flavor quark mass, m phys s /80 < ∼ m q < ∼ m phys s /27 then corresponds to a light pseudoscalar Goldstone mass in the range 80 MeV < ∼ m π < ∼ 140 MeV in the continuum limit of (2+1)-flavor QCD.

C. Chiral observables
On each saved gauge field configuration, we calculated the chiral condensate, ψ ψ , and its susceptibility, χ tot , which are obtained from the free energy density of threeflavor QCD, f (T, m q ) = −(T /V ) ln Z(T, V, m q ), as first and second derivative with respect to the quark mass m q , The chiral susceptibility is represented in terms of disconnected, χ disc , and connected, χ con contributions, χ tot = χ disc + χ con . These chiral observables are given in terms of the inverse of the staggered fermion matrix, D q (m q ), and its higher powers, We evaluate these chiral observables using up to 100 Gaussian random vectors to evaluate the trace of the inverse of D q . After calculating the traces on each configuration, the observables were calculated by dividing the total number of configurations into 10 bins and using the jackknife procedure.
In the continuum limit the chiral observables defined in Eq. 5 require additive and multiplicative renormalization. The former arises from a ultra-violet divergent contribution to the chiral condensate, ψ ψ = ψ ψ ren + (c U V /a 2 )m q . We use the strange quark mass m phys s (β) for multiplicative renormalization of the chiral observables to define the order parameter, M , and its susceptibility, χ M as, Here m phys s is our reference mass in three-flavor QCD, whose choice is motivated by the LCP determined in (2+1)-flavor QCD calculations, and H = m q /m phys s . This makes the observables dimensionless and RGinvariant up to logarithmic corrections. As we do not have direct access to the UV-divergent contribution, we instead evaluate the unsubtracted observables With this we have Here the 1/a 2 divergence manifests itself as N 2 τ , since the continuum limit is obtained by sending N τ → ∞ at fixed T . For fixed, N τ , the UV-term picks up a quadratic temperature dependence and a logarithmic correction arising from the anomalous dimension of the strange quark mass term m phys s . In fact, in a small temperature interval around the chiral phase transition temperature the temperature dependence of the UV-term is to a good approximation linear in T . This is shown in Fig. 1.
Using the parametrizations of the strange quark mass in lattice units, m phys s (β), and the kaon decay constant, f K a(β), given in [24] we evaluated the normalization factor multiplying the UV constant c U V and interpolated the result with a quadratic ansatz. For this we find with δT = (T − 100)/100. This quadratic approximation is shown as a solid line in Fig. 1. Here we also give the result of a linear fit in the temperature interval [100 : 125 MeV]. We note that this UV-term can be treated as part of the regular contributions in a (finite-size) scaling analysis. We will discuss this in more detail in the next section. In order to eliminate the divergent UV contributions explicitly in chiral observables we introduce the difference In fact, this observable can also be considered as an order parameter for the chiral phase transition. At high temperatures it vanishes, M χ ∼ H 3 , and at low temperatures it equals the order parameter M at H = 0 but receives different corrections at O( √ H).

III. SCALING AND FINITE-SIZE SCALING OF CHIRAL OBSERVABLES
As there is increasing evidence that also three-flavor QCD will have a second order phase transition in the chiral limit, or at a rather small quark mass, it is appropriate to analyze thermodynamic observables in three-flavor QCD also in terms of relevant scaling functions. As we are currently analyzing the chiral limit at fixed values of the cut-off the universality class of 3-dimensional, O(2) symmetric models would be appropriate, while the Z(2) universality class is of relevance, if a second order phase transition occurs at a small value of the quark mass, m c q . In the vicinity of a critical point, (T c , m c q ), the thermodynamic free energy, f (T, m q ), can be resolved into singular and regular contributions, f (T, m q ) = f s (T, m q ) + f r (T, m q ). The temperature and quark mass dependence of f s (T, m q ) is expressed in terms of a universal scaling function, f f (z), which is characteristic for a particular universality class. Thus we have with t and H being dimensionless variables constructed from the temperature T and quark mass m q , The constants β and δ are critical exponents of the 3D, O(2) universality class for which we use [29], while t 0 and h 0 are non-universal constants that are introduced to fix the overall normalization of the order parameter M [19]. For small quark masses, in the vicinity of the chiral transition temperature, T c , the dominant contributions to the order parameter M and its susceptibility χ M arise from the singular part of the free energy. Scaling relations for these observables are then obtained by taking derivatives of f s (z) with respect to H. We have and respectively. Here f G (z) and f χ (z) are also universal functions of the scaling variable z, that can be obtained from f f (z). Both these functions have been determined numerically using high-statistics Monte Carlo simulations for the 3D, O(2) and O(4) universality classes [29][30][31][32]. We will make use of the implicit parameterization provided for these functions for the O(2) case in Ref. [33].
In the scaling regime the difference of order parameter and chiral susceptibility, introduced in Eq. 11, is given in terms of these scaling functions From Eqs. (15) and (16), it is readily seen that irrespective of the quark mass m q . Curves for different quark masses will have a unique crossing point at T = T c . We also will analyze the related observable, obtained from the ratio of the non-subtracted chiral order parameter and its susceptibility, M b /χ M b for various H at their corresponding pseudo-critical temperature [34]. The pseudo-critical temperatures defined by the maximum of the chiral susceptibility correspond to a specific value of the scaling variable z, i.e. z = z p , with z p = 1.58(4) and f G (z p )/f χ (z p ) = 1.58 and 1.51 in the O(2) and Z(2) universality classes, respectively. We note that the RHS of Eq. 19 is solely determined by the universal contribution which is an important aspect for the analysis of the nature of the chiral transition [34]. Note that Eqs. 18 and 19 only hold provided the quark mass is sufficiently close to m c q so that the regular contributions coming from f r (T, m q ) can be ignored.
In addition to regular contributions, Eq. (18) also receives corrections when the system size L is finite. Then the scaling functions f G (z) and f χ (z) in Eqs. (15) and (16) must be replaced by the corresponding Finite-Size Scaling (FSS) functions f G,L (z, z L ) and f χ,L (z, z L ). Here z L = L 0 /Lh ν/βδ is a second scaling variable and ν = β(δ + 1)/d with d = 3 is the critical exponent controlling the divergence of the correlation length at the critical point. L 0 is another non-universal constant that can be fixed through a normalization condition, e.g. for the chiral condensate [35]. In the thermodynamic limit (L → ∞) the finite size scaling functions, f G,L (z, z L ) and f χ,L (z, z L ), go over to their infinite-volume equivalents f G (z) ≡ f G,L (z, 0) and f χ (z) ≡ f χ,L (z, 0).
The infinite and finite volume scaling functions have been determined for the 3D, O(2) and O(4) cases [33,35]. In our FSS analysis we use a rational polynomial parametrization of the FSS functions [36] similar to what has been used also for the analysis of FSS in O(4) spin models [35] As we will see, the observed transition is a crossover for all the quark masses that we studied. According to the standard picture of the phase diagram [3], this transition should turn into a first order phase transition if the quark mass is less than a certain value m q < m c q . For m q = m c q then, the transition will be second order belonging to the 3D, Z(2) universality class. If on the other hand, the expected first order region is absent, then the transition will be second order belonging to the 3D, O(2) universality class in the three-flavor chiral limit for fixed N τ . Since we did not find any evidence of a non-zero critical quark mass m c q in our study, we will assume m c q = 0 in the following and use the O(2) scaling functions for the rest of this work.

IV. CHIRAL ORDER PARAMETER AND ITS SUSCEPTIBILITY
In the analysis of the chiral transition in (2+1)-flavor QCD it has been found that in the continuum limit the pseudo-critical temperature is about T pc 156.5 MeV at physical values of the quark masses and drops in the limit of vanishing light quark masses to the value of the chiral phase transition temperature of about 132 MeV. On lattices with temporal extent N τ = 8 the corresponding pseudo-critical and chiral phase transition temperatures are larger by 5 MeV and 12 MeV, respectively. As we expect the pseudo-critical temperatures in three-flavor QCD to be smaller than those of the (2+1)-flavor theory we explored in our calculations a range of temperatures 100 MeV < ∼ T < ∼ 170 MeV.

A. Chiral observables on Nτ = 8 lattices
We present our results for the non-subtracted chiral order parameter, M b , and its susceptibility, χ M b , calculated  on lattices with temporal extent N τ = 8, in Fig. 2. The order parameter M b varies rapidly but smoothly, starting from a high value and decreasing to a low value over the temperature range considered here. For three-flavor quark masses corresponding to H = 1/27 we observe the most rapid change of M b at temperatures around T (135 − 140) MeV. In this temperature range also the chiral susceptibility, χ M b , has its maximum. This indicates that at this quark mass value the pseudo-critical temperature in three-flavor QCD is shifted by about (20)(21)(22)(23)(24)(25) MeV relative to that of (2+1)-flavor QCD. This trend also persists for the smaller quark masses examined by us.
Except for the case H = 1/60 we show in Fig. 2 results for two different volumes. While the volume dependence is negligible in M b and χ M b slightly above the corresponding pseudo-critical temperatures, evidence for a characteristic volume dependence is seen at smaller temperatures. In fact, contrary to what would be expected at or close to a second or first order phase transition, we see no increase in the peak height of χ M b with increasing volume. Instead we observe a slight decrease of the  peak height of χ M b and a shift of the peak position towards larger temperatures as the volume is increased. At the same time the peak becomes more pronounced with increasing volume. This behavior is reminiscent of the finite volume effects known from an analysis of finite size scaling functions in 3-d, O(N ) symmetric spin models [33,35]. At fixed quark mass we also see evidence in the data for a slight volume dependence of the non-subtracted chiral order parameter M b . This is clearly visible in Fig. 2 (top). In Fig. 3 we show the chiral order parameter as function of H for several values of the temperature. As can be seen for T ≥ 140 MeV the order parameter depends linearly on the quark mass and extrapolates smoothly to zero for H → 0. At lower temperatures the quark mass dependence of the location of the maximum in χ M b suggests that in the chiral limit all our calculations correspond to a range of temperatures in the chirally symmetric phase.
For a first analysis of the quark mass dependence of M b we thus fitted the data on the largest lattices available to an ansatz, Results for the exponent e are shown in the inset of this figure. Within errors it is consistent with unity for T ≥ 140 MeV and decreases continuously with decreasing temperature. At the lowest temperature, T 110 MeV, we find e 0.27, which is larger but compatible with the exponent one expects to find at a critical point belonging to the 3-dimensional O(2) or Z(2) universality classes, i.e. e ≡ 1/δ 0.21.
It is evident from Fig. 2 (top) that finite volume effects are relevant at lower temperatures and need to be treated carefully to arrive at results in the thermodynamic limit. We will analyze this volume dependence further in Sec. III.   Table I. Results for the pseudo-critical temperature obtained from the location of peaks in χ disc M b for Nτ = 6. They are obtained from Ref. [10] as discussed in the text.

B. Comparison with results from Nτ = 6 lattices
Before going into a more detailed analysis of the volume dependence and universal scaling behavior of our results, obtained on lattices with temporal extent N τ = 8, we want to compare with earlier results obtained on coarser lattices with temporal extent N τ = 6. In that case results exist only for the disconnected chiral susceptibility, χ disc M b , and critical couplings had been extracted from the peak position of χ disc M b . For better comparison with our current results we thus also show this susceptibility obtained by us on lattices with temporal extent N τ = 8, in Fig. 4. In Tab. I and II we give pseudo-critical temperatures obtained from the peak positions of χ disc M b for N τ = 6 and 8.
In both cases we used the temperature scale determined from the parametrization 1 of f K a.  Table II. Results for pseudo-critical temperatures obtained from the peak of the disconnected part of the chiral susceptibility on the largest lattice used in simulations on lattices with temporal extent Nτ = 8.
The calculations for N τ = 6 [10] have been performed at various temperature values, keeping the quark mass fixed in units of the lattice spacing while in our current calculations on N τ = 8 lattices the quark mass is varied along a line of constant physics by keeping the ratio H fixed. Using the parametrization m phys s (β) we obtained the corresponding ratio H at the pseudo-critical coupling β c (m q ) also for the N τ = 6 data sets given in [10]. The resulting ratios H −1 are given in Tab. I.
The pseudo-critical temperatures, obtained for N τ = 6 and 8 from maxima in the disconnected susceptibilities 2 are shown in Fig. 5. Here the data are plotted versus H 1/βδ , with βδ = 1.67 for the 3D, O(2) universality class. We fitted these data using an ansatz inspired by the universal scaling of pseudo-critical temperatures in the vicinity of second order phase transitions [22], continuum limit both ways of scale setting will lead to a unique temperature scale [37]. 2 Further details on the determination of these maxima are discussed in Section V.B The term proportional to b c arises from the leading regular contribution to the free energy, being proportional to H. However, as can be seen in Fig. 5 our data are not sensitive to this correction. In fact, they are well described by a straight line fit in terms of H 1/βδ . The fits shown in the figure for b c = 0 yield T disc pc = 108.3(4) MeV for N τ = 6 and 95.8 (7) MeV for N τ = 8, respectively. They both have a χ 2 /dof less than unity.
This result, as well as the quark mass dependence of the order parameter shown in the previous subsections motivated a more detailed scaling analysis of the threeflavor results for the chiral order parameter and its susceptibility, which we present in the following section.

V. FINITE SIZE SCALING ANALYSIS OF CHIRAL OBSERVABLES
We want to improve here on our determination of the chiral transition temperature obtained in the previous sections through an analysis of disconnected part of the chiral susceptibility on the largest lattices available. We will analyze the finite size dependence of the unsubtracted chiral order parameter M b and its susceptibility As can be seen in Fig. 2 M b (T, m q , N σ ) and χ M b (T, m q , N σ ) show a sizeable volume dependence for all the quark masses and for temperatures T < T pc (m q ). Above the pseudo-critical temperature this volume dependence is significantly weaker. We define the pseudocritical temperature, T pc (m q , N σ ), on lattices with spatial extent N σ as the location of the maximum of χ M b for a given quark mass and volume. We also define T pc (m q ) as the infinite volume pseudo-critical temperature for a given quark mass.

A. FSS analysis of the chiral condensate and its susceptibility
We start with an analysis of the difference of the chiral order parameter and its susceptibility, introduced in Eq. 11. As discussed there this difference eliminates any dependence of the chiral observables linear in H and thus also is independent of leading correction to scaling relations arising from regular contributions to the chiral observables. In Fig. 6, we thus compare M χ (T, m q , N σ ) to the difference of FSS functions, (22) Aside from the critical temperature in the chiral limit, T c , a fit to this ansatz involves the three non-universal parameters, h 0 , z 0 = h a finite volume, given in the previous sections, that the chiral phase transition will be located at a temperature below 110 MeV. We thus perform fits only for temperatures close of the chiral transition temperature, i.e. for temperatures below 121 MeV. The results of such fits in different fit ranges for the quark masses, H ∈ [0, H max ] are given in Tab. III and shown in Fig. 6 . As can be seen the fits yield a good χ 2 /dof when leaving out the data sets for our largest quark mass ratio H = 1/27. However, even when including these data sets we obtain a reasonable fit result for 28 data points using only the three non-universal parameters of the 3D, O(2) scaling functions.  Table III. Non-universal parameters of the O(2) scaling functions obtained from fits of the difference of order parameter and chiral susceptibility, Mχ.
A related observable is the ratio of the chiral susceptibility and the chiral order parameter which has been used in [20] to determine the chiral phase transition temperature in (2+1)-flavor QCD, (23) Unlike the difference M χ analyzed above this ratio is sensitive to regular contributions as well as the subtraction of a UV divergent term. It, however, eliminates the explicit dependence on the non-universal scale parameter h 0 . For the contributions of the regular term we use a leading order Taylor expansion in the vicinity of the chiral transition temperature, T c , Solid: H max =1/27 Dotted: H max =1/40 In the temperature range analyzed by us this regular term, also takes care of the UV divergent contribution to the order parameter, as is apparent from Fig. 1. We note that in this fit ansatz the non-universal parameter h 0 does not explicitly appear as an independent fit parameter, but is absorbed as an overall factor h 1/δ 0 in the fit parameters a i of the regular term.
We perform a FSS analysis of the ratio given in Eq. 23. We again performed fits in different fit intervals H ∈  Table IV. Parameters of fits to the ratio HχM b /M b using the ansatz given in Eq. 23. The upper half shows fit results obtained with a regular term including terms linear in temperature, while in the lower half only results for fits with a constant term in the regular part are shown. For the smallest bound on the quark mass ratio, Hmax = 1/80, we also show the result from a fit the uses only the singular part of the fit ansatz.
We note that the fits based on differences of M b and χ M b and the fit of the ratio Hχ M b /M b are in excellent agreement with each other, although only the latter is sensitive to regular contributions to the chiral observables. The parameters of the singular part of these fits are consistent with each other within errors. The parameters of the regular term seem to be quite sensitive to the upper bound (H max ) for the set of quark masses used in the fit. They are not well determined in the small temperature interval used for these fits. In fact, the χ 2 /dof changes only little when only the leading temperature dependent term in the fit is used. The resulting parameters for this fit are shown in the lower half of Tab. IV. For the case of H max = 1/80 the regular contribution vanishes within errors. We thus also give the result of a fit that only uses the singular term in the fit ansatz. As can be seen this suffices to obtain a good fit at this small value of the quark mass and yields fit results for the non-universal parameters of the scaling functions that are in good agreement with those obtained in a larger quark mass range by including contributions from a regular term.
We also performed fits that include a quadratic correction in the regular term. This allows us to enlarge the fit range for temperatures above the transition temperature up to values that correspond to the pseudo-critical temperature for the largest quark mass contributing to the fit. This way we still obtain fits with χ 2 /dof close to unity. Results of these fits are summarized in Tab. V. We note that these fits give large coefficients with opposite sign for terms linear and quadratic in temperature. The reduced temperature factor tt 0 = (T − T c )/T c is about 1/3 at the upper end of the fit interval, T 130 MeV, which is comparable with the difference in magnitude of the coefficients a 1 and a 2 , i.e. contributions from linear and quadratic terms in the regular part compensate each other to a large extent. We thus conclude that the parameters of the singular part of our fits are well determined, while the parametrization of the regular is not strongly constrained.

B. FSS analysis of the pseudo-critical temperature
The chiral phase transition temperature T c can be determined from a finite size scaling analysis of the pseudocritical temperature T pc (m q , N σ ), which can be obtained straightforwardly from our data for χ M b (T, m q , N σ ).
To determine T pc (m q , N σ ) from the maxima of χ M b (T, m q , N σ ), 100 bootstrap samples were constructed at the level of saved gauge configurations for each volume and quark mass, followed by fitting the χ M data in the peak region to a quadratic ansatz for each such boot-  Table VI. Results for pseudo-critical temperatures obtained from the maxima of the total chiral susceptibility on lattice with temporal extent Nτ = 8 and different spatial lattice sizes.
strap sample. Starting with a minimum of three points, the fit interval was increased by adding one point at a time on either side of the peak to go up to a maximum of five points in all. In this way, we performed around 3-5 fits for each bootstrap data set. The final value of T pc (m q , N σ ) for a given quark mass and volume was obtained by taking the mean of all the fits performed on the 100 bootstrap samples and error was the standard deviation of this distribution. The resulting estimates for the finite-volume pseudo-critical temperatures for all available lattice sizes are given in Tab. VI. Similar procedure was used to locate the peak position of χ disc M b for the highest available volumes for each quark mass.
We present our results for T pc (m q , N σ ) for all the quark masses and volumes in Fig. 8. We fitted these data using the FSS scaling ansatz for T pc deduced from the FSS ansatz for χ M b (T, m q , N σ ) = ∂M b /∂H using Eq. 23, where z p (0) = z p gives the location of the maximum of the infinite volume O(2) scaling function f χ (z) and z p (z L ) is a parametrization of the finite volume dependence of this peak location [36], This parametrization holds for z L ≤ 0.9. The resulting fits for the pseudo-critical temperatures are also shown in Fig. 8, again for the cases with and without including the largest quark mass, H=1/27, used in our calculations. The chiral phase transition temperatures obtained from these fits are given in Tab Figure 8. Results for the peak location Tpc(mq, Nσ) plotted versus the finite volume scaling variable (zL/zL,0) 3 = (Nτ /(NσH ν/βδ )) 3 . The data are plotted as points whereas the bands are the fit results for two different choices of the fit interval in H.
In Fig. 9 we show extrapolation of the pseudo-critical temperatures to the chiral limit using T pc obtained on the largest volumes available for each mass. This gives T c = 103(3) MeV and 97(4) MeV with H max = 1/27 and 1/40, respectively. These numbers are in complete agreement with the previously found numbers from the FSS fits. In fact, the FSS fits shown in Fig. 8 suggest that the pseudocritical temperatures, obtained on our largest lattices, differ from the infinite volume extrapolated results by less than 2 MeV.
We took into account the systematic differences in our fits resulting from changes of the fit range for H as well as the fit ansätze that include or leave out contributions from regular terms in the different observables we fitted. Averaging over all these fit results for T c we obtain for the chiral phase transition temperature in three-flavor QCD on lattices with temporal extent N τ = 8, This result is consistent within errors with all fit results for T pc presented in Tabs. III to VII and shown in Fig. 10.
Weighting with the Akaike information criterion or by the inverse of the squared error form each fit yields estimates for T c which agree very well with the value quoted in Eq.

27.
Our result for the chiral phase transition temperature, obtained from a scaling analysis for non-zero quark masses corresponding to pseudoscalar Goldstone masses m π > ∼ 80 MeV, relies on the occurrence of second order phase transition in the chiral limit. In support of this  we show in Fig. 11 the ratio M b /χ M b , evaluated at the position of the maxima of χ M b for fixed H. The dashed curves, following the data, show the result of the fit to Hχ M b /M b given in Eq. 23, using a regular term up to quadratic order in T and evaluated at T pc obtained from Eq. 25, with z p (z L ) approximated by its infinite volume value z p (0). This ratio is compared to straight lines with slopes given by the parameter ratio of scaling functions for the 3D, O(2) universality class (red) and three lines corresponding to the Z(2) universality class, respectively. The latter is shown for three hypothetical critical masses H c , below which a region of first order transitions might exist in 3-flavor QCD. This comparison puts stringent bounds on the possible on a critical quark mass below which a first order chiral phase transition may still occur. For quark mass ratios smaller than H 1/60 results for Hχ M b /M b , evaluated at T pc (m q , N σ ) are in good agreement with the parameter free universal O(2) scaling function. This leaves little room for a non-vanishing H c and a possible approach of the data for M b /χ M b to the universal Z(2) scaling function. Our results thus suggest a continuous chiral phase transition in the 3D, O(2) universality class at least at finite values of the cut-off corresponding to N τ = 8.

VI. CONCLUSIONS
We have presented results on the chiral phase transition in three-flavor QCD. Our calculations have been performed for finite values of the lattice spacing, corresponding to N τ = 8. For the range of quark masses, corresponding in the continuum limit to light pseudoscalar Goldstone masses in the range 80 MeV ≤ m π ≤ 140 MeV we find no direct evidence for a conjectured first order phase transition. In the transition region we observe pseudo-critical behavior with a finite volume dependence that is consistent with the expected FSS behavior in the 3D, O(2) universality class. The parameter free comparison of the ratio M b /χ M b with 3D, O(2) and Z(2) scaling functions gives further support for a second order phase transition in the chiral limit of three-flavor QCD, at least on lattices with temporal extent N τ = 8. For the chiral phase transition temperature at these non-vanishing values of the lattice spacing we find T c = 98 +3 −6 MeV. In (2+1)-flavor QCD the chiral phase transition at this value of the cut-off was about 10 MeV larger than the continuum limit extrapolated phase transition temperature. Assuming that cut-off effects are of similar magnitude also in three-flavor QCD, in the continuum limit the chiral phase transition temperature in 3-flavor QCD is expected to be below 90 MeV. This needs to be confirmed in future calculation.
We also find that the pseudo-critical temperature at a pseudoscalar mass of about 110 MeV, which corresponds to H 1/40, is about (30-40) MeV larger than T c , i.e. it is about 135 MeV as shown in Fig. 8. This is indeed consistent with the estimate of the transition temperature at this value of the pseudoscalar mass obtained with Wilson fermions [11]. In this case, however, the transition tem- perature is identified as the location of the end point of a region of first order phase transitions. The difference between the two results obtained within the staggered and Wilson fermion discretization schemes, respectively, will need to be investigated further in the future.
The data corresponding to the plots in this work can be found in Ref. [38].