Gradient flow step-scaling function for SU(3) with twelve flavors

We calculate the step scaling function, the lattice analog of the renormalization group $\beta$-function, for an SU(3) gauge theory with twelve flavors. The gauge coupling of this system runs very slowly, which is reflected in a small step scaling function, making numerical simulations particularly challenging. We present a detailed analysis including the study of systematic effects of our extensive data set generated with twelve dynamical flavors using the Symanzik gauge action and three times stout smeared M\"obius domain wall fermions. Using up to $32^4$ volumes, we calculate renormalized couplings for different gradient flow schemes and determine the step-scaling $\beta$ function for a scale change $s=2$ on up to five different lattice volume pairs. Our preferred analysis is fully $O(a^2)$ Symanzik improved and uses Zeuthen flow combined with the Symanzik operator. We find an infrared fixed point within the range $5.2 \le g_c^2 \le 6.4$ in the $c=0.250$ finite volume gradient flow scheme. We account for systematic effects by calculating the step-scaling function based on alternative flows (Wilson or Symanzik) as well as operators (Wilson plaquette, clover) and also explore the effects of the perturbative tree-level improvement.


I. INTRODUCTION
The renormalization group β-function characterizes the nature of gauge-fermion systems with gauge group G and N f fermion flavors in representation R. Choosing e.g. fermions in the fundamental representation and the SU(3) gauge group, the system is chirally broken with a fast running coupling for a small number of flavors. A particularly well studied case is Quantum Chromodynamics (QCD). This is a chirally broken theory which exhibits only the Gaussian fixed point (GFP) at bare coupling g 2 0 = 0. When increasing the number of flavors, predictions based on perturbation theory suggest that a second, infrared fixed point (IRFP) develops at some definite g 2 IRFP [1]. For an even larger number of flavors, the theory becomes IR free. Theories exhibiting an IRFP are conformal. Of special interest is the lowest number of flavors, N c f , for which a gauge-fermion system (G, R) is conformal because this denotes the onset of the conformal window. 1 Perturbative methods to determine the nature of a system with N f flavors may not be reliable because the IRFP can occur at large g 2 IRFP and may lie outside the trustworthy region of perturbation theory. Inside the conformal window, the IRFP is expected to move to stronger couplings as the number flavors approaches N c f . Thus a perturbative estimate of N c f is particularly troublesome and warrants us to use nonperturbative methods [5][6][7][8].
In this work we focus on the SU(3) gauge system with twelve flavors in the fundamental representation and * oliver.witzel@colorado.edu 1 Theories just below the conformal window are promising candidates to describe physics beyond the Standard Model e.g. composite Higgs scenarios as discussed in Refs. [2][3][4] and references within.
present details of our nonperturbative, large-scale investigation using lattice field theory techniques. We evaluate the gradient flow (GF) step-scaling function, the lattice analogue of the renormalization group (RG) β-function in the infinite volume continuum limit. Perturbation theory up to the 4-loop level predicts that the system is conformal. The 5-loop β-function, however, does not predict a fixed point [5]. Since the 5-loop correction is very large, Refs. [7] suggested to improve the convergence of the perturbative series by using Padé approximation. Various forms of the Padé series predict that the system is conformal. The scheme-independent form [8] also finds an IRFP. There are several nonperturbative lattice studies of this system based on staggered fermions [9][10][11][12][13][14][15][16][17][18][19][20][21]. The results of lattice studies are controversial as Ref. [18] finds a FP in the gradient flow c = 0.25 scheme at g 2 ≈ 7.3, while Refs. [19][20][21] do not observe any sign of a FP. The origin of this disagreement is still to be resolved. Our calculation is the first based on simulations with domain wall (DW) fermions [22][23][24][25]. DW fermions are chiral and exhibit continuum-like flavor symmetry, a property which may be crucial to investigate near-conformal or conformal gauge fermion-systems. Furthermore, DW fermions are simulated with a Pauli-Villars term which reduces effective gauge terms generated by the many flavors, resulting in smoother gauge field configurations and reduced cutoff effects.
In our calculation we generate gauge field configuration with stout-smeared [26] Möbius domain wall fermions [25] and Symanzik gauge action [27,28] on which we perform gradient flow measurements. Our preferred determination of the step scaling function is based on Zeuthen flow [29,30] combined with the Symanzik operator to estimate the energy density. Other gradient flows (Wilson, Symanzik) as well as different operators (Wilson plaquette, clover) are considered, too. Our preferred combination of action, flow, and operator is fully O(a 2 ) Symanzik improved [29,30] and discretization effects are further suppressed by the perturbative tree-level normalization introduced in Ref. [31]. We demonstrate that this combination has indeed small discretization errors which result in mild continuum limit extrapolations. Using GF renormalization schemes c = 0.250, 0.275, and 0.300 we present continuum-limit results which show that SU (3) with twelve fundamental flavors exhibits an IRFP, implying that this theory is conformal. The IRFP may exhibit a mild dependence on the renormalization scheme which, however, is not statistically resolved by our data. In the c = 0.250 scheme we observe the IRFP in the range of 5.2 ≤ g 2 c ≤ 6.4. Early results of this project have been presented in Refs. [32][33][34][35]. The remainder of this paper is organized as follows: after introducing the gradient flow stepscaling function and the perturbative improvement in Sec. II, we discuss the details of our numerical simulations in Sec. III. In Section IV A we present the results for the step-scaling β-function using different renormalization schemes c for our preferred combination of Zeuthen flow with Symanzik operator with and without tree-level normalization (short: nZS and ZS). Alternative determinations using different gradient flows, operators with and without tree-level normalization are presented in Sec. IV B. Additional details are collected in the appendices: Appendix A lists the tree-level normalization factors for Symanzik gauge action and our three gradient flows with the three different operators, in Appendix B we present the renormalized couplings for our preferred (n)ZS analyses and in Appendix C we show comparison plots including additional continuum extrapolations. Finally we summarize our results and conclude in Sec. V.

II. THE CONTINUUM LIMIT OF THE STEP-SCALING FUNCTION
The finite volume step scaling function, or discrete β function of scale change s is defined by where g 2 c (L; a) is a renormalized coupling at the energy scale set by the volume µ = (cL) −1 which is related to the gradient flow time t by µ = 1/ √ 8t. The parameter c denotes the renormalization scheme corresponding to a renormalized gradient flow coupling given by with N = 3 for SU(3), E(t) the energy density at gradient flow time t, and C(c, L/a) a perturbatively computed tree-level improvement term 2 [31]. Without treelevel improvement C(c, L/a) is replaced by the term 2 Numerical values for C(c, L/a) are listed in Table III. 1/(1 + δ(t/L 2 )) that compensates for the zero modes of the gauge fields in periodic volumes [36].
The dependence on the lattice spacing a reflects cutoff effects that arise because simulations are not performed with a "perfect action," i.e. along the renormalized trajectory (RT) emerging from the Gaussian FP. As the gradient flow time t/a 2 increases, irrelevant operators die out and cutoff effects are reduced. Thus the continuum limit is approached by increasing L/a while keeping c and the renormalized coupling g 2 c fixed. In asymptotically free theories this forces the bare coupling toward zero (the Gaussian fixed point) or β ≡ 6/g 2 0 → ∞. If all but one irrelevant operators are negligible the remaining cutoff effects are proportional to ( √ t/a) α or (L/a) α where α is the scaling exponent of the least irrelevant operator. In the vicinity of the Gaussian FP α = −2, thus an extrapolation in a 2 /L 2 at fixed g 2 c predicts the continuum limit.
In practice simulations are performed at many bare coupling values, "daisy-chaining" them to cover the investigated renormalized coupling range. In slowly running ("walking") systems a very large scale change is required to cover even a moderate change in the renormalized coupling. Figure 1 shows g 2 c evaluated with Zeuthen flow, Symanzik operator and tree-level normalization in the c = 0.250 scheme as the function of the bare coupling on lattice volumes ranging from 8 4 to 32 4 . The 16 bare coupling values cover the range g 2 c ∈ (1.5, 6.5) in roughly uniform ∆g 2 c ≈ 0.4 increments. Since the N f = 12 coupling runs very slowly, prohibitively large volumes would be needed to reach the strong coupling regime as the bare coupling is tuned toward the GFP. When approaching an IRFP, β c,s → 0 and lattice volumes would have to grow without bounds. Nevertheless the L/a → ∞ limit predicts the infinite-cutoff step scaling function correctly, as long as the GF flow time is large enough to approach the vicinity of the RT. This reflects the fact that the renormalized trajectory describes continuum physics.
Cutoff effects scale with the scaling exponent of the least irrelevant operator. In the vicinity of the Gaussian FP α = −2, thus an extrapolation in a 2 /L 2 removes the cutoff effects at weak coupling. Along the RT α can however change and the a 2 /L 2 continuum extrapolation can become incorrect at stronger gauge couplings. Since it is very difficult to determine the unknown exponent α numerically, an elegant way out is to chose an improved setup where cutoff effects are suppressed and the infinite volume extrapolation is mild enough not to depend on the exact extrapolation form. We will show that our favored setup of Symanzik action, Zeuthen flow, and Symanzik operator has this property in the range of couplings covered by our simulations.
After extrapolating the discrete β-function β c,s (g 2 c ; L) to the continuum limit (i.e. L/a → ∞) at fixed g 2 c , the continuum β-function β c,s (g 2 c ) depends only on the renormalized coupling g 2 c . Therefore it is expected that β c,s (g 2 c ) is free of effects from irrelevant operators introduced by the lattice regularization and depends only on the renormalization scheme c and scale change s.

III. NUMERICAL SIMULATION DETAILS
Our nonperturbative determination of the gradient flow β function requires the generation of dynamical gauge field configurations on (L/a) 4 hypercubic volumes. Choosing tree-level improved Symanzik (Lüscher-Weisz) gauge action [27,28] and Möbius domain wall fermions (MDWF) [25] with three levels of stout-smearing [26] ( = 0.1) for the fermion action, 3 we generate ensembles of gauge field configurations using the hybrid Monte Carlo (HMC) update algorithm [39] with six massless two-flavor fermion fields and trajectories of length τ = 2 in molecular dynamics time units (MDTU). MDWF provide a prescription to simulate chiral fermions with continuum-like flavor symmetries by adding a fifth dimension to separate the chiral, physical modes of four dimensional space-time. In practice, the extent of the fifth dimension, L s , is finite which results in a residual chiral symmetry breaking, conventionally parametrized by an additive mass term am res . To study the β function for twelve flavors in the chiral limit, we set the input quark mass to zero and choose L s such that am res < 5 · 10 −6 . As can be seen in Fig. 2, the residual mass increases when approaching strong coupling but the numerically determined values of am res are largely independent of the four  32 ---32 ---dimensional volume. In our step-scaling calculation we therefore increase L s from 12 at weak coupling up to 32 for simulations at our strongest couplings to ensure that any effect from nonzero am res is negligible. The specific values of L s for each value of β and L/a are listed in Table I. At the strongest gauge couplings we performed additional simulations with alternative choices of L s to verify that the L s values listed in Table I Table V. In all cases we find that renormalized couplings for the same choice of L/a and β but different choices of L s agree at the 1σ level. We further found consistent results for β c,s (g 2 c , L) when substituting ensembles with different L s .
We set the domain wall height M 5 = 1 and have Pauli-Villars terms of mass one. Simulations are performed with the same boundary conditions (BC) in all four directions: periodic BC for the gauge field and antiperiodic BC for the fermion fields. The latter triggers a gap in the eigenvalues of the Dirac operator and thus allows simulations with zero input quark mass. In order to explore renormalized couplings up to g 2 c ∼ 6.5, we create ensembles for a set of bare couplings starting in the weak coupling limit with β = 7.00 and going down to β = 4.15 as our strongest coupling. 4 We choose hypercubic (L/a) 4 volumes with L/a = 8, 10,12,14,16,20,24,28, and 32 and typically generate 6-10k (2-4k) thermalized MDTU for most ensembles with L/a ≤ 24 (L/a = 28, 32). All ensembles are generated using GRID [40,41] and a gauge field configuration is saved every five trajectories (10 MDTU). In Fig. 3 we demonstrate that our combination of actions leads to sufficiently smooth gauge field configurations with the average plaquette, the smallest 1 × 1 Wilson loop (normalized to 1), ranging from 0.78 to about 0.59. Such values are comparable to those in QCD simulations and hence we expect the gauge fields to be sufficiently smooth to extrapolate to the correct continuum limit.
We perform gradient flow measurements on every fifth trajectory (10 MDTU). With the intention to investigate systematics of the gradient flow, we measure Wilson (W), Symanzik (S), and Zeuthen (Z) flow on each gauge field configuration using QLUA [42,43]. For all three flows, we estimate the energy density using three different operators: the Wilson-plaquette (W), clover (C), and Symanzik (S) operator. As described in Sec. II, the energy densities are proportional to the renormalized coupling which in turn lead to the step-scaling β function. In Appendix B we show the full details of our determinations of the energy densities for our preferred analyses based on Zeuthen flow with Symanzik operator (with and without tree-level normalization). Besides presenting g 2 c values in the renormalization schemes c = 0.250, 0.275, 4 For L/a = 8 and 16 we also simulated at β = 4.13. and 0.300, we also list the total number of measurements for each ensemble as well as the integrated autocorrelation time determined using the Γ-method [44].
In addition we use the gradient flow measurements to determine the topological charge and confirm it vanishes as expected for massless simulations.

IV. GRADIENT FLOW β FUNCTION FOR TWELVE FUNDAMENTAL FLAVORS
As discussed above, we have calculated energy densities using different gradient flows and operators for all of our ensembles of gauge field configurations. This allows us to determine the gradient flow β function in multiple ways and check for possible systematic effects.

A. Preferred (n)ZS analysis
Our preferred analysis is based on Zeuthen flow and the Symanzik operator. Since our ensembles are generated with Symanzik gauge action, this combination is fully O(a 2 ) improved and we indeed find small discretization effects. In the weak coupling limit, discretization effects can be further suppressed by applying the perturbatively calculated tree-level normalization factors (see Eq. (2)). A priori the range of validity in g 2 c for perturbatively calculated coefficients is not known. Therefore, we present the results for our preferred analysis with and without tree-level normalization. We refer to the two analysis as nZS and ZS, respectively. Our results are presented for the renormalization schemes c = 0.250, 0.275, and 0.300 in Figs. 4-6 and are obtained from the renormalized couplings listed in Appendix in B Table IV. Following Eq. (1), we calculate the discrete β c,s functions for five different volume pairs with scale change s = 2. The resulting values are denoted by the colored data symbols in the plots in the top panels of Figs. 4-6. By simulating a set of bare couplings β for all volumes, we directly obtain these statistically independent data points.
We interpolate values for each lattice volume pair using a polynomial Ansatz motivated by the perturbative expansion In practice we find that n = 3 is sufficient to describe our data and obtain fits with good p-value. In the case of nZS data, discretization effects are sufficiently small that we omit the constant term with coefficient b 0 . Looking at the data, it does not seem justified to force a zero intercept for ZS although numerically our fit does not resolve b 0 . We list the results of our interpolation in Table II where we also quote the χ 2 /dof as well as the pvalue of the fits. The p-values reflect that the fitted lines including statistical uncertainties pass through almost all   data points. Hence the fit Ansatz Eq. (3) leads to a good description of our data. Only the 8 → 16 volume pair corresponding to the smallest volumes used exhibits a low p-value around 5%. This data set suffers most from discretization effects and also exhibits a few "outliers" in the strong coupling limit which cause a larger χ 2 (lower p-value).
The interpolating fits predict finite volume discrete step-scaling functions β c,s (g 2 c ; L). These discrete stepscaling functions are shown in the top row panels of Figs. 4-6 by dashed lines with shaded error band in the same color as the corresponding data points. The inter-polation for nZS starts at g 2 c = 0 due to the constraint b 0 = 0 of the interpolating polynomial, whereas for ZS it begins with our data at the weakest coupling. The upper end of the interpolation depends on the range where we have data and varies slightly depending on the chosen scheme c.
Subsequently, we perform infinite volume continuum limit extrapolations using the interpolated finite-volume β-functions β c,s (g 2 c ; L), which are continuous in g 2 c . We consider two Ansätze for the extrapolation   • a quadratic fit in a 2 /L 2 to all five volume pairs, which are motivated by the expected O(a 2 ) discretization effects in the weak coupling limit. Similar to the use of tree-level normalization, the validity of an extrapolation proportional to a 2 /L 2 is limited to sufficiently weak couplings. Moreover correction terms of higher order may only be resolved if the statistical uncertainties are small enough. We therefore monitor the p-value of the extrapolation fits as a criteria to judge the validity of the continuum limit extrapolations. p-values are shown as a function of g 2 c in the plots in the second row panels of Figs. 4-6. In addition we show details of the continuum limit extrapolations at g 2 c = 2.0, 3.4, 4.8, and 6.2 in the four smaller plots at the bottom of the three figures. The resulting continuum step-scaling functions are shown in the top row plots by the solid black line with gray error band (linear fit) and the gray dash-dotted line (quadratic fit).
Overall most extrapolation exhibit excellent p-values over the full range in g 2 c covered by our simulations and corrections to an a 2 /L 2 extrapolation seem to be small. Performing a linear extrapolation in a 2 /L 2 to our three largest volume pairs is well justified for all data sets although a noticeable drop of the p-value can be observed   for nZS data around g 2 c ∼ 4.5 for all three c values considered. However, even for our strongest couplings the linear fit still gives a good p-value greater than 10%. In case of the ZS data, the p-values for linear extrapolations are mostly well above the 20% level. They do however show a similar behavior like in case of the nZS fits with lower values around g 2 c ∼ 2.9 and 5.5. The quality of the quadratic extrapolation based on all five lattice volume pairs is strongly affected by our smallest 8 → 16 volume pair. Discretization effects for the 8 → 16 data set are quite significant for c = 0.25 but decrease for c = 0.275 and 0.3. This corresponds to the well-known property that smaller c values require larger lattice volumes [45]. Using tree-level normalization, this effect is attenuated. For nZS, the p-value of the quadratic extrapolation even at c = 0.250 remains at the acceptable 4-5% level for the full range in g 2 c covered, whereas the quadratic extrapolation for ZS at c = 0.250 exhibits vanishing p-values around g 2 c ∼ 5.5.
Looking more closely at the details of the extrapolation plots presented at the bottom of Figs. 4-6, we observe that the nZS data for the three largest lattice volume pairs exhibit very little a 2 /L 2 dependence for all c values and in fact can also be described by fitting only a con-  stant [33]. For the smaller volumes at stronger couplings, corrections to a 2 /L 2 arise and the data show an upward curvature. In contrast to that the large volume ZS data clearly exhibit a slope in a 2 /L 2 for all c-values. An explanation could be that corrections at weak coupling are caused by higher order terms which can essentially be removed by the tree-level normalization, whereas the corrections at strong coupling indicate that perturbative improvement and extrapolations in a 2 /L 2 are subject to nonperturbative corrections.

B. Alternative flow/operators
In order to estimate systematic effects, we consider the step-scaling β function obtained with different flows (Wilson and Symanzik flow) as well as different operators (Wilson plaquette and clover). We repeat the same steps of the analysis as for our preferred (n)ZS data i.e. we first obtain data points for the discrete β functions for our five lattice volume pairs and next perform an interpolating fit. Again we constrain b 0 to be zero when using tree-level improvement but fit b 0 otherwise. Finally we again carry out an infinite volume continuum limit extrapolation using a linear Ansatz for the three largest volume pairs and a quadratic Ansatz using all five volume pairs. As examples for our alternative determinations, we present in Fig. 7 the determination for the step-scaling β function using Symanzik flow and the clover operator and in Fig. 8 Wilson flow with the Wilson plaquette operator. The plots on the left show the analysis with tree-level normalization factors, the plots on the right the analysis without tree-level improvement. Again we show results (top to bottom) for c = 0.250, 0.275, and 0.300. When using tree-level improvement the discretization effects be-tween the different lattice volume pairs are significantly reduced resulting also in a much smaller spread (smaller range of the ordinate) of the data sets. The continuum limit extrapolations are again shown by the solid line with gray error band (linear fit to the three largest volume pairs) and the gray dash-dotted line (quadratic fit to all five volume pairs). Details on the continuum limit extrapolations are presented in Appendix C in Figs. 12-14 where we show extrapolations for our preferred (n)ZS data in comparison to (n)SC and (n)WW.
While continuum limit predictions for the renormalized β c,s (g 2 c ) obtained from different flow/operator combinations are expected to agree, the finite volume predictions of β c,s (g 2 c , L) are subject to discretization effects and hence depend on the specific flow/operator combination. At finite volume, different flow/operator combinations predict different renormalized couplings even on the same set of gauge field ensembles evaluated using the same scheme c. Only after taking the continuum limit, the obtained β c,s (g 2 c ) are expected to be free of discretization effects and can be meaningfully compared. The accessible range of g 2 c and β c,s (g 2 c ) varies for each flow/operator combination at each bare coupling β ≡ 6/g 2 0 and (L/a) 4 volume. How different flow/operator combinations approach the same continuum limit can be seen by comparing Figs. 4-8. In total we obtain 36 different continuum limit predictions for the step-scaling function based on nine different flow/operator combinations, analyzing the data with and without tree-level normalization, and performing two different extrapolations to the 18 data sets. Choosing four representative couplings for the range of g 2 c simulated, we present in Fig. 9  At the weakest coupling shown (g 2 c = 2.0), the different analysis show very little scatter and in particular the tree-level improved results (bottom half) are well "aligned." Without improvement, differences can be mostly observed between linear and quadratic extrapolations. These are likely due to discretization effects on the smallest 8 → 16 volume pair. Looking at stronger coupling (g 2 c = 3.4 and 4.8), we observe that the "scatter" of the results grows. Zeuthen and Symanzik flow with tree-level improvement as well as Wilson flow without improvement are still very consistent. Only extrapolations of unimproved SC or WC combinations at c = 0.250 dif-  fer by more than 1σ from either of our preferred (n)ZS determinations. In both cases, the continuum extrapolation is less reliable because it covers a very large range (more than an order of magnitude) as can be seen e.g. in the upper right plot of Fig. 7. The larger scatter for the intermediate couplings could be related to the observation that most of our extrapolations around g 2 c ∼ 5 show a lower p-value, hence continuum predictions may be less reliable. At the strongest coupling, g 2 c = 6.2, statistical errors are significantly larger and within these all (available) results agree at the 1σ level.
Comparing the different step-scaling results for different c values, the overall consistency improves significantly when c increases. While at c = 0.250 our preferred determinations based on nZS and ZS exhibit a small tension at weak and intermediate couplings, both determinations are consistent at the 1σ level for c = 0.275 and 0.300. Being aware that c = 0.250 might be affected by discretization effects on our smaller lattice volume pairs and perturbative improvement may not be fully justified at the stronger couplings, we therefore choose to quote the envelope covering both, nZS and ZS, determinations as our final result which also accounts for systematic effects. Systematic effects due to different flow/operator combinations are not resolved within the combined nZS+ZS uncertainty. In fact only a few extrapolations for alternative flow/operators at c = 0.250 lead to results differing by more than 1σ. Such extrapolations cover a large range and are less reliable because extrapolating all five points using the quadratic Ansatz tends to differ by several sigmas from extrapolating the three largest volume pairs using a linear Ansatz. While the analyzed choice of flow/operator combinations is by no means complete, in fact an arbitrary number of operators and flows could be considered, we do expect to have studied a representative set. Finally, it is noteworthy that the agreement at weakest and strongest couplings is better than at the intermediate range.
Our final continuum limit results for the gradient flow step-scaling function obtained in renormalization schemes c = 0.250, 0.275, and 0.300 are shown in Fig. 10 where we also show perturbative predictions [5]. Our results are consistent with the perturbative predictions obtained at 3-or 4-loop in the MS scheme but differ significantly from the 5-loop, which grows rapidly for g 2 c > 4.0. Our prediction in the c = 0.250 scheme indicates a fixed point near the 4-loop MS value, whereas for c = 0.300 the 4-loop MS value coincides with the lower bound of our prediction of a possible fixed point Unfortunately, our uncertainties are too large to establish or rule out a dependence on the renormalization scheme c. Our predicted range of couplings for an IRFP of an SU(3) gauge theory with twelve flavors is consistent with predictions by Ryttov and Shrock who improve perturbative convergence using Padé approximants [7,48].

V. CONCLUSION
We have presented details of our gradient flow stepscaling calculation for SU(3) with twelve dynamical flavors. Our calculations are based on gauge field ensembles generated with Symanzik gauge action and three times stout-smeared Möbius domain wall fermions. Using Zeuthen flow combined with the Symanzik operator, we determine renormalized couplings and predict the step-scaling β-function. We assign a systematic uncertainty to our numerical predictions that covers the difference between the continuum limit extrapolated results of various flows and operators, different extrapolation forms, and also the effect of the tree-level normalization. For the c = 0.250 scheme we can identify a sign change of the β-function (infrared fixed point) in the range of 5.2 ≤ g 2 c ≤ 6.4, while for schemes c = 0.275 (c = 0.300) we can only name a lower limit for a possible fixed point 5.7 ≤ g 2 c (5.9 ≤ g 2 c ). It appears that the step scaling function and the value of the fixed point exhibit a mild dependence on the renormalization scheme c, however, further data at even stronger couplings with greater precision are required to resolve such a dependence. The set of available gauge field ensembles restricts the maximum g 2 c value a given flow/operator combination can reach. Less than half of the 36 combinations we consider reach g 2 c = 6.2, which unfortunately limits the control of our systematic error estimate around g 2 c ∼ 6.0 where we observe the IRFP for c = 0.250. Since the reach in g 2 c depends on the flow/operator combination, the converse argument implies that on a given set of gauge field ensembles, an IRFP can be identified for some flow/operator combinations but missed for others.
Our presented result are consistent with perturbative predictions obtained at 3-or 4-loop in the MS scheme. Finally, we compare our continuum limit results to other, nonperturbative determinations published in the literature. So far only calculations based on staggered fermions have been performed [17][18][19][20][21]. The work by Hasenfratz and Schaich [18] extends from the weak couplings at g 2 c ∼ 2.3 to strong couplings g 2 c ∼ 8.2 identifying an IRFP around g 2 c ∼ 7.3 for both c = 0.25 and c = 0.3. In the weak coupling their prediction is consistent with the findings of Ref. [17] which however uses larger c values; for g 2 c 6.0 it is in tension with the determination by the LatHC collaboration [19][20][21]. In Fig. 11 we compare our DWF result 6 to the prediction of the stepscaling function by Hasenfratz/Schaich [18] and LatHC [19][20][21]. Reference [18] presents results for the renormalization schemes c = 0.250 and 0.300 and accounts for systematic effects due to the extra-and interpolations. References [19][20][21] only use c = 0.250 without quantifying systematic effects. Extending what has so far been done in the literature, we emphasize that our results account for the systematic effect estimated by considering tree-level normalization, different flow/operator combinations, and continuum limit extrapolations.
Our results exhibit a similar shape as the β-function predicted by Hasenfratz and Schaich in Ref. [18]. While for c = 0.250, the two predictions differ by ∼ 2σ for intermediate and strong couplings, the differences almost vanish for c = 0.300. This might suggest that uncertainties in particular for c = 0.250 are underestimated. In contrast, the results of Refs. [19][20][21] predict a qualitatively different step scaling function that is nearly constant in a wide g 2 c coupling range, without a sign of an IRFP. Further investigations are needed to track down the origin of this disagreement. Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy and the RMACC Summit supercomputer [49], which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236), the University of Colorado Boulder, and Colorado State University. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562 [50] through allocation TG-PHY180005 on the XSEDE resource stampede2. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. We thank Fermilab, Jefferson Lab, NERSC, the University of Colorado Boulder, Texas Advanced Computing Center, the NSF, and the U.S. DOE for providing the facilities essential for the completion of this work.    (14) 3 (1) 3.433 (17) 3.434 (17) 3 (1) 3.440 (21) 3.441 (21) 3(1) 32 5.00 339 3.062 (17) 3.064 (17) 5 (2) 3.080 (22) 3.081 (22)  5(2) 3.099 (28)  For easier comparison we also list and highlight in bold face the ensembles used in our main analysis. Ensembles are specified by the spatial extent L/a, the bare gauge coupling β, and the value of the fifth dimension Ls. As above we list N , the number of measurements, as well as the renormalized couplings g 2 c for the analysis with (nZS) and without tree-level improvement (ZS) for the three renormalization schemes c = 0.250, 0.275 and 0.300. In addition the integrated autocorrelation times determined using the Γ-method [44] are listed in units of 10 MDTU.  (7) 5.5108(98) 5.923 (11) 0.58 (7)       FIG. 14. Details of the continuum extrapolation for our preferred (n)ZS data in comparison to alternative determinations based on (n)SC and (n)WW for c = 0.300. Continuum limit values of βc,s and corresponding p-values of the extrapolation are quoted in the legend.