QCD transition at the physical point, and its scaling window from twisted mass Wilson fermions

We study the scaling properties of the finite temperature QCD phase transition, for light quark masses ranging from the heavy quark regime to their physical values. The lattice results are obtained in the fixed scale approach from simulations of $N_f=2+1+1$ flavours of Wilson fermions at maximal twist. We identify an order parameter free from the the linear contributions in mass due to additive renormalization and regular terms in the Equation of State, which proves useful for the assessment of the hypothesized universal behaviour. We find compatibility with the 3D $O(4)$ universality class for the physical pion mass and temperatures $120$ MeV $ \lesssim T \lesssim 300$ MeV. We discuss violation of scaling at larger masses and a possible cross-over to mean field behaviour. The chiral extrapolation $T_0 = 134^{+6}_{-4}$ MeV of the pseudocritical temperature is robust against predictions of different universality classes and consistent with its estimate from the $O(4)$ Equation of State for the physical pion mass.


Introduction
The appearance of pseudo-Goldstone bosons in the spectrum signals the spontaneous breaking of chiral symmetry of strong interactions. It is well known that chiral symmetry is restored at high temperatures, and imprinting of this phenomenon for physical masses is visible in the behaviour of the order parameters and in the spectrum. If the symmetry restoring transition is continuous, the link between the genuine critical behaviour in the chiral limit and the observations at finite masses is made transparent by the universal Equation of State (EoS). This is true only within a limited region around criticality -the scaling window.
Despite substantial progress [1,2], the nature of the continuum limit and the extent of its scaling window are still open issues. Chiral and axial symmetries play a pivotal role here: in this study we make use, as in our previous work, of Wilson fermions at maximal twist [3][4][5][6][7], a lattice formulation with good chiral properties and an alternative to the more used staggered fermions.
The global symmetry of the QCD Lagrangian U(n) L ×U(n) R S U(n) × S U(n)×U(1) V × U(1) A , valid at classical level, is broken by topological quantum fluctuations. In the limit of an infinite strange mass there are at least two possible scenarios for the high temperature transition depending on the fate of the axial symmetry [8][9][10]: if the axial symmetry breaking is not much sensitive to the chiral restoration, the breaking pattern is S U(2) L × S U(2) R → S U(2) V . Due to the associate diverging Email addresses: a.kotov@fz-juelich.de (Andrey Yu. Kotov), lombardo@fi.infn.it (Maria Paola Lombardo), amtrnn@gmail.com (Anton Trunin) correlation length, the theory is effectively 3D, in the universality class of the classical four-component Heisenberg antiferromagnet: O(4) → O(3) [8]. If instead axial symmetry is correlated with chiral symmetry, the relevant breaking pattern would be U(2) L ×U(2) R → U(2) V , hinting either at a first or at a second order transition with different exponents [10].
We are interested in properties of the transition for physical values of the strange quark mass m s . If the two flavor transition is of first order, we are likely dealing with a first order transition for any m s . Alternatively, if the m s = 0 first order transition ends at m s = m s crit , m s merely renormalizes the coefficients of the effective action for m s > m s crit [8,9], without altering the critical behaviour. In this case one conventionally assumes that there is a line of second order transition T c (m l = 0, m s ) > T c (m l = 0, m s ) for m l = 0, ∞ > m s > m s crit . Away from the true critical point dimensional reduction may fail, resulting in a 4D theory and a mean field behaviour. The extent of the scaling window and the threshold of dimension reduction are non-universal features, which should be investigated by ab-initio methods.
In this study we define an ad-hoc order parameter, which we dub ψ ψ 3 , free from contributions linear in mass [48]. We explore a range of masses and temperatures, in an attempt to identify the limits of the scaling window, and its possible crossover to a mean field behaviour. A very preliminary account of some of the results has been presented in Ref. [3].

Observables, magnetic Equation of State, and scaling
A quantitative way to describe a critical system with a breaking field h relies on the EoS for the order parameter M (1) where m l is the light quark mass and T 0 ≡ T c (m l → 0) is the critical temperature in the chiral limit, δ and β are critical exponents. Note that there are two independent normalizations for M and for t, which, like T 0 , depend on the strange quark mass. In our range of masses it is legitimate to trade m l m 2 π , with a suitable adjustment of normalizations. f (x), with x = t/h 1/βδ , is a universal curve for a given breaking pattern.  For the 3D O(4) scaling function we rely on the parameterization from Ref. [49] as well as on its polynomial interpolation [50,51], which is smooth in the critical region. A successful scaling of results for a finite mass would identify the critical temperature in the chiral limit, besides confirming universality. One less stringent approach to scaling relies on pseudo-critical temperatures associated with features of the order parameter. Consider the two susceptibilities: χ L = ∂ ψ ψ ∂m l and χ ∆ = ∂ ψ ψ ∂T , which peak at t/h 1/βδ = 1.35 (3) and t/h 1/βδ = 0.74(4), respectively, for 3D O(4) universality class with the critical exponents δ = 4.8(1) and β = 038(1). The peak positions define pseudocritical temperatures A is a mass independent parameter, and T c for different observables should scale with the same exponent 2/βδ, but with different z p s.

A new order parameter
In the same spirit as Ref. [52], we consider the transverse χ T = ψ ψ m l and longitudinal χ L = ∂ ψ ψ ∂m l susceptibilities and their difference O LT ≡ χ T − χ L . O LT is an order parameter for the transition, since χ L = χ T in the chiral limit in the symmetric phase. Moreover, by construction, any linear contribution in m l -either coming from the regular part of the EoS or from additive renormalization -drops from the difference. However O LT is divergent in the chiral limit in the broken phase. The singularity is avoided by defining: This is the order parameter, whose Taylor expansion in m lwhen defined -starts at a third order, as noted in Ref. [48].  Deriving the Equation of State for ψ ψ 3 is straightforward:   The high temperature leading term is ψ ψ 3 ∝ t −γ−2βδ rather than ψ ψ ∝ t −γ : the decay is rather fast. The inflection point that drives the behaviour of the pseudo-critical temperature associated with ψ ψ 3 is x infl = 0.55 (1). Note that the pseudocritical temperature for this observable follows the same scaling in Eq. (2), but with a smaller z p : this means that the pseudocritical temperature associated with the inflection point of ψ ψ 3 is smaller than the pseudo-critical temperatures associated with chiral condensate and susceptibility. This also implies that it is closer to the true critical one in the limit m l → 0.
In Fig. 1 we compare the EoS for ψ ψ 3 with the one for ψ ψ for the O(4) Universality class and for the mean field. We note the sharper decrease of ψ ψ 3 , very understandable given that it is closer to the chiral condensate in the chiral limit, followed by the high temperature behaviour just described. For either observable we also show the mean field result: it is indeed very close to the 3D O(4), so the transition from the scaling window to a regime with small fluctuations is expected to be very smooth.

Numerical results
In this work we present new results for the physical pion mass, together with an extended analysis of the results for higher pion masses [5]. Simulations [53], renormalization factor Z p for these ensembles was measured in [54]. All parameters of the M ensemble are from [55]. The strange and charm quark masses are close to their physical values. approach, where we keep the bare lattice parameters fixed and vary the temperature by varying the number of lattice points in the temporal direction N t . Gauge fields ensembles used in this paper are summarized in Table 1. New simulations for the physical pion mass use Wilson-clover twisted mass fermions, with the parameters of the action corresponding to the zero temperature ensemble cB211.072.064 of ETMC [56]. The number of configurations generated for these ensembles is given in Table 2. First results of our simulations with physical pion mass were presented in [3]. The chiral condensate is an (approximate) order parameter    for the chiral symmetry in the light quark sector: In numerical simulations the trace of M −1 was measured using noisy stochastic estimator with 24 random volume sources. It is important to note that the bare chiral condensate (5) should be renormalized. However, in the fixed scale approach, applied in this paper, both additive and multiplicative renormalization factors are equal for all points. Thus, the renormalization procedure has no effect on the pseudo-critical temperatures extracted from the peak of χ ∆ , or, equivalently, from the inflection point of the chiral condensate.
In Fig. 2 we present the dependence of the chiral conden-sate (5) on the temperature for the physical pion mass. We extract the pseudo-critical temperature T ∆ from the inflection point of this dependence. For this purpose we fitted the chiral condensate with several functions, varying fitting interval: We find that the quality of the fit worsens when the upper limit approaches 300 MeV. The final estimation is obtained by averaging over three functional forms and two fitting intervals [120:210] and [120:250] MeV (apart from the polynomial fit, which describes the data poorly at large interval [120:250] MeV). The difference between various functions/intervals was used to estimate the systematic uncertainty. The final results for the pseudo-critical temperature T ∆ are quoted in Table 3. The results for heavier pion masses are from our previous analysis [5]. The chiral susceptibility, which is defined as the mass derivative of the chiral condensate χ L = ∂ ∂m l ψ ψ , consists of con- nected χ conn and disconnected χ disc contributions: Also, the chiral susceptibility suffers from additive and multiplicative renormalizations, which, again, are not affecting the estimate of the pseudo-critical point. In Fig. 3 we show the chiral susceptibility for all ensembles. We note that this observable has a strong regular contribution, in addition to an additive renormalization. These features are qualitatively clear in the plots: rather than the simple symmetric shape predicted by the EoS, the curves are skewed and have a long high temperature tail. For this reason, and given a small number of points in the transition region, we decided to use a model independent estimate of the pseudo-critical temperature, based on cubic spline interpolation, instead of using an explicit functional form. To estimate the statistical uncertainty we added random Gaussian noise to each point, with the amplitude given by the statistical uncertainty of our data points. Typically we used O(2000) splines for error estimate, and we find this a rather robust procedure. The dispersion of the different results can be appreciated from the plots, and we show the statistical errors in a few selected cases -the others behave similarly.
Finally we compute ψ ψ 3 as in (3). As discussed above in Section 2, ψ ψ 3 is free from linear additive renormalization as well as from linear correction to scaling. It still needs a multiplicative renormalization, which is obviously the same Z p as for the chiral condensate, Table 1.
The temperature dependence of ψ ψ 3 is shown in Fig. 4. We used several functional forms to fit the temperature depen-dence of this observable, which shows a fast fall-off with temperature. We then find that the logistic curve is able to capture the behaviour for a sizeable range, while polynomial fits are limited to a few points. Upon the fit we tried to keep the fit interval fixed in units T/T c : for logistic fit the interval is roughly [0.8, 1.5], for polynomial [0.8, 1.25]. We used both functional forms to estimate the pseudo-critical temperature and its systematic uncertainty. The final results for the pseudo-critical temperature extracted from the chiral susceptibility and for ψ ψ 3 for all ensembles are presented in Table 3.

The magnetic Equation of State and ψ ψ 3
In Fig. 5 we show the results for ψ ψ 3 as a function of temperature, with superimposed logistic fits. The results are converted to physical units using the lattice spacing and the multiplicative renormalization given in Table 1 and divided by m 2/δ π , with δ fixed at the O(4) value. At the critical temperature ψ ψ 3 ∝ m 2/δ π , hence ψ ψ 3 m 2/δ π should not depend on the mass at the critical point. The logic is similar to the one adopted in Refs. [33,52] when searching for the fixed point of χ L /χ T . The crossing point of the curves for the lightest masses identifies a candidate critical temperature in the chiral limit, T δ 0 = 138(2) MeV, where the superscript δ refers to T 0 being estimated from the crossing point. Possible scaling violations would result in a mass dependence of T δ 0 . Next, we fit to the 3D O(4) EoS with open critical temperature. We rely on a smoothing spline to define the EoS and we introduce appropriate (pion mass dependent) scaling parame-  (2) MeV, but there are obvious violations for larger masses. In the right-hand side of Fig. 6 we compare the mean field fits with the O(4) EoS: as it was already clear from Fig. 1, they are very close to each other. In the smaller interval they basically coincide, only for the larger interval shown here, and for the lightest mass, O(4) fits may be slightly favored.
Finally, in Fig. 7 we show high temperature fit of the condensate ψ ψ 3 , constrained to the O(4) behaviour: ψ ψ 3 ∝ t −γ−2βδ for T 0 = 138 MeV 1 . The results are rescaled by m 3 l m 6 π , the anticipated high temperature leading behaviour. In the interval of temperatures 160 MeV < T < 300 MeV (marked bold) the O(4) prediction fares nicely through the data. For T 300 MeV the behaviour is distinctly different, and the data collapse on a single curve, as expected of the leading mass   term according to Griffith analyticity, regardless the critical behaviour. This suggests that the temperature extent of the scaling window above T 0 extends up to about 300 MeV. Interestingly, in a previous study [5] we have found that this is also the threshold for a dilute instanton gas behaviour.

The scaling of the pseudo-critical temperatures and the chiral limit
In the summary Table 3 we present all pseudocritical temperatures. The results on the B ensembles help monitoring the finite spacing effects, and we confirm that they are small also for the new observables introduced here. A more detailed discussions on the continuum limit can be found in Ref. [5] as well as in ETMC papers, see e.g. Ref. [57], and references therein.
As discussed, we fit the pseudo-critical temperatures to a power law in the pion mass, see Eq. (2), where z p is tabulated for the different observables in Table 4, A is a universal constant depending on the normalization scales, and T 0 ≡ T c (m π → 0) corresponds to the critical temperature in the chiral limit. We kept 2/(βδ) = 1.083 fixed to the value given by O(4) critical exponents. In Table 4 we present T 0 and ratio of z p extracted from the fit as well as the prediction of O(4) scaling. There is a semi-quantitative agreement with O(4) when considering the relative slopes of ψ ψ and ψ ψ 3 . When considering only the  Figure 9: Pseudo-critical temperatures with their chiral extrapolations: comparison with the results from the HotQCD Collaboration [58], FASTSUM Collaboration [28,29], Wuppertal-Budapest Collaboration [59]. The purple diamond at m π = 0 marks the critical temperature [30] which compares well with our result T 0 = 134 +6 −4 MeV (light-green cross, slightly shifted for better readability). Ensemble 157.8 (7)   The average of the extrapolated values of T c in the chiral limit is compatible with the estimated from the crossing point T δ 0 = 138(2) MeV, and only slightly below the estimate from the EoS for physical pion mass. This last point is consistent with the physical pion mass being close to the scaling window. Note that the error in Eq. (7) contains also systematic uncertainty, coming from extrapolation of different observables.
Similarly to Ref. [7] we also consider the breaking pattern corresponding to an effective restoration of the axial symmetry U(2) L ×U(2) R → U(2) V , leading either to a second order transition with δ = 4.3(1), β = 0.40 (4) or of first order transition, in which case the pseudo-critical behaviour would be driven by an endpoint in the Z 2 universality class. In the first case the exponents are close to the O(4) ones and within our errors the results are indistinguishable. For the first order scenario we rely on the experience gathered with three-flavor QCD, where it was found that the mixing associated with the critical endpoint could well be low [60]. Hence, as in Ref. [7], we ignore the mixing and fit to T c (m π ) = T 0 + B(m 2 π − m 2 c ) 1/βδ (8) with 1/βδ = 0.64 for the Z 2 3D universality class [61]. Also in this case, the critical exponents for O(4) and Z 2 behaviour are close to each other, and it is very difficult to distinguish between these two scenarios. One can clearly see from Fig. 8 that the data can be also described with Z 2 behaviour and critical mass up to m c ∼ 100 MeV. If, artificially, we constrained m c = 0, the lines would be indistinguishable from O (4). All in all, we concur with other observations [1,33] that extrapolation alone does not suffice to discriminate among different scenarios. A possible way to distinguish between different universal behaviours would include simulations with lower pion masses [30]. On the positive side, the extrapolated values of the pseudo-critical temperatures T are robust (up to the critical pion mass m π cr ) against the choice of critical exponents.
In Fig. 9 we compare our estimations of the pseudo-critical temperature with the results obtained by staggered [58,59] and Wilson [28,29] fermions. The nice agreement of our results with others can be appreciated. We also superimpose O(4) fits to our pseudo-critical temperatures and present our estimation of T 0 in the chiral limit which compares well with the estimations of HotQCD Collaboration [30].

Summary
The chiral behaviour of the QCD transition is a difficult, much studied problem. Direct tests of universality are hampered by scaling violations, and by the similarity among different scenarios. Here we have tackled these issues by introducing an order parameter, whose behaviour is closer to the critical one, hence with a reduced contribution from scaling violating terms; and by exploring a larger parameter region, with the goal to follow the evolution from a scaling behaviour till a regular one, unrelated with criticality.
Particularly interesting is the behaviour with physical quark masses. We have measured the pseudo-critical temperature for three different observables, finding critical temperatures in agreement with staggered estimates for the chiral condensate and the chiral susceptibility. The new order parameter ψ ψ 3 gives a lower pseudo-critical temperature T ∆ 3 = 146.2(21)(1) MeV, in accordance with its being closer to the true critical behaviour. We have used three different strategies for the computation of the critical temperature in the chiral limit: the first one is based on the scaling of ψ ψ 3 at the critical temperature T 0 (m π = 0), the second one is based on the fit of ψ ψ 3 to the magnetic Equation of State, which should give the value of T 0 for each single pion mass, barring scaling violations, the last one is based on extrapolation of the pseudo-critical temperatures extracted from different observables to the chiral limit. A search for the anomalous scaling of ψ ψ 3 at the chiral transition picks the critical temperature of T 0 = 138(2) MeV. For all the masses the results are apparently well described by the universal EoS, however only for the physical pion mass the estimated critical temperature T 0 = 142(2) MeV approaches the results obtained by extrapolation, T 0 = 134 +6 −4 MeV. We note that the pseudo-critical temperature T ∆ 3 = 146.2(21)(1) MeV for a physical pion mass serves as an upper bound for the critical temperature. As a final result for the critical temperature we quote T 0 = 134 +6 −4 MeV. In Fig. 9 we present our results for the pseudo-critical temperatures, their chiral extrapolation and comparison to the results of other groups, with which we find a nice agreement.
For larger masses the 3D O(4) EoS fits are indistinguishable from mean field, while a mild tension appears for the physical mass, with O(4) slightly favoured. The high temperature behaviour shows a clear threshold at temperatures about 300 MeV, above which a trend consistent with O(4) scaling gives way to a simple leading order analytic scaling dictated by Griffith analyticity. As a side comment, in a previous study [5] we have found that this is also the threshold for a dilute instanton gas behaviour. A threshold around the same temperature was also observed in other lattice studies [62][63][64]. Our observations are compatible with an onset of scaling around the physical value of the pion mass, and temperatures below 300 MeV. In Fig. 10 we sketch the critical region, highlighting the compatibility with O(4) scaling. Our results seem to be in contradiction with the results of [21,46], where scaling is observed only for very tiny pion masses. Possible sources of this discrepancy and connection to the presented results are the subject of future work and will be discussed elsewhere.
We also studied the possibility of first order phase transition in the chiral limit, in which case the observables are expected to scale according to Z 2 universality class with some critical pion mass m c . Our data can be also described with Z 2 behaviour and critical mass up to m c ∼ 100 MeV. Thus we conclude that extrapolation alone does not suffice to discriminate among different scenarios. Similar conclusions were reached in the two flavor model [7]. A possible way to distinguish between different universal behaviour would include simulations with lower pion masses [30,33], hoping to observe directly a first order transition. While a direct observation of such a first order transition would settle the issue, a lack of observation makes it difficult to exclude that such a transition appears at very low masses. We hope that by exploiting more the universal properties of the Equation of State, possibly using some of the strategies outlined here, and in comparison with analytic approaches, may provide a more solid answer.
It would also be interesting to repeat a similar analysis using staggered results for pion masses smaller than physical [33] as well as with FRG results [21]. Also, Fig. 1 suggests that results below the critical temperature in the chiral limit may be useful to further discriminate true critical behaviour from mean field. This work should be ameliorated and extended along several lines: we have not performed a continuum limit extrapolation. The indirect checks performed by the ETMC collaboration as well as by ourselves [5], and the good consistency of the results at the larger mass, see Table 3, give some confidence that residual spacing corrections should not exceed a few per-cent. The fixed scale approach has several advantages, but we have to rely heavily on interpolations. Results on finer lattices should also help in this respect, by producing results at intermediate temperatures. Since the ETM Collaboration has recently released a new set of zero temperature tuned parameters [57], we hope to return to this point in the future. Finally, we have completed an analysis of correlators and screening masses on the configurations for physical pion mass, to complement the discussion on chiral and axial symmetry, and the results will appear soon [65].