Lattice Monte-Carlo study of pre-conformal dynamics in strongly flavoured QCD in the light of the chiral phase transition at finite temperature

We study the thermal phase transition in colour SU(3) Quantum Chromodynamics (QCD) with a variable number of fermions in the fundamental representation by using lattice Monte-Carlo simulations. We collect the (pseudo) critical couplings for N_f=(0, 4, 6,8), and we investigate the pre-conformal dynamics associated with the infra-red fixed point in terms of the N_f dependence of the transition temperature. We propose three independent estimates of the number of flavour N_f^* where the conformal phase would emerge, which give consistent results within the largish errors. We consider lines of fixed N_t in the space of (N_f, bare lattice coupling), and locate the vanishing of the step scaling function for N_f^*\sim 11.1\pm 1.6. We define a typical interaction strength (g_TC) at the scale of critical temperature T_c, and we find that g_TC meets the zero temperature critical couplings estimated by the two-loop Schwinger Dyson equation or the IRFP coupling in the four-loop beta-function at N_f^*\sim 12.5\pm 0.7. Further, we study the N_f dependences of T_c/M where M is a UV N_f independent reference scale determined by utilising the coupling at the scale of the lattice spacing. Then, T_c/M turns out to be a decreasing function of N_f, and the vanishing T_c/M indicates the emergence of the conformal window at N_f^* \sim 10.4 \pm 1.2.


Introduction
The analysis of the phases of strong interactions presents many fascinating aspects -mechanisms of confinement, different realisations of the chiral symmetry, the nature of the symmetric phase, the emergence of conformality, and many others. All these topics are under active scrutiny both theoretically and experimentally [1]. While strong interactions spontaneously break chiral symmetry in ordinary QCD at zero temperature, the chiral symmetry is realised either at high temperatures -in the so-called quark-gluon plasma (QGP) phase -and at a large number of flavours N f > N * f (even at zero temperature) [2,3,4,5]. In the latter case, the theory is expected to become not only chirally but also conformally invariant. This is due to the emergence of an infra-red fixed point (IRFP) for N f > N * f at a coupling which is not strong enough to break chiral symmetry. Both physics intuition and phenomenological analysis based on functional renormalisation group [6] and finite temperature holographic QCD [7] indicate that the conformal phase of cold, many flavour QCD and the high temperature chirally symmetric phase are continuously connected. In particular, the onset of the conformal window coincides with the vanishing of the transition temperature, and the conformal window appears as a zero temperature limit of a possibly strongly interacting QGP.
The analysis of the finite temperature phase transition is a well-established line of research within the lattice community, and our approach will be completely conventional here. According to the Pisarski-Wilczek scenario [8], the most likely possibility for N f ≥ 3 is a first order chiral transition in the chiral limit, turning into a crossover above a critical mass endpoint, and/or on lattices which are not large enough. We will identify such crossover with confidence for a number of flavours ranging from four to eight, and we will complement these results with those of the deconfinement transition in the quenched model. Then, we study the approach to the conformal phase in the light of the chiral phase transition at finite temperature with variable number of flavours.
One problem of this approach is the setting of a common scale among theories which are essentially different. We will propose two alternative possibilities to handle this problem, one evolving from our previous work [9], and the other from a recent analysis [10]. Interestingly, this latter approach analyses the dependence of the confinement critical couplings at each N f . In Section 5, we investigate the N f dependences of the chiral crossovers, and estimate the lower edge of the conformal window N * f . Finally in Section 6, we provide concluding remarks. The Appendix is devoted to the summary tables of the simulation parameters and the numerical results obtained by analysing the simulation outputs.

Simulations' setup
We investigate finite temperature QCD with different number of flavour N f = (0, 4, 6, 8) by utilising the publicly available MILC code [51]. The temperature T is related to the inverse of the lattice temporal extension, and we control it by varying β L at fixed N t . The number of lattice points in the spacial directions N s is chosen such that the aspect ratio N s /N t ≥ 2 in all our runs. For each N f , we use a single bare fermion mass ma = 0.02. The simulation parameters used in this study are summarised in Appendix A.

Action and algorithm
The setup for the action explained below is the same as the one used for N f = 8 in Ref. [28] up to the number of flavour. We use an improved version of the staggered action, the Asqtad action [52], with a one-loop Symanzik [53,54] and tadpole [55] improved gauge action, where g L is the lattice bare coupling, and β i are defined as The plaquette coupling β p = 10/g 2 L ≡ β L is a simulation input. The M[am, U, u 0 ] in Eq. (2) denotes the matrix for a single flavour Asqtad fermion with bare lattice mass am, and U C i represents the trace of the ordered product of link variables along C i , for the 1 ×1 plaquettes (i = p), the 1 ×2 and 2 ×1 rectangles (i = r), and the 1 ×1 ×1 parallelograms (i = pg), respectively -all divided by the number of colours. The tadpole factor u 0 is determined by performing zero 3 temperature simulations on the 12 4 lattice (the second column of Tables A.4 -A. 14), and used as an input for finite temperature simulations.
To generate configurations with mass degenerate dynamical flavours, we have used the rational hybrid Monte Carlo algorithm (RHMC) [56], which allows to simulate an arbitrary number of flavours N f through varying the number of pseudo-fermions. The quenched (N f = 0) system has been realised by using massive bare fermion mass ma = 1.0 in the four flavour system. The six flavour system has been realised by using two pseudo-fermions in the rational approximation with a quarter root technique, N f = 4 · 2 · 3/4 = 6. Then, we have assumed the rooting does not affect the results within the accuracy of our simulation. For the other number of flavour (N f = 0, 4, 8), we do not use the rooting.
We have adjusted the micro-canonical step length δτ and the step length of a single trajectory ∆τ = 20 × δτ to realise 75 − 80 percent Metropolis acceptances. Details are reported in the fourth column of Tables A.4 -A.14. For each parameter set, we have collected a number of trajectories ranging from a one thousand to ten thousands -the latter closer to the chiral crossover regime.

Observables
The focus of this paper is the analysis of the chiral phase transition The fundamental observable is then the order parameter for chiral symmetry, the chiral condensate: where N s (N t ) represents the number of lattice sites in the spatial (temporal) direction. We have measured a 3 ψ ψ by using a stochastic estimator with 20 repetitions. We have also measured connected and disconnected chiral susceptibilities, Here we have conveniently written the chiral condensate and its susceptibilities in terms of traces of (products of) the staggered fermion matrix M. We note that the MILC convention for the chiral condensate gives the twice of Eq. (5), as will be indicated several times in the following sections for results. We have measured the susceptibilities a 2 χ conn and a 2 χ disc separately. The disconnected chiral susceptibility is a non-local quantity which can be estimated from the variance of the bulk behaviour of the chiral condensate. Since we have used the stochastic estimator for the chiral condensate measurements, the variance would automatically include part of the connected contributions through random sources multiplying themselves. Following Bernard et al. [57], we take into account this effect in our estimate for the disconnected part a 2 χ disc by considering the only off-diagonal elements of the covariance matrix for the random sources.
The measurements of a 3 ψ ψ and a 2 χ conn,diss allow us to construct two physically relevant quantities: the scalar and pseudo-scalar susceptibilities, is a probe of the chiral symmetry [28,58]. This is owing to the fact that χ σ and χ π are related through Ward identities to the spacetime volume integral of the scalar (σ) and pseudo-scalar (π) propagators. In the chiral limit, the susceptibility ratio R π should be one in chirally symmetric regime due to the degeneracy of the chiral partners, while it should be zero in the spontaneously broken phase. Even including a finite bare fermion mass, R π still has a strong signal for the chiral transition or crossover. In particular, R π ∼ 1.0 in the chirally symmetric regime holds true till the chiral condensate is dominated by the linear mass term contribution. It turns out that R π allows the identification of a pseudocritical coupling β c L associated with the chiral crossover, which , in the cases we have studied, coincides in the error with the pseudo-critical coupling determined from the maximum of the chiral susceptibility.
In the gauge sector, we measure the Polyakov loop, where tr c denotes the trace in colour space, and U 4,tx is the temporal link variable. From the variance of L, we also evaluate the susceptibility for the Polyakov loops.

Statistics, and error analysis
In the vicinity of the chiral crossover, we have a long auto-correlation time, and thermalization checks require extra care. Here we explain our analyses by using a typical example: The left panel of Fig. 1 displays the evolution of the chiral condensate on the lattice volume 24 3 × 12 just before the chiral crossover β L = 5.45 and 5.50 at N f = 6, one of the most time-consuming examples in our simulations. (In order to shorten the simulation time, we started the evolution from thermalized configurations obtained at β L < 5.45.) We have computed the ensemble averages by using the last 2500 (2000) trajectories at β L = 5.45 (5.50) and we have confirmed that they are consistent with those obtained by using last 2000 (1500) trajectories. We have then used the latter trajectories to evaluate the average. In the cases we are considering, this corresponds to the data found in the right-hand side of the vertical green (dashed) lines in the left panel of Fig. 1.
We divide the obtained data set into several bins and utilise the jackknife method in order to take into account the auto-correlation effect in the error estimate. As a bin-size s bin becomes larger, the jackknife error increases (the right panel of Fig. 1), which is due to the decrease of the effective number of (uncorrelated) data (n ave /s bin , n ave = the number of trajectories to calculate the average). For a sufficiently large s bin , the jackknife errors at β L = 5.45 and 5.50 level off, giving a reliable error estimate.
Here is the result obtained from the above procedures: We have performed the analyses explained here for all the various β L , N f , and the lattice volumes. The results are summarised in Appendix A.

Results on the lattice thermal transition
In this section, we show our simulation results on the chiral and deconfinement crossover for the different number of flavours N f .
We have used a common bare fermion mass ma = 0.02 for all simulations at finite N f . According to the Pisarski-Wilczek scenario [8], the most likely possibility for N f ≥ 3 is a first order chiral transition in the chiral limit. Introducing a bare fermion mass, the first order phase transition will eventually turn into a crossover for masses larger than some critical mass. Since the chiral condensate looks smooth in our results, we are most likely above the critical endpoint in all the cases we have studied, and we use the terminology of "chiral crossover" in the following.
The finite bare mass ma = 0.02 might have a different physical relevance at each N f , as well as for different bare coupling for a fixed N f . It remains then to be seen how our results would change in the chiral limit, and we hope to come back to this point in a future study. Since we have noted that at strong coupling and small masses the improvement term in the Action might be responsible for the spurious phases [23,22] observed also in Ref. [32,25], we might also consider an unimproved action for this study. Before entering into details, let us summarise our main results, i.e. the critical lattice couplings β c L associated with the chiral crossover in Table 1. For N f = (4, 6, 8) we have observed that the peak position of the chiral susceptibility a 2 χ σ , whenever clearly defined, coincides within the errors with the inflection point of R π defined in Eq. (9), as well as with the inflection point of the chiral condensate and that of the Polyakov loop. This indicates that the crossover region is rather narrow, as different indicators give consistent pseudo-critical points. We then quote the common pseudo-scalar coupling, with a conservative error estimate. For the quenched (N f = 0) case, we have extracted the pseudo-critical coupling from the deconfinement crossover by evaluating the peak position of the Polyakov loop susceptibility.
In the following subsections, we present these results in detail, starting from N f = 6 and 8 in the first two subsections, and continuing with the N f = 4 and N f = 0. The reader who is not interested in these technical details is advised to skip the rest of this Section and proceed directly to the next one.

Chiral crossover at N f = 6
We show the N f = 6 results for a fixed bare fermion mass ma = 0.02. In Figs. 2, the chiral condensate a 3 ψ ψ (PBP, red ) the real part of Polyakov loop L (Re[PLOOP], blue ), the chiral susceptibility (χ σ , red +), and the chiral susceptibility ratio (R π , blue ×) are displayed as a function of a lattice coupling β L = 10/g 2 L . The first, second, third, and fourth lines in the figure show the results obtained by using temporal extensions N t = 4, 6, 8, and 12, respectively. We shall now extract the critical lattice couplings β c L associated with the thermal chiral crossover from these results.
As shown in the left panel of the first line in Fig. 2, the largest decrease of chiral condensates (PBP, red ) (as well as a drastic increase of the real part of Polyakov loops (Re[PLOOP], blue )) is found between β L = 4.65 and 4.70. Thus, we expect the chiral crossover in this region. As shown in the right panel of the first line in Fig. 2, the chiral susceptibility a 2 χ σ (red +) has a clear peak at β L = 4.65. In order to have a practical and coherent procedure to estimate the maximum, we have performed Gaussian fits: The Gaussian fit for the susceptibilities in the range [4.4, 4.9] leads to a maximum at a slightly larger β L (red dashed line). Further, the susceptibility ratio R π (blue ×) has an inflection point around β L = 4.65 − 4.70. For larger β L , the increasing rate of R π significantly reduces, and gets to almost unity. Thus, all observables consistently indicate the pseudo-critical coupling to be β c L = 4.675 ± 0.05 for (N f , N t ) = (6, 4). The error is determined to include the next-to-neighbour data and the maximum of the Gaussian fit.
The second line in Fig. 2 displays the results for N t = 6. As shown in the left panel, the largest decrease of chiral condensates (PBP, red ) (as well as a drastic increase of the real part of Polyakov loops (Re[PLOOP], blue )) is found between β L = 5.00 and 5.05, and we expect the chiral crossover in this region. As shown in the right panel, the chiral susceptibility a 2 χ σ (red +) has a peak at β L = 5.05, and the Gaussian fit for the susceptibilities in whole range of β L has a maximum at a slightly smaller β L = 5.0 (red dashed line). The susceptibility ratio R π (blue ×) has an inflection point around β L = 5.00 − 5.05, and then, it goes into the plateau domain. All observables consistently indicate the pseudo-critical coupling to be β c L = 5.025 ± 0.05 for (N f , N t ) = (6, 6). The error is determined to include both β L = 5.0 and 5.05 enough.
The third line in Fig. 2 shows the results for N t = 8. As indicated by the left panel, the chiral condensates as well as the Polyakov loops look smooth at almost everywhere, and it is difficult to locate the crossover point from them. As shown in the right panel, the chiral susceptibility a 2 χ σ (red +) has a peak at β L = 5.2, and the Gaussian fit for the susceptibilities in whole range of β L has a maximum at a slightly smaller β L = 5.17 (red dashed line). The  , blue ), and the right panel displays the chiral susceptibility (χ σ , red +) and the chiral susceptibility ratio (R π , blue ×), as a function of β L . For N t = 4, the Gaussian fit for the chiral susceptibility has been performed in the range [4.4, 4.9] to capture the peak structure.  susceptibility ratio R π (blue ×) exhibits the largest variation between β L = 5.15 and 5.2, after which the increasing rate of R π reduces and eventually evolves into almost unity. From the peak position of a 2 χ σ , we estimate the critical coupling to be around β c L = 5.20 ± 0.05 for (N f , N t ) = (6, 8). The error is determined to include the next neighbour data, the maximum of the Gaussian fit for χ σ , and the R π inflection point.
Finally, we analyse the results for N t = 12, the largest temporal extension in our N f = 6 simulations. The β L dependences of chiral condensates are found to be particularly smooth for whole range of β L = 4.7 − 5.7. Note that in this case, the aspect ratio (N s /N t ) is only two, and larger volumes would be required to reach a comparable clarity in the signal. As shown in the left panel of final line in Fig. 2, the onset for the Polyakov loop at β L = 5.525 is still appreciable (blue ). We here notice that the increase of Polyakov loops so far has been found just before the chiral crossover in the case of smaller temporal extensions N t = 4 − 8, though the Polyakov loop itself is not associated with the chiral dynamics. Based on such an experience, we assume that the chiral crossover at N t = 12 is in the vicinity of the onset of the Polyakov loop, and carefully investigate the corresponding region β L = 5.35 − 5.60. The chiral condensates do not have any clear signal (red in the left panel). As shown in the right panel, the chiral susceptibility a 2 χ σ (red +) has a small peak-like structure at β L = 5.575, and a bump-like structure at β L = 5.50. The chiral susceptibility ratio R π (blue ×) has the largest increase between β L = 5.55 and 5.575, and tends to be flat in β L ≥ 5.575. Thus, the critical lattice coupling would be in the range 5.50 ≤ β c L ≤ 5.575. Here, we employ a conservative estimate β c L = 5.55 ± 0.1, which sufficiently covers the whole candidate range. Our β c L collection at N f = 6 is found in the third line of Table 1. As will be shown in the next subsection, the N t dependent nature of β c L at N f = 6 (a thermal scaling) is associated with the uniqueness of the physical critical temperature, indicating the chiral (non-conformal) dynamics at N f = 6.

Chiral crossover at N f = 8
In our previous paper [28], we have studied the chiral phase transition at N f = 8 by using two lattice temporal extensions: N t = 6 and 12. One of the main results was that the chiral phase transition at N f = 8 still showed a thermal scaling property, which indicated the existence of a typical scale associated with the chiral dynamics rather than the conformality. We here add additional data computed at N t = 8, and confirm the thermal scaling at N f = 8, for this largish mass.
The left panel of Fig. 3 shows ensemble averages of chiral condensates a 3 ψ ψ (PBP, red ), the real part of Polyakov loop L (Re[PLOOP], blue ) as a function of β L . We observe the largest decrease of the chiral condensate between β L = 4.25 and 4.30, while the real part of the Polyakov loop stars growing around β L = 4.25. Although the error is huge, the chiral susceptibility ratio R π seems to have a larger increase between β L = 4.25 and 4.30. The large error of R π at β L = 4.25 comes from a very long auto-correlation in the Monte Carlo trajectories, which would have required a much larger statistics. The long correlation hints at a criticality. All observations consistently indicate the critical coupling to be β c L = 4.275 ± 0.05. Combining with N t = 6 and 12 data [28], we summarise the critical coupling β c L at N f = 8 in the final line of Table 1. 8 Here we should put some caveats on the N f = 8 results: First, we have not observed a peak-like structure in the chiral susceptibility χ σ . We should probably study larger spatial volumes, with similar aspect ratio as those studied in other cases. The location of the pseudo-critical point might change. For the time being, we rely on the experiences with the other systems, and infer the pseudo-critical coupling from the other observables. Second, even in the strong coupling region β < β c L , R π shows relatively large value ∼ 0.7 − 0.8. We should go to even stronger coupling and smaller masses before observing a clear mass gap. And third, the N t dependence of β c L shows a large deviation from the two-loop asymptotic scaling law as will be shown in the next subsection. This is not surprising, given that the couplings explored for N f = 8 are larger than in other cases. This could imply something which cannot be captured at two-loop, for example, the pre-conformal dynamics. Apparently, these caveats call for more detailed and quantitative lattice studies with a larger lattice size and a smaller bare fermion mass before drawing definite conclusions on N f = 8. We note a recent study claiming the conformality emerges for N f = 8 for small enough quark masses [31].

Chiral crossover at N f = 4
In Figs. 4, we show the chiral condensate a 3 ψ ψ (PBP, red ), the real part of Polyakov loop L (PLOOP, blue ), the chiral susceptibility (χ σ , red +), and the chiral susceptibility ratio (R π , blue ×) are displayed as a function of a lattice coupling β L = 10/g 2 L . The first, second, and third lines in the figure show the results obtained by using temporal extensions N t = 4, 6, and 8, respectively.
As shown in the left panel of the first line in Fig. 4, the largest decrease of chiral condensates (PBP, red ) (as well as a drastic increase of the real part of Polyakov loops (Re[PLOOP], blue )) is found between β L = 5.60 and 5.65, and we expect the chiral crossover in this region. As shown in the right panel of the first line in Fig. 4, the chiral susceptibility a 2 χ σ (red +) gets to a maximum at β L = 5.65. The Gaussian fit for the susceptibilities in the range [5.4, 5.8] leads to a maximum at a slightly larger β L (red dashed line). Further, the susceptibility ratio R π (blue ×) has an inflection point around β L = 5.60 − 5.65. For larger β L , the increasing rate of R π significantly reduces, and eventually evolves into almost unity. Thus, all observables consistently indicate the pseudo-critical coupling to be . The error is determined to include the next-to-neighbour data and the maximum of the Gaussian fit.
The second line in Fig. 4 displays the results for N t = 6. As shown in the left panel, the chiral condensates (PBP, red ) are found to be smooth, and it is difficult to locate the chiral crossover. The real part of Polyakov loops (Re[PLOOP], blue ) starts increasing around β L = 5.95. Based on our previous experiences, the chiral crossover could be around this region. As shown in the right panel, the chiral susceptibility a 2 χ σ (red +) has a maximum at β L = 6.00, and the Gaussian fit for the susceptibilities in the range of [5.8, 6.2] has a maximum at β L = 6.00 (red dashed line). From the maximum position of the chiral susceptibilities, we estimate the pseudo-critical coupling to be β c L = 6.00±0.05 for (N f , N t ) = (4, 6). The error is determined to include the next-to-neighbour data. The susceptibility ratio R π (blue ×) has a significant increase in β L > 5.9, and then, it flattens at β L = 6.1. This behaviour would be consistent to the above estimate β c L = 6.00. The third line in Fig. 4 shows the results for N t = 8. As indicated by the left panel, the chiral condensates as well as the Polyakov loops look smooth at almost everywhere, and it is difficult to locate the crossover point from them. As shown in the right panel, the chiral susceptibility a 2 χ σ (red +) has a peak at β L = 5.95. Indeed the susceptibility ratio R π (blue ×) also shows a bump structure around β L = 5.95, and implies some kinds of an instability of the system. However, the value of R π at β L = 5.95 turns out to be at most 0.4, which indicates a large remaining of the chiral symmetry breaking. In turn, R π keeps increasing till it approximately reaches to unity at β L = 6.30, and thus, the first peak at β L = 5.95 would not well capture the position of the chiral crossover. We here postpone the precise determination of the chiral crossover, and just provide a rough estimate of the pseudo-critical coupling: The susceptibility ratio R π has a large increasing rate in the range [6.0, 6.3]. We adopt the intermediate value as the pseudocritical coupling with the error covering whole range of [6.0, 6.3], β c L = 6.15 ± 0.15. This also includes the maximum of the Gaussian fit for the chiral susceptibility β L = 6.04.
The second line of Table 1 provides a summary of β c L for N f = 4.

Deconfinement at N f = 0
In this subsection, we estimate the critical lattice coupling β c L for deconfinement in the quenched (N f = 0) system. We note that both deconfinement and chiral transitions are associated with the thermal phase transition from the hadronic phase to the non-Abelian plasma phase with a drastic increase of the pressure (degrees of freedom). Then, our interest is the probe of the system with various N f in light of such a thermal phase transition. In this sense, we regard the deconfinement crossover at N f = 0 as a continuation to the chiral crossover at finite N f . In our setup, this connection is made explicit by the fact that we are realising a quenched system by use of a large mass in the four flavour system. It should be noted, anyway, that our result for the (pre-)conformal dynamics does not crucially depend on the quenched data. In Figs. 5, the thermalized ensemble averages of the absolute of Polyakov loop (|L|, blue ), its susceptibility (χ |L| , red +), are displayed as a function of a lattice coupling β L = 10/g 2 L . The first, second, and third lines in the figure show the results obtained by using temporal extensions N t = 4, 6, and 8, respectively.
As shown in the left panel of the first line in Fig. 5, the largest increase of the absolute value of Polyakov loops (|L|, blue ) is found between β L = 7.30 and 7.35, and we expect the deconfinement crossover in this region. As shown in the right panel, the susceptibility for |L| (red +) has a clear peak at β L = 7.35, hence we estimate the pseudo-critical coupling to be β c L = 7.35 ± 0.05 for (N f , N t ) = (0, 4). The error is determined to include the next-to-neighbour data. The second line in Fig. 5 displays the results for N t = 6. As shown in the left panel, the largest increase of the absolute value of Polyakov loops (|L|, blue ) is found between β L = 7.80 and 7.90, and we expect the deconfinement crossover in this region. The maximum of the susceptibility evaluated from |L| (red +) is observed at N f = 7.9. The large error indicates the long correlation time of the Monte Carlo trajectories. The Gaussian fit for the susceptibility has a maximum at β L = 7.97. From this, we estimate the pseudo-critical coupling to be β c L = 7.97 ± 0.07 for (N f , N t ) = (0, 6). The error is determined to include the next-to-next-to-neighbour data from the maximum point.
The third line in Fig . The error is determined to include the next-to-next-to-neighbour data from the maximum point.
The first line of Table 1 provides a summary of β c L for N f = 0. The N t dependent nature of β c L reflects the thermal nature of the crossover.

Asymptotic scaling analyses for chiral phase transition
In this section, we investigate the asymptotic scaling of the pseudo-critical temperatures T c where β c L have been computed in the previous section, and discuss the connection to the continuum physics. In the first subsection 4.1, we introduce the normalised critical temperature T c /Λ L/E (see e.g. [60]) where Λ L (Λ E ) represents the lattice (E-scheme) Lambda-parameter defined in the two-loop perturbation theory with or without a renormalisation group inspired improvement [61]. Then in the subsections 4.2, the asymptotic scaling will be assessed by studying the N t (in)dependence of T c /Λ L/E for each N f .

Normalised critical temperature T c /Λ L/E
We consider the two-loop beta-function The coupling g can be either the lattice bare coupling g L = 10/β L or the E-scheme renormalised coupling g E = 3(1 − P (g L )), where P (g L ) is the zero temperature plaquette value. If the one-loop perturbation theory exactly holds, the E-scheme coincides the lattice scheme.
Integrating Eq. (12), we obtain the well-known two-loop asymptotic scaling relation, where Λ L (Λ E ) is the Lattice (E-scheme) Lambda-parameter. To take into account higher order corrections, we have also considered the renormalisation group inspired improvement [61] where β L/E = 10/(g L/E ) 2 . The coupling β 0 can be arbitrarily set and the parameter h is adjusted so as to minimise the scaling violation. Note that h = 0 reproduces the standard asymptotic scaling law Eq. (15). The asymptotic scaling as described above is valid in the massless limit. In the following, we will use it to analyse results obtained at finite bare fermion mass ma = 0.02 by assuming that the shift of the (pseudo) critical coupling induced by a non-zero mass is smaller than other errors. This assumption should ultimately be tested in future studies by performing simulations with different masses and extrapolating to the chiral limit.
We now substitute β c L/E into the temperature definition Eq. (1), and insert the scale Λ L/E : The left-hand side is a given number, and Λ L/E a(β c L/E ) in the right-hand side is evaluated by using Eq. (15) and our critical couplings β c L/E . Thus, Eq. (17) allows us to convert the critical couplings into the (normalised) critical temperature T c /Λ L/E . When we adopt the improvement Eq. (16),    where g L/E denotes either the bare lattice coupling or the coupling defined in the E scheme. In addition, we consider the renormalisation group inspired definition,

Results for T c /Λ L/E and T c /Λ
where R imp is given by Eq. (16). The numerical results for T c /Λ L/E and T c /Λ imp L/E are collected in Table 2 and Table 3.

N f = 6
The left panel of Fig. 6 shows T c /Λ L as a function of N t for N f = 6, without ( red (×) symbols) and with (blue ( ) symbols) improvement. The improvement parameter h = 0.03 is adjusted to minimise the N t dependence. We have checked that the results are stable against small variation of h. Moreover β 0 is adjusted to match the results at β 0 = β c L (N t = 12) = 5.55. Similarly, the right panel of Fig. 6 shows T c /Λ E , which is nearly constant with N t . The overall behaviour suggests that the residual scaling violations at N t = 12 are small. Table 3: Summary of T c /Λ E and T c /Λ imp L/E for N f = 6 and N f = 8. The first (second) line at fixed (N f , N t ) shows the value of T c /Λ E (T c /Λ imp L/E ), and the last two columns give the parameter h and β 0 appeared in the improved asymptotic scaling Eq. (16). For N f = 6, the improvement was not necessary.  Fig. 7, we show the N t dependence of the normalised critical temperature in the lattice scheme, without ( red (×) symbols) and with (blue ( ) symbols) improvement. Again we adjust the improvement parameter h so to minimise the N t dependence, and we find that a larger h ≃ 1.08 is needed. This is consistent with the larger scaling violation observed between N t = 6 and N t = 12. Similar observations can be made on the E-scheme, where the scaling violations, as seen in the N t dependence, appear to be larger than those in N f = 6 case. Introducing the improvement in the E-scheme again requires a largish h = 0.4 (β 0 = β c E (N f = 8, N t = 12) = 5.94). In short summary, the system with N f = 8 shows much larger and less controllable deviations from the two-loop asymptotic scaling, than the ones observed for N f = 6 (and for N f = 4 and N f = 0, to be discussed next). This might be natural, in view of the largish values of the coupling involved. These observations confirm that in this case the β function at two loops cannot offer a quantitative guidance to strongly coupled pre-conformal dynamics.

N f = 0 and N f = 4
For N f = 0, we find about 10 percent variation of T c /Λ L in the whole range N t ∈ [4,8] (In this case T c represents the critical temperature associated with the deconfinement transition rather than the chiral transition). The N t dependence can be reduced to less than 10 percent by using the improved scaling with a small h = 0.05. Turning to N f = 4 we find about 20 percent variations of T c /Λ L between the N t = 4 and the N t = 8 results. The improved asymptotic scaling Eq. (16) works well, and the variation reduces into 10 percent level in whole range of N t = 4 − 8.

Scale separation?
In summary, T c /Λ computed using different schemes (Λ = Λ L or Λ E ) consistently shows an increase with N f , confirming and extending the findings of our early work [9]. As discussed in [9] this indicates that Λ L/E vanishes faster than T c upon approaching the critical number of flavour. Within the various uncertainties discussed here, this can be taken as a qualitative indication of scale separation close to the critical number of flavors.

Onset of the conformal window
In this section, we study the emergence of the conformal phase with increasing N f . In the first Subsection 5.1, we consider the phase diagram in the space spanned by the bare coupling g L and the number of flavor N f . We consider the (pseudo)critical thermal lines which connect the lattice (pseudo)critical couplings for a fixed N t . We argue that the critical number of flavor N * f can be identified with the crossing point of the pseudo-critical thermal lines obtained for various N t 's.
In the second Subsection 5.2, we introduce the thermal critical coupling g c T (N f ) as a typical interaction strength at the scale of the critical temperature T c (N f ). Since T c approaches zero when the number of flavor approaches the lower edge of conformal window N * f , g c T (N f = N * f ) should be equal to the zero temperature critical coupling g c (to be specified and estimated in the following). The relation g c T (N * f ) = g c thus defines implicitly the critical number of flavor N * f . In the final Subsection 5.3, we develop an improved version of the approach used in our early paper [9]. We introduce a UV N f independent reference scale, and compute the critical temperature for each N f in units of this scale. The infra-red dynamics affecting the critical temperature T c is then clearly exposed, and N * f can be estimated from the vanishing of T c .
Before turning to the details, we summarise our results for N * f : 11.1 ± 1.6 (from the lattice thermal lines) , 12.5 ± 1.6 (from the strength of the coupling at T c ) , 10.4 ± 1.2 (from the vanishing of T c ) . (20)

The critical number of flavor from the lattice thermal lines
We plot the lattice critical couplings g c L (N f , N t ) = 10/β c L (N f , N t ) ( Table 1) in the space spanned by the bare coupling g L and the number of flavor N f , and we consider the lines which connect g c L with N t fixed: g c L (N f )| N t =fix (see Fig. 8). These pseudo-critical thermal lines separate a phase where chiral symmetry (approximately) holds from a phase where chiral symmetry is spontaneously broken 1 . The resultant phase diagram may be seen as an extension of the well-known Miransky-Yamawaki phase diagram [3] to finite temperature.
We here argue that the critical number of flavor N * f can be read off from the crossing point of thermal lines obtained for different N t . To see this, we consider the well-known step-scaling function: where β L and β ′ L give the same physical scale ξ: Here,ξ is the dimension-less lattice correlation length, andξ/ξ ′ = s. In our case, ξ = T −1 c ,ξ = N t ,ξ ′ = N ′ t , and the above relation Eq. (22) reads As discussed in the previous study Ref. [34], ∆β s L = 0 holds at the IRFP regardless the scale factor s. In principle, we could then compute the step-scaling function from our numerical results, and try to see where it vanishes. Alternatively, we can look for the intersection of pseudo-critical thermal lines: obviously, ∆β s L = 0 holds at the intersection point regardless the value of the scale factor s.
To demonstrate this procedure, we consider the pseudo-critical lines obtained for N t = 6 and N t = 12 as shown in Fig. 8. Note their positive slope: the lattice critical coupling g c L is an increasing function of N f . This is a consequence of enhanced fermionic screening for a large number of flavor, as noted first in Ref. [62]. Interestingly, the slope decreases with increasing N t , which allows for a crossing point at a larger N f . Thus, we estimate the intersection at (g c L , N * f ) = (1.79 ± 0.12, 11.1 ± 1.6). We underscore at this stage that the above analysis is merely a qualitative discussion: the precise shape of the pseudo-critical thermal lines with fixed N t is dictated by the beta-function, which is unknown. Hence, a linear extrapolation which only uses two values of N t has only the meaning to illustrate a viable strategy which we plan to further pursue in the future. This caveat issued, the agreement of the results found here with the estimates presented below is rather gratifying.

The critical number of flavor from the interaction strength at T c
In this second subsection, we will follow the approach of a recent paper [10], and compute the coupling g c T (N f ) at the scale of the critical temperature for each N f . L at N f = 6 and 8, and considered "constant N t " lines with N t = 6, 12. If the system is still described by one parameter beta-function in this range of coupling, the IRFP could be located at the intersection of the fixed N t lines -or equivalently, in the region where the step-scaling function vanishes. To demonstrate the procedure -as a preliminary example -we have considered the intersection of the N t = 12 and N f = 6 lines.
To obtain the coupling at the scale of the temperature, we evolve the coupling at the scale of the lattice spacing a up to the temperature inverse scale N t a. To this end, we make use of two-loop expressions. Consider the renormalisation flow:R Since we are interested in the thermal critical coupling g c T , we set the reference mass to be the critical temperature itself: M(g ref L ) = 1/N t a(g c L ). Inserting this into Eq. (24), we see that g c T is implicitly given as: where we use the following g c L from Table 1: Alternative choices corresponding to smaller N t produce results for g c T suffering from the modest scaling violations discussed above.
The red ( ) symbol in Fig. 9 shows g c T as a function of N f . We superimpose a fit obtained by using the ansatz proposed in Ref. [10]  with A and B fit parameters, which describes well the data.
Since the critical temperature is zero in the conformal phase, the thermal critical coupling g c T should equal a zero temperature critical coupling g c when N f = N * f . Of course, g c is not known exactly and we have to rely on approximations.
The first estimate is based on the best available value g c SD obtained by using the two-loop Schwinger-Dyson equation [50]. In this case, the lower edge of the conformal window N * f is defined by the condition g c T (N * f ) = g c SD (N * f ). In Fig. 9 g c SD is plotted as a blue solid line. We then estimate the intersection of g c T and g c SD -hence the onset of the conformal window as well as the IRFP coupling at N * f -at (g * , N * f ) = (2.79, 13.2) ± (0.13, 0.6). One second possibility for estimating N * f is the following: the conformal phase would emerge when the coupling at IRFP (g IRFP ) is not strong enough to break chiral symmetry, i.e. g IRFP ≤ g c T . Here, we utilise the four-loop result for g IRFP 4l [49] as the best available. In Fig. 9, we show g IRFP 4l as magenta , with superimposed a linear interpolation. In the plot, we use the results for g IRFP 4l in theMS scheme. The errors are estimated by considering the scheme dependence [49], which turns out to be rather mild at four loops. We can then locate the intersection of g c T and g IRFP 4l and obtain (g * , N * f ) = (2.51, 11.8) ± (0.15, 0.9). Ideally, the three lines in Fig. 9 should meet at a (single) IRFP fixed point, if all the quoted results -including the analytic ones -were exact. Indeed the intersections we have estimated are consistent within the largish errors. We then quote the average of the above two estimates as our final result from this analysis, N * f ∼ 12.5 ± 1.6. In addition, and on a slightly different aspect, we note that g c T is an increasing function of N f . This indicates that the quark-gluon plasma is more strongly coupled at larger N f , as discussed in Ref. [10].

The critical number of flavor and the vanishing critical temperature
Finally, we present an estimate of the onset of the conformal window which closely follows our previous approach [9], based on the analysis of the N f dependence of the pseudo-critical temperature in units of a UV dominated scale. In this subsection, we introduce and exploit a new UV reference scale M.
Before going to details, we first explain the basic idea which follows the FRG analysis by Braun and Gies [6]. They used the τ lepton mass m τ = 1.777 (GeV) as an N f independent UV reference scale for theories with any number of flavours. The initial condition of the renormalisation flow has been specified via the strong coupling constant in an Starting from the common initial condition Eq. (28), the N f dependence of the critical temperature T c (N f ) emerges from the N f dependent renormalisation flow at the chiral phase transition scale µ ∼ Λ QCD ≪ m τ . The N f dependence of T c as well as its novel non-analytic behaviour in the pre-conformal region becomes free from the choice of the reference scale [6] by using an N f independent UV reference scale much larger than T c . Ideally, we would like to set our UV scale by measuring on the lattice some physical quantity insensitive to IR dynamics -for instance by fixing the value of α s in the V-scheme to some appropriate value, as done in the computation of r 0 or variations thereof, following Ref. [59] and related applications. These large scale simulations are now starting [65]. Here we design a simplified procedure.
In order to determine the reference coupling g ref L which appears in Eq. (24), we utilise our plaquette results P (equivalently, the tadpole factor u 0 = P 1/4 ) shown in Fig. 10.
Let us consider a constant u 0 , for instance u 0 = 0.8 in figure, and read off the corresponding bare lattice couplings at each N f . The obtained g L (N f ) is used as a reference coupling g ref L and the corresponding mass scale M(g ref L ) is again computed according to Eq. (24).
Some remarks on the aforementioned scale setting are in order: First, we recall the scale setting procedure in the potential scheme, where the measured normalised force r 2 F(r) is proportional to the renormalised couplingḡ, and the specificationḡ 2 ∝ r 2 X F(r X ) = ∃ X sets a scale r −1 X . In short, we use our u 0 (or equivalently plaquettes) to defineḡ, and u 0 = X is regarded as the analog of the potential scheme scale setting. Second, in the leading order of the perturbative expansion, the renormalised coupling is N f independent, and proportional to the Wilson loop [63] -a property that we have already exploited in the E-scheme calculation. Hence the use of an N f independent u 0 approximately gives an N f independent scale setting, similarly to the FRG scale setting method Eq. (28). And third, such an N f independent scale setting can be performed in a sufficiently UV regime T c (N f ) ≪ M(g ref L ) by adjusting the value of u 0 to satisfy the condition g ref L ≪ g c T ( ∀ N f ). Note that the coupling at the lattice cutoff a −1 (g c L ) is N t ≫ 1 times larger than T c . Then, the scale hierarchy T c (N f ) ≪ a −1 (g c L (N f )) allows us to consider a reference scale much larger than critical temperature but smaller than the lattice cutoff T c (N f ) ≪ M(g ref L ) < a −1 (g c L (N f )). We find that u 0 ∼ 0.8 meets this requirement.
where b 0,1 has been defined in Eqs. (13) and (14) We have further investigated the stability against different choices of u 0 : As shown in the right panel of Fig. 11, N * f is relatively stable within the range 0.79 ≤ u 0 ≤ 0.94. The scale cannot be pushed further towards the UV because of discritization errors. On the other hand, a small u 0 0.7 leads to M(g ref L ) ∼ T c or smaller. In such a case, the reference scale M(g ref L ) is affected by infra-red physics and cannot be used to study the vanishing of T c . Despite these limitations, the window of relative stability is however reasonably large, and suffices to define an average value for N * f . We quote the average among the three results obtained for u 0 = (0.79, 0.80, 0.81), i.e. N * f = 10.4 ± 1.2.

Summary
We have investigated the phase transition (crossover) at finite temperature T in colour SU(3) QCD-like theories with various number of flavours N f = 0 (quenched), 4, 6, and 8 by using lattice Monte Carlo simulations. We have used a single bare fermion mass ma = 0.02. For all number of flavours, we have used the Asqtad action with a oneloop Symanzik and tadpole improved gauge action. The main focus in this paper is to investigate the chiral crossover at finite T as a function of N f , and discuss the possible implication for the (pre-)conformal dynamics at large N f . In Eq. (20), we provide the summary for the number of flavour N * f where the conformal window would emerge. The observables in our simulations were the chiral condensate, the Polyakov loop, and their susceptibilities for various lattice couplings β L , lattice sizes, and the number of flavours N f . We have collected the (pseudo) critical lattice coupling β c L as a function of (N f , N t ). Table 1 provides the summary for β c L . Our β c L are consistent with enhanced fermionic screening at larger N f . The use of several N t allows us to study the asymptotic scaling of the critical temperature. Further, by utilising β c L , we have discussed a possible implication for the (pre-)conformal dynamics at large N f . We have estimated the N * f from the vanishing thermal scaling by extrapolating our critical couplings g c L to larger N f . This gives N * f ∼ 11.1 ± 1.6. We have extracted a typical interaction strength g c T at the scale of critical temperature T c by utilising our g c L and the two-loop beta-function, and compared g c T to the zero temperature critical couplings (g c SD ) estimated by the two-loop Schwinger-Dyson equation [50] as well as the IRFP position (g IRFP 4l ) of the four-loop beta-function [49]. The coincidence between g c T and g c SD or g IRFP 4l indicates the vanishing critical temperature with the emergence of the conformal phase. Based on this reasoning, we have estimated the onset of the conformal window as N * f ∼ 12.5 ± 1.6. We have also confirmed the increasing of g c T at larger N f which has been discussed in Ref. [10] and indicates more strongly interacting non-Abelian plasma at larger N f . Further, we have examined the N f dependence of T c /M by introducing a UV N f independent reference scale M which is determined by utilising the tadpole factor u 0 in analogous ways to the potential scheme scale setting. Then, T c /M turns out to be a decreasing function of N f consistently to the FRG observations [6], and the vanishing T c /M indicates the emergence of the conformal window around N * f ∼ 10.4 ± 1.2. As a future perspective, we plan to perform more rigorous scale settings, by exploiting state-of-art measurements of lattice potential. It is also mandatory to investigate the chiral limit and the thermodynamic limit at large N f . This, together with a more extended set of flavour numbers, will allow a quantitative analysis of the critical behaviour in the vicinity of the conformal IR fixed point. We expect that our thermodynamic lattice study for the large N f non-Abelian gauge theory plays an important role as a new connection between the lattice and the Gauge/Gravity duality [64,7].

Acknowledgements
We enjoyed discussing these topics with George Fleming, Edward Shuryak, Anna Hasenfratz, Philippe de Forcrand, and Marc Wagner. We also wish to thank Carleton DeTar and Urs Heller for sharing their notes on the normalization of chiral observables in the MILC code. This project is one step in our study of the phase diagram of strong interactions, and we warmly thank Elisabetta Pallante, Albert Deuzeman and Tiago Nunes da Silva for many interesting conversations and a pleasant collaboration. We also thank them for granting access to some of u 0 data at N f = 8 used in this study, and in particular Tiago Nunes da Silva for related communications. This work was in part based on the MILC Collaboration's public lattice gauge theory code [51]. The numerical calculations were carried out on the IBM-SP6 and BG/P at CINECA, Italian-Grid-Infrastructures in Italy, and the Hitachi SR-16000 at YITP, Kyoto University in Japan.

Appendix A. Summary table for simulation results
In this appendix, we summarise our results and the parameters used in the simulations and analyses. We have used an improved version of the staggered action, the Asqtad action [52], with a one-loop Symanzik [53,54] and tadpole [55] improved gauge action. All our simulations used the same bare fermion mass ma = 0.02.
In Tables A.4 -A.14 we will quote the lattice bare couplings β L = 10/g 2 L , tadpole factors u 0 , step lengths for a single trajectory ∆τ = δτ × 20, the number of total trajectories n traj , the number of thermalized trajectories n ave which have been used to evaluate ensemble averages, bin-sizes for jackknife error estimates s bin , ensemble averages for the chiral condensate (2a 3 ψ ψ with the definition Eq. (5)) the Polyakov loop (Re L or |L| with the definition Eq. (10)), and the chiral susceptibility (χ σ , defined in Eq. (7)) and/or the susceptibility (χ |L| ) associated with the absolute value of Polyakov loop.
Here are several technical comments for the symbols: In the second column in each table, some values for the tadpole factors u 0 have been estimated by using the fit with the ansatz u 0 = 1.0 − A/(1.0 + B β 2 L ), and are shown inside of the parentheses. Some of the total trajectories in the third column have a symbol +, for which the configurations obtained in the other simulations has been utilised as the inputs. The Monte Carlo step ∆τ in the fourth column has been adjusted to give about 75 − 80 percent Metropolis acceptances in the pre-thermals domain. In the fifth column, some entries have a symbol * , which indicates that additional simulations would be preferable to confirm the thermalization though the trajectories seem to have reached a stable domain. In the seventh column, a pre-factor 2 appears in front of a 3 ψ ψ following the difference between the the MILC code convention and the standard definition Eq. (5).  (7) 2.51 ± 0.69