Nonperturbative beta function of twelve-flavor SU(3) gauge theory

We study the discrete beta function of SU(3) gauge theory with Nf=12 massless fermions in the fundamental representation. Using an nHYP-smeared staggered lattice action and an improved gradient flow running coupling $\tilde g_c^2(L)$ we determine the continuum-extrapolated discrete beta function up to $g_c^2 \approx 8.2$. We observe an IR fixed point at $g_{\star}^2 = 7.3\left(_{-2}^{+8}\right)$ in the $c = \sqrt{8t} / L = 0.25$ scheme, and $g_{\star}^2 = 7.3\left(_{-3}^{+6}\right)$ with c=0.3, combining statistical and systematic uncertainties in quadrature. The systematic effects we investigate include the stability of the $(a / L) \to 0$ extrapolations, the interpolation of $\tilde g_c^2(L)$ as a function of the bare coupling, the improvement of the gradient flow running coupling, and the discretization of the energy density. In an appendix we observe that the resulting systematic errors increase dramatically upon combining smaller $c \lesssim 0.2$ with smaller $L \leq 12$, leading to an IR fixed point at $g_{\star}^2 = 5.9(1.9)$ in the c=0.2 scheme, which resolves to $g_{\star}^2 = 6.9\left(_{-1}^{+6}\right)$ upon considering only $L \geq 16$. At the IR fixed point we measure the leading irrelevant critical exponent to be $\gamma_g^{\star} = 0.26(2)$, comparable to perturbative estimates.


Introduction
SU(3) gauge theory with N f = 12 flavors of massless fermions in the fundamental representation has been considered by many independent lattice studies in recent years. This effort is motivated by the expectation that the 12-flavor system exhibits conformal or nearconformal dynamics qualitatively different than QCD. That is, N f = 12 is likely either within or close to the lower boundary of the SU(3) conformal window N (c) f ≤ N f < 16.5, where the theory flows to a chirally symmetric conformal fixed point in the infrared (IRFP) [1,2]. Should the system undergo spontaneous chiral symmetry breaking (i.e., 12 < N (c) f ), then it provides an example of a strongly coupled theory in which lattice calculations have observed a light 0 ++ scalar [3,4]. In this case investigations of N f = 12 are relevant to explore possible strongly coupled new physics beyond the standard model (BSM), in which such a light composite scalar could be consistent with the observed SM-like Higgs boson [5,6]. Alternatively, if the 12-flavor system is within the conformal window, as our results indicate, it provides a useful testbed in which to develop and apply nonperturbative methods to investigate IR-conformal systems. Even in this case there can be connections to BSM phenomenology, in models where the mass of some of the fermions is lifted to guarantee spontaneous chiral symmetry breaking. Lattice investigations of this situation have shown that this system follows hyperscaling, a highly non-QCD-like behavior, exhibiting natural large scale separation and UV dynamics dominated by the 12-flavor IRFP [7,8].
For example, step-scaling studies of the discrete β function directly search for an IRFP within a particular range of renormalized couplings. The exactly massless fermions typically employed by such studies make it more difficult for them to explore spontaneous chiral symmetry breaking, which finite-temperature or spectral techniques are better suited to investigate. If no IRFP is observed by step-scaling studies (as in recent work on N f = 8 [77,78]), then additional computations with am > 0 are needed to investigate chiral symmetry breaking in the considered range of couplings. Without identifying spontaneous chiral symmetry breaking in the am = 0 limit it remains possible for there to be an IRFP at some stronger coupling beyond the range in which the discrete β function was explored. As spontaneous chiral symmetry breaking is an inherently non-perturbative phenomenon we wish to probe it using lattice calculations rather than relying on imprecise estimates of the critical coupling strength g 2 MS ∼ 10 [16,17]. In the case of N f = 12, the pioneering step-scaling study of refs. [23,24] identified an IRFP at g 2 SF ≈ 5 in the Schrödinger functional scheme (with purely statistical uncertainties 10%). Subsequent investigations [25][26][27][28][29][30][31][32][33][34] have attempted to improve upon this result by considering larger lattice volumes, different schemes for the running coupling, and improved lattice actions with smaller discretization artifacts. Two recent large-scale projects are of particular note. Ref. [33] explores the discrete β function up to g 2 c 6 in the c = √ 8t/L = 1 A recent five-loop β function computation [11] appears to change this trend, although the subsequent refs. [12][13][14][15] argue that all systems with 9 ≤ N f ≤ 16 exhibit IRFPs at the five-loop level. We address this development in section 6. 2 At the perturbative g 2 = 0 fixed point staggered lattice fermions are equivalent to continuum fermions.
At a non-trivial IRFP this is not necessarily the case; instead, the different chiral symmetry properties of different lattice fermion formulations could correspond to different fixed points. Such behavior has been studied in three-dimensional spin systems [75,76].
0.45 and 0.5 gradient flow schemes with color-twisted boundary conditions (BCs) and an unimproved lattice action. Although the resulting step-scaling function approaches zero it does not vanish in the accessible range of couplings, and a bulk transition into a lattice phase obstructs progress to larger g 2 c . Ref. [34] employs very large lattice volumes and an improved action to explore the very narrow region 6 g 2 c 6.4 in the c = 0.2 gradient flow scheme, also obtaining a non-zero discrete β function. 3 As we show in figure 5, both of these investigations are consistent with our full c = 0.25 and 0.3 results that predict an IRFP at g 2 = 7.3 +5 −3 (despite the slightly different renormalization schemes considered). In addition to the step-scaling studies summarized above, most other N f = 12 investigations offer further evidence supporting the existence of a conformal, chirally symmetric IR fixed point. Investigations of the phase diagram both at zero and finite temperature have observed a first-order bulk phase transition that extends from the am = 0 chiral limit to non-zero mass [35,38,[41][42][43][44][45][46]. At finite temperature T = 1/(aN t ), where a is the lattice spacing and N t is the temporal extent of the lattice, the chiral transition lines run into the bulk phase at non-zero mass [42,45]. This is a necessary condition for IRconformality, where the finite-temperature transitions in the chiral limit must accumulate at a finite coupling as N t → ∞, and remain separated from the weak-coupling conformal phase. No lattice investigations of the 12-flavor phase diagram have been able to identify spontaneous chiral symmetry breaking in the form of chiral transitions that remain in the weakly coupled phase upon extrapolation to the chiral limit.
Spectral studies offer another means to explore the IR dynamics, by confronting nonzero-mass lattice data with expectations based on either chiral perturbation theory or conformal finite-size scaling. While refs. [53,56,58] observe consistency with conformal hyperscaling for N f = 12, ref. [49] reported a very low level of confidence in conformality. However, subsequent re-analyses of the data published by ref. [49] suggest that this conclusion is sensitive to the details of the analyses [50,51,56]. In particular, by taking into account corrections to scaling arising from the nearly marginal (i.e., slowly running) nature of the gauge coupling, in ref. [56] we were able to carry out consistent finite-size scaling analyses that included both our own spectrum data as well as those published by refs. [49,53].
Our finite-size scaling study predicted the scheme-independent mass anomalous dimension γ m = 0.235(15) at the 12-flavor IR fixed point. A similar result γ m = 0.235 (46) was reported by ref. [58]. 4 In addition, our studies of the massless Dirac operator eigenmodes independently predict γ m ≈ 0.25 [61,62]. These results are quite close to the four-loop perturbative prediction γ m = 0.253 in the MS scheme [9], and the new five-loop result γ m = 0.255 [13], though a recent scheme-independent series expansion [81] obtains a larger 3 This particular range of g 2 c was chosen based on some results in our earlier publication [31], which identified a 12-flavor IRFP at g 2 = 6.2 (2). In appendix B we discuss the reasons for the discrepancy between our full results and that previous work. 4 Finite-size scaling analyses without corrections to scaling typically obtained larger values that often varied non-universally depending on the observables analyzed: 0.2 γm 0.4 [54], γ m = 0.403(13) [50], γ m 0.35 [51] and γ m = 0.4-0.5 [53]. These results are all consistent with an upper bound γ m ≤ 1.29 from the conformal bootstrap program [79], though not with the perturbative γ m ≈ 1.3-1.5 reported by ref. [80]. γ m = 0.400 (5) [14,15]. This small, potentially perturbative mass anomalous dimension, in combination with the assumption that γ m 1 around the lower edge of the conformal window, may suggest that N f = 12 is quite deep within the conformal regime.
Despite the many high-quality, large-scale investigations of the 12-flavor system summarized above, there is still progress to be made in resolving its IR properties. In this work we report our final results on the step-scaling calculation of the discrete β function for N f = 12. These results supersede the partial analysis included in ref. [31], and predict a conformal IR fixed point at g 2 = 7.3 +3 −2 in the gradient flow scheme with c = 0.25. We also investigate the slope of the step-scaling function at the IRFP, both directly and via finite-size scaling as in refs. [24,33]. This slope is related to the leading irrelevant critical exponent γ g , for which we find γ g = 0.26 (2), consistent with the four-loop perturbative prediction γ g = 0.282.
Compared to ref. [31] we have accumulated significantly more data, in particular generating several new lattice ensembles at relatively strong couplings β F 4 on each lattice volume up to 36 4 . This allows us to explore the discrete β function up to g 2 c 8.2, extending well past the IRFP. We now compare multiple discretizations of the energy density E(t) in the gradient flow renormalized coupling, obtaining consistent results. Finally, we add two new lattice volumes, 20 4 and 30 4 , that allow us to omit the 12 4 volume used in ref. [31]. As we show in appendix B, analyses that include 12 4 volumes in the c = 0.2 gradient flow scheme suffer from particularly large systematic uncertainties that were not comprehensively considered in ref. [31].
Although our 12-flavor results are qualitatively different than those we previously obtained for the 8-flavor discrete β function [77], much of our analysis follows the same procedure as that work, and the next three sections are organized in the same way. We begin by reviewing gradient flow step scaling in the next section, including the improvement of the gradient flow running coupling. In section 3 we describe our numerical setup and lattice ensembles. We use an nHYP-smeared staggered fermion lattice action [82,83], with both fundamental and adjoint plaquette terms in the gauge action [38,42,45,65]. We employ this same action in our 12-flavor finite-temperature [38,42,45], spectral [56] and eigenmode [61,62] studies summarized above, which can therefore be consistently compared. On each of eight L 4 volumes with 12 ≤ L ≤ 36 we generate between 14-35 ensembles at different bare couplings in the range 3 ≤ β F ≤ 9.
Our step-scaling analyses and results are presented in section 4, including discussion of systematic uncertainties from the stability of the (a/L) → 0 extrapolations, the interpolation of g 2 c as a function of the bare coupling, and the improvement of the gradient flow running coupling. We compare the clover and plaquette discretizations of E(t) as another consistency check, obtaining agreement in all cases we consider. Finally, we also confirm the consistency of our results with those recently reported by refs. [33,34]. In section 5 we investigate the leading irrelevant critical exponent from the slope of the step-scaling function at the IRFP, observing γ g = 0.26 (2), comparable to perturbative estimates. We check this result by carrying out a finite-size scaling analysis. We conclude in section 6 with some brief discussion of how our new results affect the broader context of 12-flavor lattice investigations summarized above, and highlight a few directions that merit further study in the future.
We include three appendices collecting some supplemental checks of our results. In appendix A we briefly consider the discrete β functions resulting from two scale changes s = 2 and s = 4/3 different from the s = 3/2 considered in the body of the paper. In contrast to s = 3/2, for both of s = 2 and 4/3 we are forced to include small-volume 12 4 lattice ensembles in our analyses. We obtain consistent results from all three scale changes, as summarized in table 2. However, as we show in appendix B, systematic uncertainties increase dramatically when combining smaller c 0.2 with smaller L ≤ 12. These systematic uncertainties were not comprehensively considered in the partial analysis we included in ref. [31], which reported g 2 = 6.2(2) with c = 0.2 and L ≥ 12, compared to the central value g 2 = 5.88 and upper bound g 2 ≤ 7.13 we now obtain with this choice of c and L min (table 2). Finally, appendix C provides a subset of our data.

Gradient flow step scaling and its improvement
We investigate a renormalized coupling defined through the gradient flow, which is a continuous transformation that smooths lattice gauge fields to systematically remove shortdistance lattice cutoff effects [84]. The demonstration that the gradient flow is mathematically well defined and invertible [85] inspired its use in a wide variety of applications (recently reviewed by ref. [86]). Here we consider the coupling [87] where the energy density E(t) is evaluated after 'flow time' t, corresponding to the energy scale µ = 1/ √ 8t. We will compare two lattice operators that can be used to define the energy density, first E(t) = − 1 2 ReTr [G µν (t)G µν (t)] with the symmetric clover-leaf definition of G µν , and second E(t) = 12(3 − (t)) where is the plaquette normalized to 3. The overall normalization N is set by matching g 2 GF (µ) with the continuum MS coupling at tree level. To carry out step-scaling analyses we tie the energy scale to the lattice volume L 4 by fixing the ratio c = √ 8t/L, as proposed by Refs. [88][89][90]. Each choice of c defines a different renormalization scheme, producing different results for the renormalized coupling g 2 c (L) and for the discrete β function in the continuum limit. When periodic BCs are used for the gauge fields, these β functions are only one-loop (and not two-loop) universal [88].
Extrapolating (a/L) → 0 is required to remove cutoff effects in the gradient flow renormalized couplings g 2 c . 5 These cutoff effects depend on the lattice action used to generate the configurations, on the gauge action used in the gradient flow transformation, and on the lattice operator used to define the energy density E(t). It is possible to systematically remove lattice artifacts by improving all three quantities simultaneously [91,92]. Here we take a simpler approach, using the Wilson plaquette action in the gradient flow transformation (i.e., the "Wilson flow") and combining two improvements that suffice to greatly reduce -and often essentially remove -cutoff effects. First, following ref. [93], we modify the definition of the renormalized coupling to perturbatively correct for cutoff effects, In this expression C(L, c) is a four-dimensional finite-volume sum in lattice perturbation theory, which depends on the action, flow and operator. We use the tree-level computation of C(L, c) from ref. [93], including a term that accounts for the zero-mode contributions allowed by the periodic BCs for the gauge fields. As we will see in figure 2, even this perturbatively improved gradient flow coupling can exhibit significant cutoff effects. While larger values of c 0.3 reduce these artifacts to some extent, this is accomplished only at the price of increased statistical uncertainties [90]. A better option, introduced in ref. [31], is to slightly shift the flow time at which the energy density is computed: with |τ 0 | t/a 2 . This t-shift τ 0 can be either positive or negative. Its effects vanish in the continuum limit where τ 0 a 2 → 0 so that g 2 GF (µ) = g 2 GF (µ). For O(a)-improved actions like those we use, choosing an optimal τ 0 value τ opt allows the removal of all O(a 2 ) corrections of the coupling g 2 GF (µ; a) defined in eq. 2.3. Although this optimal τ opt changes as a function of g 2 GF (µ), in this work we observe that τ opt depends only weakly on g 2 GF (µ), as in our previous studies of 4-, 8-and 12-flavor SU(3) systems [31,77]. Therefore we simply use a constant value of τ opt for all g 2 GF (µ), which suffices to remove most observable lattice artifacts throughout the ranges of couplings we explore.
Since we optimize τ 0 after applying the tree-level perturbative corrections discussed above, these two improvements do not interfere with each other. Nor do either of them require any additional computation, since the numerical integration through which we evaluate the gradient flow already provides all the data needed to shift t → t + τ 0 a 2 . Using the resulting g 2 c gradient flow running coupling, we will investigate the 12-flavor discrete β function corresponding to scale change s, . (2.4) We will also refer to this quantity as the step-scaling function σ s (u, L) with u ≡ g 2 c (L; a). To obtain our final results for the continuum discrete β function β s (g 2 c ) = lim (a/L)→0 β s ( g 2 c , L) we extrapolate (a/L) → 0.
We emphasize that different values of τ 0 should all produce the same β s (g 2 c ) in the continuum limit [31]. In appendix B we will show that this requirement is not satisfied for the lattice volumes we can access when c 0.2. In this case continuum extrapolations with different t-shifts disagree by statistically significant amounts, which likely contributes to the discrepancy between refs. [31] and [34]. In this work, when such sensitivity to the t-shift is present we will account for it as a source of systematic uncertainty, which was not done in ref. [31].
The different discretizations of E(t) should also produce the same β s (g 2 c ) in the continuum limit. We will separately analyze the plaquette and clover definitions of E(t), and find that they produce consistent results within uncertainties when c ≥ 0.25 and L ≥ 16. In appendix A we note that reducing L ≥ 12 requires increasing c ≥ 0.3 in order to maintain the good agreement between these two sets of results. When identifying the location of the IR fixed point, we will include the predictions of both discretizations in our determination of the total uncertainties on g 2 .

Numerical setup and lattice ensembles
Our numerical calculations use nHYP-smeared staggered fermions [82,83] with smearing parameters α = (0.5, 0.5, 0.4), and a gauge action including fundamental and adjoint plaquette terms with couplings related by β A /β F = −0.25 [38,42,45,65]. The fermions are exactly massless (am = 0), which freezes the topological charge at Q = 0. We impose anti-periodic BCs for the fermions in all four directions, while the gauge fields are periodic. Previous studies of this lattice action observed an " S 4 " lattice phase in which the singlesite shift symmetry (S 4 ) of the staggered action is spontaneously broken [38,42,45]. At am = 0 a first-order transition into the S 4 phase occurs at β We generate ensembles of gauge configurations with eight different L 4 lattice volumes with L = 12, 16, 18, 20, 24, 30, 32 and 36. Depending on L we study 14-35 values of the bare coupling in the range 3 ≤ β F ≤ 9. The 158 resulting ensembles are summarized in tables 3-10 in appendix C. These volumes allow us to consider three scale changes s = 2, 3/2 and 4/3, each with at least three pairs of volumes for continuum extrapolations as listed in table 1. In the body of the paper we focus on s = 3/2 where we can retain three points with L ≥ 16; we will see in the next section that the L = 12 ensembles exhibit potentially significant cutoff effects. Even so, we obtain comparable results for s = 2 and 4/3 analyses including L = 12 data, which are collected in appendix A.
We use the hybrid Monte Carlo (HMC) algorithm to generate configurations. Even at the strongest bare couplings we investigate we retain good HMC acceptance and reversibility in the am = 0 chiral limit with unit-length molecular dynamics trajectories and step sizes δτ ≈ 0.1 at the outer level of our standard multi-timescale Omelyan integrator. While the performance of the HMC algorithm is not a robust means to monitor the phase structure of the system, this behavior indicates that none of our ensembles exhibit chiral symmetry breaking. This conclusion is supported by our observation of a gap in the Dirac operator eigenvalue spectrum on many of these ensembles, including the strongest couplings β F ≥ 3 that we consider [61,62].
In figure 1 we show the gradient flow renormalized coupling g 2 c (L) measured on each ensemble for c = 0.25 and 0.3 (using the clover discretization of the energy density). These data use the optimal t-shift value τ opt = 0.08 that we discuss in the next section, and also include the tree-level perturbative correction factor C(L, c) in eq. 2.2. The perturbative corrections are fairly mild for our lattice action, Wilson flow, and clover or plaquette discretization of the energy density. The largest is C(12, 0.25) ≈ 1.12 for the plaquette discretization, with all others smaller than 6.2% effects. From these plots we can already see that the 12-flavor coupling runs very slowly, with little change in g 2 c (L) as the volume increases by a factor of three, especially for β F > 4.0. This feature of the system was mentioned in section 1, as the reason that finite-size scaling analyses need to account for the corresponding corrections to scaling.
The lines in figure 1 are interpolations using the rational function form in eq. 4.1. The plots in the bottom row zoom in on narrow regions of width ∆β F = 0.5 where the interpolations from different lattice volumes all cross each other. At the weak-coupling edge of these plots, β F = 5 (4.9) for c = 0.25 (0.3), the interpolated g 2 c monotonically increase with L from 12 to 36. At the strong-coupling edge, β F = 4.5 (4.4), the order has completely reversed and the interpolated g 2 c monotonically decrease as the lattice volume increases. Of course there are statistical uncertainties in the data that make the full analysis more complicated: To reduce clutter in these figures we don't display the uncertainties on the interpolations, within which most of the interpolations remain consistent with each other throughout much or all of this range.
The finite-volume crossings visible in these plots could be extrapolated to the infinitevolume limit to predict a 12-flavor IRFP, as in the c = 0.2 analysis of ref. [31]. 6 With c = 0.25 the crossings occur at g 2 (L) 7 but extrapolate to a slightly larger value g 2 ≈ 7.3 in the continuum limit. With c = 0.3 the crossings all cluster around g 2 (L) 7.3, with a nearly constant continuum extrapolation. Instead of taking this approach, however, in this work we construct the full continuum-extrapolated discrete β function across a broad range of couplings, the topic to which we now turn.

Step-scaling analyses and results
Following the standard procedure for lattice step-scaling analyses, for each L we first fit the renormalized couplings g 2 c (L) to some interpolating function in the bare coupling β F ≡ 12/g 2 0 , then use those interpolations to determine the finite-volume discrete β functions β s ( g 2 c , L) from eq. 2.4, which we extrapolate to the (a/L) → 0 limit. We will refer to the last step as the 'continuum extrapolation', although this is strictly true only for couplings weaker than the g 2 of the IR fixed point. While the choice of interpolating function is essentially arbitrary, typically some functional form motivated by lattice perturbation theory is used. For example, refs. [24,78,88] to polynomials in g 2 0 . Following refs. [77,95] we instead use the rational function which also produces the expected g 2 c ∝ g 2 0 at weak coupling. These interpolations are shown in figure 1. Most of the fits shown are of good quality, although there are some outliers with χ 2 /d.o.f. 1 corresponding to confidence levels CL 0.1. For reference we collect all this information in tables 11-13 in appendix C. Notably, the worst-quality interpolations are for the L = 12 data that we omit from our s = 3/2 step-scaling analyses.
To investigate potential systematic effects from our choice of interpolating function we also carry out analyses using [78] 1 where we include five terms to produce the same number of fit parameters as eq. 4.1. Although these interpolations appear satisfactory upon visual inspection, they generally produce much larger χ 2 /d.o.f. than the rational function in eq. 4.1 (tables 11-13). Therefore we will use the rational function for our final results, and treat any statistically significant differences between these two analyses as another source of systematic uncertainty.
Turning to the (a/L) → 0 extrapolations, we show several representative extrapolations in figure 2, for c = 0.25 and 0.3 at two values of the renormalized coupling u = 4 and 8 on either side of the IR fixed point. Since staggered fermions are O(a) improved, we extrapolate linearly in (a/L) 2 . In each figure we compare results from the clover discretization of the energy density E(t) in eq. 2.3 for several values of the t-shift improvement parameter τ 0 , including τ 0 = 0 and the optimal τ opt = 0.08. We also include one set of results from the plaquette discretization of E(t), at the corresponding optimal τ (plaq) opt = 0.02. We use the same vertical scale for both c = 0.25 and 0.3, to illustrate how the larger value of c reduces the size of cutoff effects for fixed τ 0 , as expected [90].
The unshifted (τ 0 = 0) results in figure 2 all show significant dependence on (a/L) 2 , despite the tree-level perturbative correction discussed in section 2. We optimize τ 0 by finding the value τ opt for which these cutoff effects are minimized. Since we use constant τ opt for all couplings, at most values of u the O(a 2 ) effects are only reduced and not entirely removed. For both c = 0.25 and 0.3 we find that τ opt = 0.08 (0.02) for the clover (plaquette) discretization of E(t) is satisfactory for the full range of couplings we consider. Figure 2 demonstrates the resulting reduction of cutoff effects on both sides of the IR fixed point.
At u = 4 the expected linear dependence on (a/L) 2 provides a good description of the data for L ≥ 16, with average confidence levels of 0.70 for c = 0.25 and 0.58 for c = 0.3. However, the L = 12 → 18 points clearly deviate from this linear scaling, which is our motivation for omitting these data from our main analyses. Figure 3 illustrates the effects of the L = 12 → 18 data on the quality of the (a/L) 2 → 0 extrapolations, by plotting the resulting χ 2 /d.o.f. for the full range of u that we access. While the larger c = 0.3 improves the quality of the extrapolations as expected [90], for most u the dominant contribution to the χ 2 comes from the L = 12 → 18 point. The exception is the region at stronger couplings u 7, where figure 2 suggests that the L = 16 → 24 points start to deviate from the larger-volume results. To account for this effect we repeat all continuum extrapolations with only the two points involving L ≥ 20, and include any differences between these results and the full L ≥ 16 prediction as another systematic uncertainty. From figure 2 we note that dropping L = 16 at strong coupling will produce (a/L) 2 → 0 extrapolations farther below zero, reinforcing the existence of the IR fixed point.
Ref. [33] comments that 'Symanzik-type' continuum extrapolations of the form shown in figure 2 -employing polynomials in (a/L) 2 -are guaranteed to be valid only in the basin of attraction of the gaussian UV fixed point, and not necessarily in the vicinity of the non-trivial IR fixed point. Our improvement of the gradient flow running coupling, discussed in section 2, addresses this issue. First, for any u we can find a value of the t-shift τ 0 for which the extrapolation is independent of L and therefore insensitive to the power of (a/L) in the extrapolation. Then, by demanding that all τ 0 produce the same result upon extrapolating (a/L) 2 → 0 we can check the validity of these extrapolations, and include any deviations as a systematic uncertainty. In this context, it is interesting to note that the resulting systematic uncertainties often increase significantly at couplings comparable to and stronger than g 2 (cf. figure 11 in appendix B), which may be related to this underlying issue.
So far we have discussed three potential sources of systematic error that we account for in our analyses. For convenience we briefly summarize them here: Interpolation: We interpolate g 2 c (L) as functions of β F on each lattice volume, fitting the data to both a rational function (eq. 4.1) and a polynomial (eq. 4.2). We take our final results from the rational function interpolations, and include any discrepancies between the two approaches as a systematic error. For c = 0.25 and intermediate u ≈ 5-6 this is the source of the largest systematic uncertainty, which is comparable to the statistical uncertainty. For c = 0.3 the different interpolations are much more consistent. Extrapolation: To assess the stability of the linear (a/L) 2 → 0 extrapolations we repeat all analyses without including the smallest-volume L = 16 → 24 data, considering only 20 → 30 and 24 → 36 points. We take our final results from the three-point extrapolations, with another systematic uncertainty defined by any disagreement between the two-and three-point analyses. This systematic uncertainty is largest at our stronger couplings u 7, where it can be approximately 2.5 times the statistical uncertainty, for both c = 0.25 and 0.3. As we emphasized in figure 2, the larger volumes produce extrapolated results for β s farther below zero, reinforcing the existence of the IR fixed point. Optimization: Finally, we account for any sensitivity to the t-shift improvement parameter τ 0 . Recall from section 2 that different values of τ 0 should all produce the same β s ( g 2 c ) in the continuum limit. Whenever our final results using the optimal τ opt differ from the results we would have obtained from unshifted (τ 0 = 0) analyses, we include the difference as a third systematic error. This is a conservative prescription, because we introduced the t-shift improvement to remove these cutoff artifacts, by enabling more reliable (a/L) → 0 extrapolations. Even so, this systematic uncertainty vanishes for all the s = 3/2 analyses considered in the body of this paper, which involve c ≥ 0.25 and L ≥ 16. In appendices A and B we report that this is not the case for some supplemental checks that include L = 12 data. Including L = 12, this systematic uncertainty vanishes only for c ≥ 0.3, and can even be the largest source of uncertainty if we consider the small c = 0.2 analyzed by refs. [31,34]. In all three cases, to ensure that statistical fluctuations are not double-counted as both statistical and systematic errors we take the latter to vanish when the results being compared agree within 1σ statistical uncertainties. We carry out separate error analyses for each of the clover and plaquette discretizations of the energy density E(t) in eq. 2.3. Additional systematic effects from the choice of E(t) discretization can be assessed by comparing the two sets of numerical results that we include in figure 4. We now present our final results for the 12-flavor system in figure 4, which shows the continuum-extrapolated s = 3/2 discrete β function for two different renormalization schemes, c = 0.25 and 0.3. In each panel we include our non-perturbative results for both the clover and plaquette discretizations of the energy density E(t) in eq. 2.3. Statistical uncertainties are shown by the darker error bands, while the lighter error bands indicate the total uncertainties, with statistical and systematic errors added in quadrature.
Along with our numerical results, figure 4 also shows the two-and four-loop perturbative predictions for the s = 3/2 discrete β function. These perturbative predictions are based on for N f fermions transforming in representation R of the gauge group. For the fundamental representation of SU(3) gauge theory we have 4) so that N f = 12 gives b 0 = 3 and b 1 = −50. Higher-order coefficients b i depend on the renormalization scheme. In the MS scheme, ref. [9] reports numerical values b 2 ≈ −1060 and b 3 ≈ 6808 for 12-flavor SU(3) gauge theory (see also ref. [10]). For most g 2 c our results in figure 4 lie in between the two-and four-loop perturbative curves, both of which predict an IR fixed point. At the weakest couplings we explore our results agree with the fourloop prediction, which remains slightly below the two-loop value. Since the discrete β function is scheme dependent these various results do not need to agree at non-zero u. Our comparisons with perturbation theory are for illustration only.
Finally, in figure 5 we compare our new results with the two recent large-scale stepscaling projects discussed in section 1 [33,34]. We overlay our c = 0.25 and 0.3 results from figure 4, adding c = 0.45 results from ref. [33] and c = 0.2 results from ref. [34], all using the clover discretization of E(t). Both of the latter analyses employ scale change s = 2 rather than the s = 3/2 that we use. Considering that all four sets of numerical results in figure 5 use different renormalization schemes, they are in good agreement throughout the range of couplings where they overlap. Had refs. [33,34] been able to explore the stronger couplings u 8 that we reach, we expect that they would have observed the same IR fixed point that we report. By coincidence, our IRFP is located at the same g 2 = 7.26 for both c = 0. 25

The leading irrelevant critical exponent
Now that we have observed a robust IR fixed point at g 2 = 7.26, we will extract the universal critical exponent related to the slope of the discrete β function at this IRFP. Linearizing β(g 2 ) ≈ γ g g 2 − g 2 around the fixed point, eq. 4.3 implies  Raw data for finite-size scaling analyses of the critical exponent γ g . The scaling relation in eq. 5.3 corresponds to straight lines on this log-log plot of | g 2 c − g 2 | vs. L. For both the plaquette discretization of E(t) at c = 0.25 (left, with τ (plaq) opt = 0.02) and the clover discretization at c = 0.3 (right, with τ opt = 0.08) we see g 2 c increase towards g 2 = 7.26 as the bare coupling increases from β F = 5.5 to 4.75 (empty symbols), then move to even stronger renormalized couplings for 4.25 ≤ β F ≤ 3.75 (filled symbols). Around β F ≈ 4.5 the signal effectively vanishes since g 2 c is so close to g 2 for all L. The other combinations of c and E(t) discretizations produce similar figures.
where ∆ ≡ g 2 (sL) − g 2 (L) = β s log(s 2 ) from eq. 2.4. Solving for the discrete β function allows us to relate its slope at the IRFP to γ g , Our convention in eq. 4.3 of considering the RG flow from the UV to the IR, L → sL, produces both β s < 0 and γ g < 0. We omit this negative sign to simplify comparisons with continuum predictions. Figure 4 already shows that we should obtain results comparable to four-loop perturbation theory in the MS scheme, which predicts γ g = 0.282 about 20% smaller than the two-loop result γ g = 0.360. A recent scheme-independent estimate γ g = 0.228 from ref. [15] is somewhat smaller still. Directly fitting the data shown in figure 4 to a linear form around g 2 = 7.26 produces c = 0.25 c = 0.3 Clover γ g = 0.251(4) γ g = 0.278(3) Plaquette γ g = 0.240(4) γ g = 0.270 (3) The uncertainties shown above include only the systematic effect of varying the fit range from g 2 ± 0.1 to g 2 ± 0.5. The high degree of correlation evident in figure 4 makes it challenging to determine meaningful statistical uncertainties from these fits. Since both observables as well as the c = 0.25 and 0.3 renormalization schemes should produce the same universal critical exponent, we can estimate additional systematic uncertainties from the spread in the numbers above. If we assume that these systematic effects dominate over the statistical uncertainties, then we end up with γ g = 0.26 (2). Alternately, we can carry out a finite-size scaling analysis to determine γ g , as in refs. [24,33]. The basic scaling relation is = 0.02) and clover discretization at c = 0.3 (right, with τ opt = 0.08). As expected, the signal vanishes around β F ≈ 4.5 where g 2 c is close to g 2 for all L, but the results are clearly consistent with the more precise predictions for γ g from the slopes of the discrete β functions in figure 4 (dashed lines). In each plot the three curves correspond to the central value of g 2 = 7.26 (green crosses) plus the minimum and maximum values of g 2 consistent with the combined statistical and systematic errors (red circles and blue squares, respectively).
for fixed bare coupling β F . In principle we could attempt to extract both g 2 and γ g from these fits, but to simplify the analysis we will use as input our determination of g 2 from figure 4. In figure 6 we show some of the data available to be analyzed, plotting | g 2 c (β F , L) − g 2 | vs. L on log-log axes for the c = 0.25 plaquette discretization and c = 0.3 clover discretization. The other two data sets (c = 0.25 clover and c = 0.3 plaquette) are similar. In all cases we can see g 2 c (β F , L) passing through the fixed-point g 2 = 7.26 around β F ≈ 4.5, causing the signal to vanish.
The finite-size scaling analysis amounts to linear fits of these data, the slopes of which correspond to γ g . Several significant systematic effects are visible in figure 6. First we can see that the slopes of linear fits will change slightly for different bare couplings β F . The scaling relation becomes more accurate closer to the IR fixed point, but the slow evolution of the coupling with L (figure 1) means that near the IRFP the signal in | g 2 c − g 2 | effectively vanishes for all L. Next, the slopes also depend on the range of L included in the fits. Empirically, we find that omitting the L = 12 data significantly increases the confidence levels of the fits. Additionally omitting L = 16 also tends to improve fit quality, while there are no obvious trends upon omitting larger L. Therefore we fit only L ≥ 18, and should account for any sensitivity to the fit range as a systematic uncertainty. We can also expect some systematic dependence on c and the E(t) discretization, as in the inline table above, which should be included in the final uncertainties as well. Finally, and most significantly, we obtain figure 6 by fixing g 2 = 7.26. Allowing g 2 to vary within the total uncertainties determined in the previous section leads to very wide variations in the resulting γ g .
In combination, these systematic uncertainties only allow us to use the finite-size scaling analysis as a consistency check on the value γ g = 0.26(2) determined directly from the slopes of the discrete β functions. This is shown in figure 7, where we plot finite-size scaling results for the critical exponent vs. the bare coupling β F , considering the same data sets shown in figure 6. In order to fill in more values of β F we interpolate these data, using the rational function discussed in section 4 (eq. 4.1). We see that the finite-size scaling results for fixed g 2 = 7.26 are clearly consistent with the γ g obtained from the corresponding β s (shown as dashed lines). As expected, the fit uncertainties blow up around β F ≈ 4.5 where the signal in | g 2 c − g 2 | effectively vanishes. Accounting for the uncertainties on g 2 produces the other two curves in each plot. Although the systematic spread of the results is enormous around the IRFP, the uncertainties are more manageable for β F 5, where they show a steady evolution towards the γ g = 0.26(2) determined above.

Discussion and conclusions
We have presented our final results for step-scaling calculations of the 12-flavor SU(3) discrete β function, using nHYP-smeared staggered fermions and an improved gradient flow running coupling. In the gradient flow scheme with c = 0.25 we observe an IR fixed point at g 2 = 7.3 +3 −2 , which changes to g 2 = 7.3 +5 −3 when c = 0.3. We are able to explore the discrete β function up to g 2 c 8.2, extending well past the IRFP, and account for systematic effects from the stability of the (a/L) → 0 extrapolations, the interpolation of g 2 c (L) as a function of the bare coupling, the improvement of the gradient flow running coupling, and the discretization of the energy density. These results, including systematic uncertainties, are collected in figure 4. At the IRFP we measure the leading irrelevant critical exponent to be γ g = 0.26 (2), comparable to perturbative estimates. This value for γ g comes from the slope of the discrete β function and we checked that it is consistent with a finite-size scaling analysis, even though the very slow running of the 12-flavor coupling makes finite-size scaling challenging for 12 ≤ L ≤ 36.
We have also shown (figure 5) that our results are consistent with the two recent largescale step-scaling projects discussed in section 1 [33,34], which were able to investigate only g 2 c 6.4. Ref. [33] emphasized the importance of comparing multiple discretizations of the energy density E(t) in the definition of the gradient flow running coupling (eq. 2.3), which motivated our investigation of both the plaquette-and clover-based observables. Considering L = 8 → 16, 10 → 20 and 12 → 24, ref. [33] found that c ≥ 0.45 was required to avoid systematic dependence on the choice of discretization. By moving to larger volumes L ≥ 16, we find good agreement between both discretizations for c ≥ 0.25. In appendix A we report that investigations including L = 12 need c ≥ 0.3 to obtain comparably good behavior. In particular, c = 0.2 analyses that include L = 12 data suffer from severe systematic uncertainties, which were not comprehensively considered in ref. [31] where we reported g 2 = 6.2 (2). With c = 0.2 and L ≥ 12 we now obtain a central value of g 2 = 5.88 and upper bound g 2 ≤ 7.13, which we discuss in appendix B (table 2 and figure 9).
Compared to perturbation theory, our results for the scheme-dependent g 2 lie in between the two-loop and four-loop MS values. At the weakest couplings we explore our s = 3/2 discrete β function agrees with the four-loop scheme, which remains slightly below the two-loop case. The scheme-independent critical exponent γ g = 0.26(2) that we obtain is consistent with the value 0.282 predicted by four-loop perturbation theory, which was also the case for the mass anomalous dimension γ m ≈ 0.235 found by refs. [56,58]. This close agreement with four-loop MS perturbation theory may be partly coincidental. Recent investigations of a scheme-independent series expansion [81] predict slightly different values γ g = 0.228 and γ m = 0.400 (5) [14,15], while an initial investigation of the five-loop MS β function [11] finds that the perturbative expansion breaks down at couplings weaker than g 2 , despite the apparently convergent behavior of the two-, three-and four-loop contributions. Although the extremely complicated computation of the five-loop β function is still being independently checked [96], subsequent analyses using it as input indicate that all systems with 9 ≤ N f ≤ 16 exhibit an IRFP [12][13][14][15].
The accumulating evidence for an IR fixed point in the discrete β function [23][24][25][26][27][28][29][30][31][32][33][34], in addition to further supporting evidence (summarized in section 1) from the phase diagram at zero and finite temperature [35,38,[41][42][43][44][45][46] as well as hyperscaling of the hadron masses and decay constants [53,56,58] increases our confidence in the conclusion that the 12-flavor system is conformal in the IR. The many existing investigations leave open a few directions that are particularly important to explore in the future. First, the existence of a conformal IRFP makes N f = 12 a useful basis for lattice studies of composite Higgs models in which the mass of some of the fermions is lifted to guarantee spontaneous chiral symmetry breaking [7,8]. Although there is some motivation for moving to a smaller N f 10 where the mass anomalous dimension may be larger, γ m 1, it is still advantageous to test this approach for N f = 12 where we have more information about the existence and characteristics of the IR fixed point. (There are relatively few lattice studies of the 10flavor system so far [97][98][99].) Finally, the fact that almost all 12-flavor lattice studies so far have employed staggered fermions makes it important to investigate the universality (or lack thereof) of the observed IRFP. As in three-dimensional spin systems [75,76], it is not guaranteed that different lattice fermion formulations with different chiral symmetry properties will produce identical predictions at a non-trivial fixed point. This provides particular motivation for studies using Ginsparg-Wilson (overlap or domain wall) fermions that possess continuum-like chiral symmetries, despite their increased computational cost.

A Results with different scale changes
As shown by table 1 in section 3, our data also allow us to carry out step-scaling analyses with scale changes s = 2 and 4/3 in addition to the s = 3/2 considered in the body of the paper, if we are willing to include the smallest lattice volume 12 4 . Following the same procedures described in section 4 produces the continuum-extrapolated discrete β function results shown in figure 8 for c = 0.25 and 0.3. While all of these analyses predict an IR fixed point consistent with that found for s = 3/2, the inclusion of the L = 12 data increases the systematic uncertainties, especially for the smaller s = 4/3 where the slow flow of the coupling is more difficult to resolve.
In particular, it is interesting to note that in the s = 2 case (L ≥ 12) where the uncertainties are better controlled, we need c ≥ 0.3 in order to obtain good agreement between results from the plaquette vs. clover discretizations of the energy density E(t) in eq. 2.3. This contrasts with the good agreement we observe even for c = 0.25 in figure 4 when considering only L ≥ 16. That is, larger lattice volumes improve the agreement between these two discretizations, which is consistent with expectations and with the results reported by ref. [33]: Considering L ≥ 8, ref. [33] found that c ≥ 0.45 was needed to obtain comparable agreement. One other notable change from the L ≥ 16 results in the body of the paper is that the systematic uncertainty due to t-shift optimization discussed in  figure 4 and also predicting an IR fixed point consistent with the s = 3/2 analyses considered in the body of the paper. The inclusion of L = 12 data in the analyses leads to larger systematic uncertainties, especially for the smaller s = 4/3 where the slow flow of the coupling is more difficult to resolve. section 4 no longer vanishes for c = 0.25. However, this systematic uncertainty continues to vanish for c = 0.3, suggesting that it -like the effect of E(t) discretization -is also sensitive to the combination of c and lattice volume.
From figure 8 we can again estimate the leading irrelevant critical exponent γ g from the slopes of the discrete β functions at the IRFP. (The finite-size scaling consistency check discussed in section 5 already included all of the data going into the s = 2 and 4/3 analyses.) Following the same procedure described in section 5 (i.e., neglecting statistical uncertainties and setting systematic uncertainties by demanding agreement for c = 0.25 and 0.3 with both plaquette and clover discretizations) produces γ g = 0.23(4) from s = 2 and γ g = 0.18(4) from s = 4/3. The former value agrees with our result γ g = 0.26(2) from s = 3/2 with L ≥ 16. The latter is 1.8σ smaller than that result, but its errors may need to be increased to reflect the larger systematic uncertainties in the s = 4/3 discrete β function. Both are smaller than the four-loop perturbative value 0.282, but are consistent with the scheme-independent 0.228 from ref. [15]. In summary, all scale changes s that we can consider with our data set consistently predict a 12-flavor IR fixed point and a leading irrelevant critical exponent comparable to perturbative estimates.
B Results with smaller c = 0.2 One advantage of the gradient flow running coupling is that it is straightforward to re-run analyses for an entire family of renormalization schemes parameterized by c = √ 8t/L. In general the renormalized coupling has smaller statistical uncertainties for smaller c, while larger c can help to reduce systematic effects [90]. We have already seen in figures 2 and 3 that c = 0.3 reduces cutoff effects and improves the quality of (a/L) 2 → 0 extrapolations compared to c = 0.25. In appendix A we discussed how analyses including L = 12 require c ≥ 0.3 in order to obtain good agreement between results employing the clover vs. plaquette discretizations of the energy density E(t) in eq. 2.3. This agreement persists even with c = 0.25 when L ≥ 16 as in the body of the paper, motivating our choice to focus on c = 0.25 and 0.3 for our main analyses.
However, since some previous works [31,34] used c = 0.2, here we consider what results our current data and analyses would produce in this scheme. Following the same procedures described in section 4 leads to the continuum-extrapolated discrete β function results shown in figure 9 for scale changes s = 3/2 (top), 2 (bottom left) and 4/3 (bottom right). In the top row of plots we contrast s = 3/2 analyses with L ≥ 16 as in the body of the paper (left), or L ≥ 12 as is required for the other scale changes (right). While the L ≥ 16 plot is well behaved and predicts an IR fixed point at g 2 = 6.93 +12 −56 , adding statistical and systematic uncertainties in quadrature, the combination of L = 12 and c = 0.2 dramatically increases the systematic uncertainties. Even though the other three analyses still produce an IR fixed point with the clover discretization, they prefer significantly smaller g 2 = 6.04, 5.88 and 5.60 for s = 2, 3/2 and 4/3, respectively, with significantly larger systematic uncertainties. In particular, the results based on the plaquette discretization of E(t) appear unreliable (giving β s < 0 for all couplings we can access), further reducing our confidence in these analyses. This is relevant since the result g 2 = 6.2(2) from ref. [31] came from using  Since we account for such discrepancies as a source of systematic error, an easier way to assess them is to inspect the 'error budgets' shown in figure 11. For each renormalized coupling u these plots show the statistical uncertainties and the three systematic uncertainties summarized in section 4, along with their combination in quadrature. (Recall from section 4 that we take systematic errors to vanish when their effects are indistinguishable from statistical fluctuations, to avoid double-counting the latter.) The top-right plot corresponds to one of the main analyses discussed in the body of the paper, with s = 3/2, c = 0.25 and L ≥ 16. As described in section 4, the optimization uncertainties vanish for all u, the interpolation uncertainties are comparable to the statistical uncertainties for intermediate u ≈ 5-6, and the extrapolation uncertainties dominate for stronger couplings u 7 (where the larger volumes L ≥ 20 would produce β s farther below zero). When we move to c = 0.2 in the top-left plot we see that the optimization uncertainties are now non-zero, in accordance with figure 10.
However, thanks to L ≥ 16, in the top-left plot of figure 11 the optimization uncertainties remain comparable to the statistical uncertainties. This is what changes when The key conclusion is that c = 0.2 was a poor choice in our previous L ≥ 12 study [31], which we have now corrected in this work by using both larger c and larger L.

C Data sets and interpolations
Tables 3-10 summarize the lattice ensembles considered in this work, with a separate table for each L = 12, 16, 18, 20, 24, 30, 32 and 36. In all cases we use exactly massless fermions with anti-periodic BCs in all four directions, while the gauge fields are periodic. For each ensemble specified by L and the bare coupling β F , the tables report results for the renormalized gradient flow couplings g 2 c (L; a) with τ 0 = 0 and g 2 c (L; a) with the optimal τ opt = 0.08, in both cases considering the clover discretization of E(t) for two values of c = 0.25 and 0.3. These results are obtained from the stated number of thermalized measurements, all of which are separated by ten molecular dynamics time units (MDTU) generated with the HMC algorithm, and combined into ten-measurement (100-MDTU) jackknife blocks to reduce autocorrelations. We also list the average plaquette (normalized to 3), to illustrate the roughness of the gauge fields.