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

We calculate the step scaling function, the lattice analog of the renormalization group $\beta$-function, for an SU(3) gauge theory with ten fundamental flavors. We present a detailed analysis including the study of systematic effects of our extensive data set generated with ten 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. In an accompanying paper we discuss that gradient flow can promote lattice dislocations to instanton-like objects, introducing nonperturbative lattice artifacts to the step scaling function. Motivated by the observation that Wilson flow sufficiently suppresses these artifacts, we choose Wilson flow with the Symanzik operator as our preferred analysis. We study systematic effects by calculating the step-scaling function based on alternative flows (Zeuthen or Symanzik), alternative operators (Wilson plaquette, clover), and also explore the effects of the perturbative tree-level improvement. Further we investigate the effects due to the finite value of $L_s$.


I. INTRODUCTION
Strongly coupled gauge-fermion systems play a central role in different beyond the Standard Model scenarios.
Two prominent examples are models for composite dark matter [1,2] and composite Higgs models [2][3][4][5][6]. The latter require a large scale separation to accommodate that only a light 125 GeV Higgs boson has been experimentally observed but no other heavier resonances. Hence composite Higgs models favor gauge-fermion systems featuring near-conformal dynamics, i.e. a slowly running (walking) gauge coupling. To identify promising candidate systems, it is essential to understand the nature of near-conformal models. Gauge-fermion systems are characterized by the gauge group G and the number N f of fermion flavors in representation R. Given these three characteristics, the infrared properties of gauge-fermion systems can be derived by determining the renormalization group (RG) β-function. Calculations based on perturbation theory predict that systems with SU(3) gauge group and fermions in the fundamental representation undergo a transition: as the number of flavors is increased, a chirally broken system with fast running coupling changes to a conformal system where the gauge coupling is irrelevant [7]. For even larger number of flavors, gauge-fermion systems become infrared free and possibly asymptotically safe [8][9][10].
While perturbation theory provides guidance about the N f dependence, the inherently nonperturbative nature of gauge-fermions systems limit the reliability of perturbative predictions. Particularly challenging is to identify the lowest, critical number of flavors N c f where a gauge-fermion system of gauge group G and representation R becomes conformal. It is known that for a system inside the conformal window, the infrared fixed point * Corresponding author: oliver.witzel@colorado.edu moves to stronger coupling as N f decreases towards N c f . Nonperturbative methods are required to determine N c f and gain field theoretical insight into how a chirally broken system changes to a conformal system.
Lattice field theory provides a nonperturbative framework to determine the RG β-function from first principles. We have studied the finite volume step-scaling function [11] for SU(3) with ten and twelve fundamental flavors [12][13][14][15] using domain wall (DW) fermions. Our work complements earlier studies of the eight and twelve flavor systems with staggered fermions [16][17][18][19][20][21]. Recently, we explored a new method, the continuous gradient flow β-function [22][23][24] presenting results for SU(3) with two and twelve fundamental flavors. In the infinite volume continuum limit both methods determine the renormalization group (RG) β-function, though in different renormalization schemes. Our recent DW results reveal that the two flavor system exhibits a fast running β function close to the perturbative 1-loop prediction, whereas for N f = 12 our step-scaling calculation shows that the β function is small in magnitude and identifies an IRFP in the range 5.2 ≤ g 2 c ≤ 6.4 using the c = 0.250 renormalization scheme.
In this work we present a detailed analysis of our stepscaling calculation for ten fundamental flavors. Compared to our results published in Refs. [14,15], we performed additional simulations at stronger bare couplings and added further volumes to improve the infinite volume continuum limit extrapolation. The additional simulations allowed us to increase the explored coupling range for c = 0.300 from g 2 c ≈ 6.5 in Refs. [14,15] to g 2 preferable Symanzik and Zeuthen flows introduce many more of these artifacts. In order to minimize this artifact, we choose Wilson flow as our preferred analysis.
Comparing the different nonperturbative lattice predictions in the c = 0.300 scheme, we find that our result is in perfect agreement at weak coupling (g 2 .0 range where we observe a down-turn of the β function pointing to a possible IRFP around g 2 c ∼ 13. Our nonperturbative results suggest that N f = 10 is likely conformal.
The bottom two panels of Fig. 1 show our continuum limit predictions in the c = 0.275 and 0.250 schemes. The results reveal that the GF step-scaling β-function exhibits a dependence on the renormalization scheme parameter c. However, cut-off effects on the finite volume step scaling function are more severe at smaller c. Unfortunately, our available data set does not allow to rigorously scrutinize our findings for c = 0.250 but we nevertheless include it in this publication for future reference.
In the remainder of this paper we present the details of our calculation and describe how we arrived at our final results shown in Fig. 1. Section II gives a brief description of the step-scaling function and Sec. III summarizes the numerical details of our simulations. We continue by discussing novel nonperturbative lattice artifact of the gradient flow in Sec. IV with a detailed investigation of this topology-related artifact presented in our companion work [25]. There we also demonstrate that nonzero topology is correlated with an increased value of the gradient flow coupling g 2 GF and the corresponding step-scaling βfunction. Afterwards we explain in Sec. V how we arrive our final results and how we check for systematic effects including a discussion on the finite value of the fifth dimension of domain wall fermion simulations. Finally, in Sec. VI we summarize our results and comment on the comparison with perturbative predictions. 1 We estimate the values of the LatHC result (gray band) for s = 2 based on Fig. 5 [26][27][28][29] and other lattice determinations [21,[30][31][32][33][34] in the c = 0.300 scheme. Lattice corrections due to small flow time in the c = 0.250 scheme could be significant, affecting the continuum limit shown on the last panel.

II. STEP-SCALING FUNCTION
The finite volume gradient flow coupling in volumes V = L 4 is where E(t) is the energy density, N = 3 for SU(3), and C(t, L/a) is a perturbatively computed tree-level improvement term 2 [35]. Without tree-level improvement C(c, L/a) is replaced by the term 1/(1 + δ(t/L 2 )) that compensates for the zero modes of the gauge fields in periodic volumes [11]. In the finite volume renormalization scheme the flow time t is set by renormalization scheme parameter c and the lattice size L t = (cL) 2 /8.
The discrete step scaling β function of scale change s is defined as in Ref. [11] β c,s (g 2 c ; L, β) = where g 2 c (L, β) = g 2 GF (t = (cL) 2 /8; L, β). The gradient flow transformation can be performed with different flow actions. In this work we consider Wilson, Symanzik and Zeuthen flows [36,37]. Similarly, the energy density E(t) at gradient flow time t can be approximated by different lattice operators. We consider the Wilson plaquette (W), Symanzik (S) and clover (C) operators. We use the shorthand notation [flow][operator] to refer to the various combinations. When the tree-level improvement term C(c, L/a) is included we add 'n' to our shorthand notation, e.g. nWS for Wilson flow, Symanzik operator, and tree-level improved coupling.
The renormalized coupling g 2 c is defined at a bare coupling β, therefore it is contaminated by cutoff effects. The infinite cutoff continuum limit requires t/a 2 → ∞ or equivalently L/a → ∞. At fixed value of g 2 c this means tuning the bare coupling g 2 0 = 6/β → 0, the Gaussian fixed point.
In practice we perform simulations at many values of the bare coupling β and combine them to cover the investigated range of the renormalized coupling. At fixed g 2 c we take the continuum limit (L/a → ∞) of the discrete step scaling function β c,s (g 2 c ; L) and obtain the continuum step scaling β-function β c,s (g 2 c ) in the renormalization scheme c.

III. NUMERICAL SIMULATION DETAILS
We determine the gradient flow step-scaling function on dynamical gauge field configurations with (L/a) 4 hypercubic volumes. As for our project with N f = 12 flavors, we choose to generate ten flavor configurations using tree-level improved Symanzik (Lüscher-Weisz) gauge action [38,39] and Möbius domain wall fermions (MDWF) [40] with three levels of stout-smearing [41] ( = 0.1) for the fermion action. 3 Our ensembles of gauge field configurations are generated using the hybrid Monte Carlo (HMC) update algorithm [47] with five massless two-flavor fermion fields and trajectories of length τ = 2 in molecular dynamics time units (MDTU). Simulations with domain wall fermions are performed by creating a five dimensional Dirac operator where the fifth dimension, L s , separates the chiral, physical modes of 4-d space-time. We choose L s = 12 for β ≥ 4.40 and L s = 16 for β ≤ 4.30. We set the domain wall height M 5 = 1 and simulate with bare fermion mass am f = 0.
The finite extent of the fifth dimension leads to residual breaking of chiral symmetry which conventionally is parametrized by an additive mass term am res . Although am res grows toward the strong coupling, we demonstrate in Sec. V C that our choice of L s is sufficient .
We set periodic boundary conditions (BC) for the gauge field and anti-periodic BC for the fermion fields in all four directions, i.e the fermion BC trigger a gap in the eigenvalues of the Dirac operator enabling simulations with zero input quark mass. To determine the step scaling β function we use hypercubic (L/a) 4 volumes with L/a = 8, 10,12,14,16,20,24,28, and 32 and generate for each volume a set of 17 bare couplings, starting in the weak coupling limit with β = 7.00 and choosing 4.02 as the value of our strongest bare coupling. Typically our ensembles consist of 6-10k for L/a ≤ 14, 3-6k for L = 16, 20, 24, and 2-4k for L/a = 28, 32 thermalized MDTU and we perform measurements separated by 10 MDTU. The simulations are carried out using GRID 4 [48] and we perform gradient flow measurements using Qlua 5 [49]. The subsequent data analysis is performed using the Γ-method [50] to estimate the integrated autocorrelation time and account for that as part of the statistical data analysis.

IV. TOPOLOGICAL ARTIFACTS OF THE GRADIENT FLOW
Before presenting details of our gradient flow calculation, we discuss a so far not considered lattice artifact related to topology. In this section we provide a brief description of the issue and refer for further details to our companion work [25].
Our numerical simulations are performed using massless, chirally symmetric domain wall fermions. With strictly massless fermions all configurations would have net topological charge zero. Configuration with nonvanishing topological charge would have zero Boltzmann weight because unpaired instantons create zero modes in the eigenvalue spectrum of the Dirac operator. If the topological charge is defined via the Dirac operator, Q ferm = 0. On the lattice topology is not a conserved quantity and the geometric definition of the charge can differ from Q ferm . Gradient flow removes the ultraviolet fluctuations of the gauge field and at large enough flow time Q geom is expected to be a good estimate of the topological charge.
In our simulations we monitor Q geom and observe Q geom ≈ 0 on configurations at relatively weak gauge coupling. However as the simulations are pushed into the strong coupling regime, Q geom = 0 configurations appear even at large flow time. We use the clover operator to estimate FF . In Fig. 2 we show the Monte Carlo (MC) histories of the topological charge Q geom determined after gradient flow time t/a 2 = 32 on our 32 4 ensembles at bare couplings β = 4.10, 4.05, and 4.02. Each of the three plots has three panels showing the results with Wilson flow (W, green), Zeuthen flow (Z, blue), and Symanzik flow (S, red). Even though the fermions cannot see unpaired instantons, apparently the gradient flow can promote vacuum fluctuations (dislocations) to instanton-like objects. This is an artifact of the gradient flow and not due to the action. In fact we have observed this effect in simulations with domain wall as well as with staggered fermions. Different gradient flows find different values for the topological charge on the same configuration, further showing that Q geom = 0 is an artifact of the flow. The fraction of nonzero topological charge configurations grows toward the strong coupling i.e. decreasing values of the bare coupling β. We also observe that Wilson flow has many fewer Q geom = 0 configurations than Symanzik or Zeuthen flow.
In Ref. [25] we correlated Q geom = 0 with the value of the gradient flow coupling and found that both the renormalized coupling and the step scaling function increases when Q geom = 0. Since the number of instantons is proportional to the square root of the volume, this nonperturbative artifact contributes a term to the step scaling function that increases with the volume. This invalidates the perturbatively motivated (a/L) 2 → 0 continuum limit extrapolation. The only way to control this artifact is to limit simulations to sufficiently weak coupling or use a gradient flow that suppresses Q geom . As is evident from In Sec. V B we revisit the effect of topology on the step scaling function when considering alternative flows.

V. GRADIENT FLOW β FUNCTION FOR TEN FUNDAMENTAL FLAVORS
We have calculated renormalized couplings g 2 c using three different gradient flows. However, we only consider Wilson flow to be acceptable at strong coupling and therefore obtain our preferred result using Wilson flow combined with the Symanzik operator. Due to too many configurations with nonzero topological charge, Zeuthen or Symanzik flows are solely used to estimate systematic effects.

A. Preferred (n)WS analysis
We present the values of the renormalized couplings  Table III. At weaker couplings, discretization effects can be reduced by applying perturbatively calculated tree-level normalization factors as included in Eq. (1). Since it is, however, not obvious if this perturbative correction is helpful at strong coupling, we carry out our preferred analysis with and without tree-level normalization and refer to the two analysis by nWS and WS, respectively. Our results show that for weaker couplings (g 2 c 6.0) tree-level normalization clearly removes discretization effects but also for stronger coupling leads to a prediction of the continuum limit step scaling function β c,s (g 2 c ) which is consistent with the WS determination without tree-level normalization. We therefore present our results for the ten flavor step scaling function in Figs. [3][4][5] showing the nWS analysis on the left and the and WS analysis on the right.
Our analysis proceeds in the following steps: • We calculate the discrete β c,s (g 2 c ; L) function defined in Eq. (3) for our five different volume pairs with scale change s = 2. The outcome is shown by the colored data symbols in the top panels of Figs. 3-5. Since simulations at different bare coupling β are statistically independent, also these data points are statistically independent.
• Next we interpolate the data for each pair of lattice volumes using a polynomial form motivated by the perturbative expansion We find that n = 4 is sufficient to describe our data well over the full g 2 c range, and the corresponding fits have all excellent p-values. Since discretization effects are sufficiently small at weak coupling for nWS, we constrain the intercept b 0 = 0 but fit b 0 for WS. We report the (correlated) fit coefficients in Table I. Some of the coefficients are not resolved statistically. This is not a concern, as long as the interpolations describe our data well. This is evident from the χ 2 /dof or p-values also given in Table I and visible by  • The interpolating fits provide us with finite volume discrete step scaling functions β c,s (g 2 c ; L) at continuous values of g 2 c . Due to the constraint b 0 = 0, nWS starts at zero but WS begins with our first data point. The range of g 2 c covered varies with the scheme c.
• To obtain β c,s (g 2 c ) in the continuum, we perform an infinite volume continuum limit extrapolation of the interpolated β c,s (g 2 c ; L) functions at fixed g 2 c values. Since in the strong coupling limit our smallest volume pair 8 → 16, and in some cases also 10 → 20 pair shows large discretization effects, we predict the continuum limit by performing a linear fit in a 2 /L 2 on the three largest volume pairs. In • The quality of the continuum limit extrapolation is monitored by the p-value shown by the second row panels in Figs. 3-5. While overall the p-values are excellent for most of the g 2 c range, the values rapidly drop for the weakest coupling. Theoretically not expected, we can only guess that poor pvalues at weak coupling are an artifact of very precise data. Since all five volume pairs sit very close to each other, there is little doubt on the continuum limit. Moreover, the weak coupling range is in good agreement with perturbative predictions which are certainly trustworthy for g 2 c 2.0. More concerning is the drop of the p-value for c = 0.250 in the range of 9 g 2 c 11. Here a low p-value around or even below 5% gives rise to concerns that the three volumes included in the continuum extrapolation may not be large enough [51].
• Further details of the continuum limit extrapolation are shown by the four panels at the bottom of Figs. 3-5. Selecting g 2 c = 3.0, 5.5, 9.0, and 10.8 we show the continuum limit extrapolation vs. (a/L) 2 . At weaker coupling (g 2 c = 3.0 and 5.5), our preferred linear fit for the three largest volume pairs is shown together with a quadratic fit using all five volume pairs. At stronger coupling (g 2 c = 9.0 and 10.8), the smallest volume pairs exhibit large cutoff effects and quadratic fits to all five volume pairs have very low p−values. Hence we show only a linear fit to the three largest data points. In all cases TABLE I. Results of the interpolation fits for the five lattice volume pairs for our preferred nWS (top half) and WS (bottom half) analysis using renormalization schemes c = 0.300, 0.275, and 0.250. Since discretization effects are sufficiently small for nWS, we constrain the constant term b0 = 0 in Eq. (5) and perform fits with 13 degrees of freedom (dof). For WS the intercept b0 is fitted and we have 12 dof. In addition we list the χ 2 /dof as well as the p-value.  (54) we observe excellent agreement between fits to the nWS and WS data.
Taking a closer look at the nWS analysis shown on the left of Figs. 3-5, it is noteworthy to point out that discretization effects are barely, if at all, visible for g 2 c 4.5. This range coincides with the value of g 2 where the perturbative predictions of the β-function at 3-, 4, and 5-loop order are also in good agreement (see e.g. Fig. 1).

B. Alternative flow/operators
Considering different operators and/or gradient flows can help to understand systematic effects. Compared to our N f = 12 calculation, our simulations for ten flavors extend to much stronger coupling. This has two consequences: on the one hand we established that Zeuthen and Symanzik flows exhibit novel lattice artifacts showing up as nonzero topological charge, on the other hand we know that discretization effects of different operators and gradient flows can also differ significantly. Hence care must be taken when considering alternative flows/operators in order to avoid that a particularly poor choice dominates the result.
To understand the effect of nonzero topological charge, we repeat our analysis for each flow/operator combination with a filtering option where we select only configurations with |Q| < 0.5, 6 i.e. we discard all configurations with topological charge |Q| > 0.5 at flow time t = (cL) 2 /8. Of course such an ad hoc measure is not theoretically sound. It gives, however, insight into the effect due to nonzero Q. For all three flows (Wilson, Zeuthen, Symanzik) we use the Symanzik operator in c = 0.275 scheme to calculate the continuum limit of the step-scaling function. The three plots in Fig. 6 show the overlay for each flow with and without filtering the topological charge. Further details of the Zeuthen flow analysis are presented in Figs. 11 and 12 in Appendix B. We would like to stress that for our preferred analysis based on Wilson flow filtering on Q has no effect, but filtering lowers the predicted step scaling function in the   In total we consider nine flow/operator combinations: Wilson (W), Symanzik (S), and Zeuthen (Z) flow each with three operators Wilson plaquette (W), Symanzik (S), or clover (C) and refer to the analysis by two capital letters indicating the flow and the operator. In addition we analyze the data without and with tree-level normalization, prefixing in the latter case 'n' to our short hand notation. Further, we perform the analysis with and without topological filtering. Choosing four selected values of g 2 c = 3.0, 5.5, 8.0, and 10.0 to cover the full range of g 2 c , we show the predicted values of β c,s (g 2 c ) in Fig. 7. The first row uses c = 0.300, the second row c = 0.275, and the last row c = 0.250. In all cases we show continuum limit results obtained from a linear fit in a 2 /L 2 to our three largest volume pairs. The shape of the symbol indicates the operator (square: Wilson plaquette, circle: Symanzik, triangle: clover) and the color distinguishes the flow (green: Wilson, blue: Zeuthen, red: Symanzik). At the strongest coupling (g 2 c = 10.0) shown, the effect of nonzero topological charge is most significant for Zeuthen and Symanzik flow. Hence the determination without filtering is shown in pale colors. In addition operators of   Focusing on determinations we consider reliable and displayed in vivid colors, we observe a very good consistency with our preferred (n)WS determinations indicated by the vertical green bands. Most determinations agree better than at the 1σ level and only the tree-level improved determinations with Symanzik flow (red) differ by more than 2σ. Moreover, Fig. 7 shows that systematic effects significantly grow at very strong coupling making a nonperturbative determination more and more challenging. The overall good consistency nevertheless bolsters confidence in our result.

C. Effect of finite Ls
As mentioned in Sec. III, domain wall fermions exhibit a small residual chiral symmetry breaking because for practical simulations the extent L s of the fifth dimension has to be finite. The residual chiral symmetry breaking is conventionally expressed in terms of an additive mass term am res . As part of our simulations we monitor am res by numerically determining it from the ra-   tio of the midpoint-pseudoscalar over the pseudoscalarpseudoscalar correlator. In Fig. 8 we show am res as a function of the bare coupling β for our simulations with ten and twelve dynamical flavors. Starting at β = 4.5 we observe a rapid growth of am res as β is decreased. While for our N f = 12 simulations we force am res < 5 · 10 −6 [12], this is not viable for N f = 10 where we intend to explore much stronger couplings. Instead we use L s = 16 for all simulations with bare coupling β ≤ 4.30 and generated additional ensembles at β = 4.05 using L s = 32 to verify that L s = 16 is indeed sufficient.
Choosing β = 4.05 has the advantage that Wilson flow still suppresses topological charges sufficiently well (see Fig. 2) and the effect of varying L s does not get obscured by nonzero topological charge contributions. Repeating our preferred WS analysis on the L s = 32 ensembles, we list in Table II the values for g 2 c (L; β) and β c,s (g 2 c ; L, β) for our four largest volume pairs using the three renormalization schemes c = 0.300, 0.275 and 0.250. Even though increasing L s from 16 to 32 leads to statistically resolved changes in g 2 c (L; β), the effect on β c,s (g 2 c ; L, β) is much less significant. The difference is even further attenuated when one compares β c,s (g 2 c ; L, β) vs. g 2 c (L; β) for L s = 32 to our default L s = 16 analysis as shown in Fig. 9. Over-   (L; β) and βc,s(g 2 c ; L, β) determined at β = 4.05 for our preferred WS analysis using ensembles with Ls = 16 and 32.

VI. CONCLUSION
Using gauge field configurations generated with stoutsmeared Möbius domain wall fermions and Symanzik gauge action, we have calculated the gradient flow stepscaling function for SU(3) with ten dynamical flavors. Our simulations explore the range of strong coupling so far not investigated in lattice calculations. Pursuing simulations in the range for g 2 c 8.0, we observe that the gradient flow occasionally promotes vacuum fluctuations (dislocations) to instanton-like objects. This is a lattice artifact that has not been described previously. The effect is more pronounced for some gradient flows than for others but always causes the gradient flow coupling to increase and run faster. Since Wilson flow does the best job in suppressing such dislocations compared to Zeuthen or Symanzik flow, we choose Wilson flow with Symanzik operator for our preferred analysis. Further we consider performing our analysis with and without treelevel normalization to reduce cutoff effects. Although justified only in the weak coupling limit, we find that our result with tree-level normalization is consistent to our unimproved result throughout the full range covered in g 2 c . Hence we quote, as shown in Fig. 10, the envelope covering our nWS and WS prediction as our final re-   the systematic effects. Using alternative flow/operator combinations to obtain a better estimate of systematic effects is however troublesome because of lattice artifacts induced by nonzero topological charge in the strong coupling regime. Discretization effects of some flow-operator combinations also grow substantially at strong coupling.
Moreover, we studied the effect due to the finite extent of L s which results in a small chiral symmetry breaking. Increasing L s from 16 to 32 at β = 4.05 we observe changes in g 2 c (L; β) which however mostly cancel in the difference β c,s (g 2 c ; L, β). In relation to our preferred analysis based on L s = 16 ensembles, the L s = 32 data are largely consistent with the interpolated L s = 16 result. This suggests that the overall effect due to the finite value of L s is negligible compared to other effects.
For all three values of the renormalization scheme c considered, we observe that the step-scaling β function exhibits a maximum in the range 8 g 2 c 9 and rapidly decreases for stronger coupling pointing to an IRFP at g 2 c 11. Although our analysis shows the β function changes sign in the c = 0.250 scheme, this result has to be taken with caution because poor p-values for the continuum extrapolation in that range may signal significant finite volume effects. Nevertheless, our results strongly suggest that the RG β function for N f = 10 exhibits an IRFP i.e. the SU(3) gauge-fermion system with ten flavors is likely conformal. This observation is fully in agreement with the spectrum analysis of the composite Higgs model featuring SU (3) with four light and six heavy flavors. There hyperscaling of ratios has been demonstrated [44][45][46] which also implies that N f = 10 is most likely conformal. Further, our data for the step-scaling β-function resolve a dependence on the renormalization scheme c.
As already shown in Fig. 1, our result is in excellent agreement with the determination by Chiu [30][31][32] for g 2 c 5.8 but significantly differ from his prediction for stronger coupling. At stronger coupling our result is however compatible with LatHC [21,33,34] who reached g 2 c ∼ 8.0. Although the nonperturbative gradient flow results correspond to a different renormalization scheme than the perturbative MS calculation, it is nevertheless instructive to compare to perturbative predictions [26][27][28][29]. Up to g 2 c ∼ 4, predictions at three, four, and five loop order are close. While 3-and 4-loop predict an IRFP around g 2 ∼ 10, the 5-loop MS result does not have a fixed point. Reference [28] suggested to improve the convergence of the perturbative series by using Padé approximation, and Ref. [52] has considered Borel re-summation to extend the perturbative range. Both references conclude that N f = 10 is so strongly coupled that these analytic calculations do not give a reliable result. Thus N f = 10 demonstrates the difficulties to obtain a perturbative prediction on the RG β-function and highlights the need for nonperturbative lattice calculations.

ACKNOWLEDGMENTS
We are very grateful to Peter Boyle, Guido Cossu, Anontin Portelli, and Azusa Yamaguchi who develop the Grid software library providing the basis of this work and who assisted us in installing and running Grid on different architectures and computing centers.  Fig. 11 shows the determination of the discrete beta function using Zeuthen flow and Symanzik operator both with and without tree-level improvement in the c = 0.300, 0.275, and 0.250 renormalization schemes. Fig 12 shows the result of the same analysis on topologically filtered |Q| < 0.5 data set. When filtering we discard all configurations with topological charge |Q geom | > 0.5 at flow time t = (cL) 2 /8. We compare the filtered and direct results in Fig. 10.