Conformality in many-flavour lattice QCD at strong coupling

It is widely believed that chiral symmetry is spontaneously broken at zero temperature in the strong coupling limit of staggered fermions, for any number of colors and flavors. Using Monte Carlo simulations, we show that this conventional wisdom, based on a mean-field analysis, is wrong. For sufficiently many fundamental flavors, chiral symmetry is restored via a bulk, first-order transition. This chirally symmetric phase appears to be analytically connected with the expected conformal window of many-flavor continuum QCD. We perform simulations in the chirally symmetric phase at zero quark mass for various system sizes L, and measure the torelon mass, the Dirac spectrum and the hadron spectrum. All masses go to zero with 1/L. L is hence the only infrared length scale. Thus, the strong-coupling chirally restored phase appears as a convenient laboratory to study IR-conformality. Finally, we present a conjecture for the phase diagram of lattice QCD as a function of the bare coupling and the number of quark flavors.


Introduction
The possibility that the Higgs boson could be a composite bound-state in a high-energy Technicolor theory [1] has generated considerable interest, especially in the lattice community. In particular, the requirement that the Technicolor theory be "walking" [2], in order to accommodate stringent bounds on flavor-changing neutral currents, has been the driving motivation behind several large-scale computer simulation efforts to determine the possible combinations of gauge groups and fermion contents leading to a conformal window.
To determine via lattice Monte Carlo simulations whether a given theory is inside the conformal window is particularly challenging, because it involves a triple difficulty: in order to identify (or not) an infrared fixed point (IRFP) which is the signature of a theory inside the conformal window, one must probe the extreme infrared properties of the theory, while at the same time taking the continuum limit of the lattice discretization, and controlling the limit when the quarks become massless. This compounded difficulty may explain why, in spite of considerable efforts, there is no consensus yet on the minimum number N f * of quark flavors needed for QCD to be inside the conformal window [3]. A numerical demonstration of walking has been provided only recently, in a toy model, the 2-d O(3) model at vacuum angle θ ≈ π [4].
Here, we relax the demand that results should be obtained in the continuum limit. On a coarse lattice, long-distance properties can be studied more economically. While such properties may differ from those of the corresponding continuum theory, it may still be instructive to consider the possible existence of an IRFP for a discretized lattice theory. The phase diagram of SU(N c ) gauge theory with N f fundamental fermions, as a function of N f and the bare gauge coupling, has been predicted in the celebrated Ref. [5], which serves as a guide to understand the results of Monte Carlo simulations performed at finite bare coupling. It is important to confront these predictions with uncontroversial numerical evidence. Therefore, we start our investigation by considering the strong coupling limit, where the lattice is maximally coarse.
Note that we consider standard staggered fermions, and (away from the strong-coupling limit) the standard plaquette action. Other discretizations could lead to different results, since only the continuum limit is universal.
The conventional wisdom for strong coupling QCD with staggered fermions is that chiral symmetry remains always broken at zero temperature, regardless of the number of colors and flavors. This belief is based on mean-field analyses performed in some of the earliest papers on lattice QCD. In particular, it was shown in [6] that at leading order in a 1/d expansion, the chiral condensate has a value independent of the number of colors N c and of the number of staggered fieldsN f = N f /4, where N f would be the corresponding number of degenerate fermion flavors in the weak-coupling limit, but depends only on the number d of spatial dimensions: Chiral symmetry may be restored by increasing the temperature T . Following the approach of [7], where explicit results are provided for a few small values of N c andN f , we calculated the chiral restoration temperature aT c and found that it is indeed non-zero for allN f , and independent of N c to leading order in 1/N f : Hence chiral symmetry will never be restored at zero temperature, according to the meanfield analysis. Since mean-field theory is expected to work well when the number of d.o.f. per site is large (e.g. providing exact results in the Gross-Neveu model for N f → ∞), there was no reason to doubt the validity of this finding. Besides, it is in accord with the intuition that the gauge field is maximally disordered in the strong-coupling limit, and that this disorder will drive chiral symmetry breaking. On the other hand, one naively should expect the above disorder to be modified by dynamical fermions, which have an ordering effect. Indeed, the loop expansion of the determinant shows that the fermionic effective action induced by dynamical fermions, S eff = − log det( / D + m q ), starts with a positive plaquette coupling ∆β, proportional to 1/m 4 q for heavy quarks, which has been studied numerically in [8]. Clearly, for N f flavors the effective action is proportional to N f , and ∆β grows proportionally. This plaquette term suppresses fluctuations in the gauge field, which suggests that chiral symmetry restoration might take place for sufficiently large N f .   Figure 2. The chiral condensate at weaker coupling, β = 5, in the (N f , am q ) plane, for 4 4 (left) and 6 4 (right) lattices.

Monte Carlo results
The only way to resolve this puzzle is to perform Monte Carlo simulations in the strong coupling limit of staggered fermions, to detect a possible chiral symmetry restoration for sufficiently largeN f . These simulations are straightforward, using the standard Hybrid Monte Carlo algorithm. As expected, the effect of increasingN f on the chiral condensate is to reduce its magnitude. But it came as a surprise to find that the chiral condensate vanishes via a strong first-order transition atN f c 13 staggered fields in the chiral limit (i.e. N f c 52 continuum fermion flavors). In the broken phase, the chiral condensate remains almost constant. It vanishes in the chiral limit due to finite-size effects only. In contrast, in the chirally restored phase the condensate is caused by explicit symmetry breaking and is proportional to the quark mass. This is illustrated Fig. 1, where the condensate is shown as a function ofN f and bare quark mass (am q ). Moreover, this N f -driven transition turns out to be a bulk, zero-temperature transition, which can be seen by the fact that finite-size effects on the phase boundary are small when comparing two different system sizes 4 4 and 6 4 , as shown Fig. 1 (left and right).
One also observes in these figures that the critical number of flavors increases with the quark mass. This is easy to understand: heavier quarks have a weaker ordering effect, so that the induced plaquette coupling ∆β decreases if one keepsN f fixed. It takes more flavors to keep the system chirally symmetric. Hence,N f c increases, and for heavy quarks One can now go back to the mean-field treatment and trace the origin of its failure. Two kinds of terms at least are neglected: (i) multiple meson hopping along a given link, and (ii) baryon loops. These terms amount to corrections O(N f /N c ) and O(N f /d 2 ), respectively, whereN f = N f /4 is the number of staggered fields, and is normally set to 1 in the mean-field treatment. Here, we considerN f 13, and the previously neglected corrections become dominant. The conventional wisdom that chiral symmetry is always broken at T = 0 in the strong-coupling limit comes from mistakenly applying the lowest-order mean-field approximation in a regime where it is invalid.
Having established anN f -driven phase transition in the strong-coupling limit, we may consider its impact on the lattice theory at non-zero lattice gauge coupling β as well. Since the transition is strongly first order, it has to persist for some range in β at least. Hence we have compared the strong coupling phase diagram with the phase diagram at weaker coupling β = 5, illustrated Fig. 2. We find a similar qualitative behavior, but withN f c drastically reduced to O(2). Finite-size effects are more pronounced, but the transition still seems to be a first-order bulk transition.
In fact, we find a smooth variation of theN f -driven transition with β at a given small quark mass am q = 0.025, as shown Fig. 3. The transition extends to weak coupling, at least to β = 5, and remains strongly first-order. Thus, it is plausible that this transition, which separates a chirally broken (small N f ) and a chirally symmetric (large N f ) phase, persists all the way to the β → ∞ continuum limit, where it is to be identified with the transition at N f = N f * between the chirally broken and the IR-conformal, chirally symmetric phase.
In other words, our chirally restored phase may be analytically connected to the conformal window in the continuum limit, because we do not observe any additional non-analyticity as β is increased. This possibility motivates our study of the properties of the strong-coupling chirally symmetric phase, looking for tests of IR-conformality.

Looking for conformality in the chirally symmetric phase
It is natural to ask whether the chirally restored phase is connected to the conformal window, i.e. whether the chirally restored phase at strong coupling is also IR-conformal. And if this is indeed the case, the next obvious question is whether this IR-conformal phase is trivial, ie. whether the IR fixed point coupling is zero or not. In this section we present measurements of gluonic and fermionic observables chosen to address these questions: the torelon mass from which we define a running coupling, the Dirac eigenvalue spectrum, and the hadron spectrum. Our results support the following conclusion: the strong-coupling chirally symmetric phase is indeed IR-conformal, and it is non-trivial.
The simulations performed here are all in the chirally symmetric phase at zero plaquette coupling, withN f = 14 and 24 staggered fields, which would correspond, in the weakcoupling limit, to N f = 56 and 96 continuum flavors, and with lattices of size 4 3 × 16,  The quark mass is set exactly to zero unless specified otherwise. We will see below that the Dirac operator has a spectral gap in the symmetric phase, which makes a study of the massless theory within reach of modest computer resources. Moreover, having one infrared scale, the system size L, rather than two scales (L and 1/m q ) is of great advantage when analyzing the results. Let us mention the average plaquette values which we measure: ≈ 0.35 and ≈ 0.52 for N f = 14 and 24, respectively (normalized to 1 for the free field). So we are very far from a plaquette value of 0, corresponding to maximally disordered gauge fields and achieved for N f = 0: the ordering effect of the dynamical fermions plays a dominant role in our case, and the vanishing of the plaquette coupling is not associated with special properties. The "torelon" is a gluonic excitation which is topologically non-trivial: it is excited by any Wilson loop which wraps around the spatial boundary in one direction, for instance, as illustrated Fig. 4, where i = 1, 2, 3 is one of the spatial directions andî the unit vector in that direction. 1 We extract the mass of this excitation from the exponential decay of the correlator To suppress excited states, we smear the links within each time-slice before constructing T i . This observable has been used for a long time to extract the string tension σ in Yang-Mills theories [9]: it can be viewed as a loop of gluonic string, whose energy m T (L) grows with its length as σL. So our initial, naive expectation was to measure a mass m T (L) growing with L, until perhaps the string would break due to fermion-pair creation. This is not at all what we observed. The dimensionless quantity which we measure, am T (L), decreases on larger lattices corresponding to a larger ratio L/a. Clearly, our theory is not confining. Moreover, as shown Fig. 5, the combination Lm T (L) is approximately constant as L is increased. So the torelon mass varies as 1/L (actually, for small L it initially decreases even faster with L as seen in the figure). Thus, there is no intrinsic mass scale which appears in this channel: the torelon mass is set by the system size L. This remarkable result is our first evidence that our theory is IR-conformal.
We actually combined the T i (t), i ∈ {1, 2, 3} into two representations of the cubic discrete rotation group: the A + 1 (corresponding to a 0 + representation of O (3)) and the E + (corresponding to a 2 + representation of O (3)). The mass of the E + seems to be slightly smaller, as observed in small-volume analytic Yang-Mills calculations [10]. Now, by relabelling the spatial direction i as the imaginary time direction, one realizes that we are measuring the correlation of two time-like Polyakov loops, whose decay rate is governed by the Debye mass, given perturbatively at lowest order by This expression allows us to define a running coupling g(L) via and we see that, in this scheme, our running coupling seems to go to a non-zero constant as L increases (although one cannot exclude, of course, that it slowly goes to zero). Therefore, we have numerical evidence supporting the view that our strong-coupling, chirally symmetric theory is IR-conformal and non-trivial.
Interestingly, the extracted value of g(L) approaches ∼ 0.95 and ∼ 0.80 forN f = 14 and 24 respectively. So the IR fixed-point coupling value decreases asN f increases. This is what one would expect: asN f keeps increasing, the ordering effect of the fermions increases, and all Wilson loops are driven towards 1, their free field value. At the same time, any definition of a running coupling will approach zero. The theory becomes trivial forN f → ∞, even in the strong-coupling limit. We will come back to this point in Sec. IV.

Characterizing the chirally restored phase: (II) Dirac Spectrum
We now turn to fermionic properties, starting with the spectrum of the Dirac operator. We have analyzed the Dirac eigenvalue spectrum of the configurations in our Monte Carlo ensembles, using a Lanczos algorithm to obtain an approximation of the whole spectral density and an Arnoldi method to extract the smallest eigenvalues to high accuracy. The observable shown in Fig. 7 is the integrated eigenvalue density, defined as: This function of λ counts the fraction of eigenvalues smaller than λ. Its derivative is simply the eigenvalue density ρ(λ). We first compare this observable on quenched configurations (N f = 0) and in the chirally symmetric phase (N f = 14). In Fig. 7 (left), 10 curves for   10 configurations are superimposed: variations in the spectrum are very small. We observe that N f = 0 andN f = 14 spectra are similar in the ultraviolet, but differ in the infrared, as illustrated in the inset. The N f = 0 curve starts linearly from the origin, reflecting an eigenvalue density approximately constant near λ = 0. On the contrary, the integrated eigenvalue density forN f = 14 shows a spectral gap for small eigenvalues, which is of course consistent with chiral symmetry restoration according to the Banks-Casher relation, since The crucial question is on which scale does this spectral gap depend. To answer this question, in Fig. 7 (right) we compare the integrated eigenvalue density atN f = 14 for various lattice volumes, L = 4, 6, 8, 10 and 12. As evidenced in the inset, we find that the spectral gap scales ∝ 1/L to a good approximation, which is a strong indication that our theory is IR-conformal 2 : There does not seem to be any length scale in the chirally restored phase other than the box size L. Actually, small deviations from 1/L scaling allow us, in principle, to extract the anomalous mass dimension γ * . We make such an attempt in Fig. 8, where the gap has been multiplied by L already: deviations from a constant are indicative of anomalous dimension, provided other corrections O((a/L) 2 ) are negligible. The effect of a finite system size L on the Dirac spectrum has not been analyzed yet. We have simply considered that the infrared conformal symmetry is explicitly broken by the infrared scale 1/L, which is the analogue of an explicit breaking by a quark mass m q . Consequently, we expect the mass gap to behave as (1/L) 1/(1+γ * ) . A crude, 2-parameter fit based on our 3 largest volumes gives γ * ∼ 0.26 and 0.38 forN f = 14 and 24, respectively. Simulations on larger volumes should be performed to bring under control the systematic error in these estimates. The true, infinite volume value of γ * seems to be approached from below.

Characterizing the chirally restored phase: (III) Hadron Masses
Finally, we turn to hadron masses measured on ourN f = 14 andN f = 24 Monte Carlo ensembles. Even though the quark mass is zero in these ensembles, we observe non-zero hadron masses. As expected, parity partners are degenerate since chiral symmetry is restored. Now, if our theory is IR-conformal, the masses which we measure are exclusively due to finite-size effects: all masses should go to zero as the lattice size L is increased. This is what we observe, as shown in Fig. 9 (left): the "pion" and "rho" masses both decrease by a factor ∼ 2 as the lattice size is increased from L = 4 to 12. Notice however that the smallest mass is still > 0.6, which is not very light.
Furthermore, one generally expects that the approach to zero should be the same for all hadrons, so that mass ratios should remain constant as L → ∞. Note that there may be exceptions to this "rule": in the 2d O(3) sigma model near θ = π, the mass of the O(3) singlet state approaches zero faster than that of the O(3) triplet as θ approaches π [11]. Here, our Fig. 9 (left) would show all data points aligned on "rays" going through the origin if mass ratios were constant. One can see deviations from this behavior, which perhaps are caused by the not-too-light masses which we measure. Another possible cause is technical: as L is increased, the groundstate masses in each channel decrease, but so do also the mass differences between groundstate and excited states. It becomes more difficult to extract the groundstate mass, and our lattices are likely too short to properly control this source of systematic error. Nevertheless, we show in Fig. 9 (right) the mass of the "pion", in which we have the most confidence, as a function of 1/L. Since 1/L breaks the conformal symmetry much like a quark mass m q would, we expect that the pion mass should scale the same way, namely as (1/L) 1/(1+γ * ) , if L is large enough. A 2-parameter fit to our three largest system sizes gives γ * ∼ 1.0 and 0.4 forN f = 14 and 24, respectively. As with the estimates of γ * from the Dirac spectrum, simulations on larger volumes are needed to bring the systematic error under control.

The conjectured phase diagram
Using both gluonic and fermionic observables, we have presented evidence that the chirally restored phase at strong coupling is IR-conformal and non-trivial, and speculated about a connection to the conformal window in the continuum. We now want to propose a phase diagram sketched in Fig. 10 (left), as a function of the plaquette coupling β = 6/g 2 0 and of the number of would-be fundamental flavors N f in the weak-coupling limit β = ∞. That is, we simply convert the numberN f of staggered fields to N f = 4N f . Moreover, we promote N f to a real, continuous parameter: while N f must be integer for a welldefined continuum theory, one may let it take any value in the statistical model defined by the lattice partition function. Our conjectured phase diagram can be compared with, e.g., those of Ref. [5,12]: one can see substantial differences. The essential feature of our phase diagram is that the β = 0 IR-conformal phase is analytically connected with the weak-coupling, continuum IR-conformal phase. This is the simplest scenario, supported by our exploratory scan in β shown Fig. 3. A single phase transition line N f c (β) separates the region of broken chiral symmetry at small N f from the chirally symmetric region at large N f . The transition is first-order, at least for some range of β starting from zero. Moreover, the numberN f c of staggered fields which bring enough order to restore chiral symmetry at (1), is remarkably close to the expected number N f * of continuum quark fields which achieve the same effect 3 . This may be more than a numerical accident. At strong coupling, taste symmetry breaking is maximal, andN f = 13 staggered fields can be viewed asN f massless fields, plus 3N f fields with mass O(a −1 ).
Only the former have a significant effect toward chiral symmetry restoration.
In the chirally symmetric phase, we see no evidence for a dynamically generated mass scale of any sort. Then, based on our results for the running coupling Fig. 6, we conjecture that for any finite β and N f > N f c (β), large-N f lattice QCD is IR-conformal, with a nontrivial fixed-point coupling g * > 0. This value changes continuously with (β, N f ), reaching the value zero for β = ∞, N f > 33/2 and for N f = ∞ ∀β, as indicated Fig. 10 (left) by a thick dotted line. g * grows as one moves away from this dotted boundary towards the phase transition line. An alternative scenario would be that the running coupling in Fig. 6 slowly approaches zero, and g * = 0 for N f > 33/2 for any β in the chirally symmetric phase. This is sketched in Fig. 10 (right). If the basin of attraction of the weak coupling trivial fixed point would extend all the way to the strong coupling limit, one should observe for the running coupling g 2 (L) ∼ 1/ log(L/L 0 ), with L 0 ∼ O(a). Whether or not this happens depends on the marginal operators induced by our lattice action. Our numerical results Fig. 6 for the running coupling are indeed consistent with this possibility. But we should also observe in this case an anomalous mass dimension γ * = 0, which is not favored by our other measurements. More careful, large-scale simulations are necessary to settle this issue.
Finally, one may consider the line g * = g 0 , with g 0 = (2N c /β) 1/2 , where the IR fixedpoint coupling has the same value as the bare coupling, so that the coupling does not run as a function of the renormalization scale. This line starts at the point (β = ∞, N f = 33/2) where g * = g 0 = 0. Its precise location depends on the chosen renormalization scheme. It is not associated with any kind of singularity of the free energy. There is no phase transition along this line: simply, on the left (resp. on the right) of that line, the coupling increases (resp. decreases) from g * as one reduces the distance scale. Since there is a lower distance cutoff a no divergence is observed as one crosses this [scheme-dependent] line.
We have determined the phase diagram in the strong-coupling region only. Studying the continuum limit is of course much more difficult, due to the large lattices that have to be used in order to control the finite size effects and the difficult control of lattice artifacts. We would like to suggest that the strong-coupling limit may represent an advantageous "poor man's laboratory" for the study of 4d IR-conformal gauge theories. In particular, as weak coupling: strong coupling: Figure 11. The ordering of scales at weak coupling (left) and strong coupling (right), showing that the range of conformal invariance is larger in the latter case.
illustrated Fig. 11, the range of scales over which conformal invariance applies for a given computing effort is greatly reduced at weak coupling: there, for a given lattice size N 4 , the scales are ordered as follows: where Λ is the scale generated by the asymptotically free gauge dynamics and Yang-Mills perturbation theory applies at distances 1/Λ. In contrast, at strong coupling, where the lattice becomes maximally coarse, there is no small distance where Yang-Mills perturbation theory applies, and the hierarchy is: Hence the dynamical range of conformal invariance, characterized by the product LΛ, is maximized at β = 0 for a given lattice size N = L/a.

Conclusion
We have shown that for β = 0, contrary to common wisdom, there exists a strongly firstorder, N f -driven bulk transition to a chirally symmetric phase. In the chiral limit, the transition occurs forN f c = 13(1) staggered fields, i.e. N f c = 52(4) continuum flavors.
This finding is in contrast to the mean-field prediction, whose failure can be traced back to approximations relying on N f being small. Clearly, the conventional, automatic association of the strong-coupling limit with confinement and chiral symmetry breaking is too naive. Furthermore, the chirally restored phase extends to weak coupling. We have also shown numerical evidence that the β = 0 chirally restored phase of "large-N f QCD" is IR-conformal, with a non-trivial, N f -dependent value of the IR fixedpoint coupling. We conclude that the strong-coupling limit is the laboratory of choice to study a 4d IR-conformal gauge theory. Simulations at large N f and zero quark mass can be performed without too much computational effort since a gap appears in the Dirac spectrum. As N f increases, the spectral gap increases, the average plaquette approaches 1, and the fixed-point coupling approaches 0. Setting the quark mass to zero eliminates one IR scale, leaving the system size L as the only remaining one. This greatly simplifies the analysis of simulation results.
Since we have not observed any evidence for an additional T = 0 phase transition as β is increased, we speculate that the strong coupling chirally symmetric, IR-conformal phase is analytically connected with the continuum IR-conformal phase.
One may ask how robust these statements are with respect to the particular discretization of the Dirac operator and the gauge action. While the quantitative details of the phase transition N f c (β) will surely change, we think that the qualitative features will remain.
Chiral symmetry breaking at strong coupling, for small N f , is a general consequence of the disorder in the gauge field. The ordering effect of fermions also is generic. So we do expect a bulk transition, generically of first-order, as a function of N f in the strong-coupling limit. Actually, such a transition was observed for Wilson fermions in Ref. [13,14]. At intermediate coupling, additional transitions may appear depending on the lattice action.
Interestingly, a first-order transition to a chirally broken phase as β is reduced has been observed many times, for various lattice actions [15][16][17][18][19]. These transitions were observed for some fixed value of N f . Here, we simply put all these earlier observations together. It is interesting that this phase transition is consistently of first-order. If the first-order nature persists all the way to the continuum limit, then walking dynamics will not be observed, and the transition to the conformal window will be characterized by "jumping dynamics", as proposed by Sannino [20]. One may also wonder what happens to the (β, N f ) phase diagram as the gauge group or the fermion content is changed. For 4d compact U (1), the change would not be dramatic, because the strong-coupling behavior is much the same as for SU (3): our first-order transition line would end at (β ≈ 1.01, N f = 0) on the horizontal axis rather than on the vertical axis, and the region of triviality would cover the whole chirally symmetric phase, except for N f = 0. For SU (2) or SU (3) with adjoint fermions, the change would be more significant: in the strong-coupling limit, increasing N f would order the plaquette in the adjoint representation, not in the fundamental. Center monopoles would likely condense [21], and might delay or prevent the restoration of chiral symmetry.
Finally, there are many directions in which to extend this exploratory study. To buttress the claim that the chirally symmetric phase is IR conformal, more observables, like the static potential and the Fredenhagen-Marcu order parameter, should be studied. Also, and to make contact with other numerical studies, a mass deformation could be introduced. As a first step in this direction, we show in Fig. 12 the quark mass dependence of the chiral condensate. This figure shows all the technical difficulties associated with extracting the anomalous dimension γ * : Heavier fermions have less ordering effect, which triggers a phase transition back into the chirally broken phase for some critical fermion mass. Finite-size effects associated with that transition should not be included in the determination of γ * . Moreover, the non-anomalous contributions ψ ψ = c 1 m q + c 2 m 3 q easily overwhelm the anomalous m or including an m 3 q term [22], favor negative values for γ * , which crucially depend on the fitting range and the included analytic contributions. We believe that extracting γ * from the Dirac spectrum provides better control of the systematics, as emphasized in [23]. In any case, larger system sizes are required before reliable estimates of γ * can be obtained. This is beyond the scope of this work. Here, we have argued that such reliable estimates can be obtained in principle from the L-dependence of the Dirac spectrum and of the meson spectrum, both measured at zero quark mass.