Extended investigation of the twelve-flavor $\beta$-function

We report new results from high precision analysis of an important BSM gauge theory with twelve massless fermion flavors in the fundamental representation of the SU(3) color gauge group. The range of the renormalized gauge coupling is extended from our earlier work {Fodor:2016zil} to probe the existence of an infrared fixed point (IRFP) in the $\beta$-function reported at two different locations, originally in {Cheng:2014jba} and at a new location in {Hasenfratz:2016dou}. We find no evidence for the IRFP of the $\beta$-function in the extended range of the renormalized gauge coupling, in disagreement with {Cheng:2014jba,Hasenfratz:2016dou}. New arguments to guard the existence of the IRFP remain unconvincing {Hasenfratz:2017mdh}, including recent claims of an IRFP with ten massless fermion flavors {Chiu:2016uui,Chiu:2017kza} which we also rule out. Predictions of the recently completed 5-loop QCD $\beta$-function for general flavor number are discussed in this context.


History of the IRFP with 12 massless fermions
A conformal infrared fixed point (IRFP) of the β-function was reported earlier with critical gauge coupling g 2 * ≈ 6.2 and interpreted as conformal behavior of the much studied BSM gauge theory with twelve massless fermions in the fundamental representation of the SU(3) color gauge group [2]. This result was claimed to confirm the original finding of the IRFP in [7,8].
In disagreement with [2,7,8], the IRFP was refuted in [1]. Recently, responding to the negative findings in [1], the authors of [2] moved the IRFP to a revised new location g 2 * ≈ 7 in [3]. The relocation of the IRFP followed the announcement of a new IRFP with ten massless fermion flavors in the fundamental representation of the SU(3) color gauge group [5,6]. The claim in [5,6] would imply that the theory with twelve flavors must also be conformal and the lower edge of the conformal window (CW) of multi-flavor BSM theories with fermions in the fundamental representation would be located below ten flavors. No trace of the reported IRFP with ten flavors was found from high precision simulations in large volumes [9]. Some related predictions of the recently completed 5-loop QCD β-function for general flavor number will be also discussed in this context.
Results are reported for the β-function from the analysis of high precision simulations in large volumes for n f = 12 flavors in two different renormalization schemes and two different implementations of the gauge field gradient flow on the lattice, providing convincing evidence for the non-existence of the IRFP in [3]. A preview from our forthcoming new publication on the n f = 10 β-function [9] is also added showing evidence for the non-existence of the IRFP published in this theory [5,6]. * Corresponding author Email address: jkuti@ucsd.edu (Julius Kuti)

Scale-dependent β-function on the lattice
The gradient flow based diffusion of the gauge fields on lattice configurations from Hybrid Monte Carlo (HMC) simulations became the method of choice for studying renormalization effects with great accuracy [10][11][12][13][14]. In particular, we introduced earlier the gradient flow based scale-dependent renormalized gauge coupling g 2 (L) where the scale is set by the linear size L of the finite volume [15]. This implementation is based on the gauge invariant trace of the non-Abelian quadratic field strength, E(t) = − 1 2 TrF µν F µν (t), renormalized as a composite operator at gradient flow time t on the gauge configurations and measured from the discretized lattice implementation, as in [11].
Following [15], we define the one-parameter family of renormalized non-perturbative gauge couplings where the volumedependent gradient flow time t(L) is set by a fixed value of c = √ 8t/L from the one-parameter family of renormalization schemes. The renormalized gauge coupling g 2 (L) is directly determined from E(t) = − 1 2 TrF µν F µν (t) on the gradient flow of the gauge field at a fixed value of c which defines the renormalization scheme. The renormalization schemes c = 0.20 and c = 0.25 used in our work are identical to what was used in [2,3] including periodic boundary conditions on gauge fields and anti-periodic boundary conditions on fermion fields in all four directions of the lattice.
A general method for the scale-dependent renormalized gauge coupling g 2 (L) was introduced earlier to probe the step β-function, defined as (g 2 (sL) − g 2 (L))/ log(s 2 ) for some preset finite scale change s in the linear physical size L of the fourdimensional volume in the continuum limit of lattice discretization [16]. In our adaptation of the step β-function staggered lattice fermions are used with stout smearing in the fermion Dirac operator. The implementation of the HMC evolution code is described in [1] together with further details on the lattice step function and its continuum limit. Identical procedures are followed here.

High precision twelve-flavor analysis in large volumes
In the continuum limit, the monotonic function g 2 (L) implies in any of the volume-dependent schemes that a selected value of the renormalized gauge coupling sets the physical size L measured in some particular dimensionful physical unit. Fixed physical size L on the lattice is equivalent to holding g 2 (L) fixed at some selected value as the lattice spacing a is varied and the fixed physical length L is held constant by the variation of the dimensionless linear scale L/a as the bare lattice coupling is tuned without changing the selected fixed value of the renormalized gauge coupling. The continuum limit of the βfunction at fixed g 2 (L) is obtained by a 2 /L 2 → 0 extrapolation of the residual cut-off dependence detected as powers of a 2 /L 2 in the step β-function at finite bare gauge couplings g 2 0 while the renormalized gauge coupling is held fixed.
In the convention we use, asymptotic freedom in the UV regime corresponds to a positive step β-function given by the perturbative loop expansion for small values of the renormalized coupling. In the infinitesimal derivative limit s → 1 the step β-function turns into the conventional one. It turns out that even at a step size as large as s = 2 the step β-function tracks the conventional continuum β-function with very good accuracy and we use step s = 2 throughout the analysis. If the conventional β-function of the theory possesses a conformal fixed point, the step β-function will have a zero at the critical gauge coupling g 2 * , independent of the scale L. Since the renormalized gauge coupling g 2 (L) is a monotonic function of L, the IRFP coupling g 2 * is reached in the L → ∞ limit.
3.1. SSC gradient flow with c=0.20 renormalization scheme Precision tuning of the bare gauge coupling 6/g 2 0 was used exclusively for each step calculation in [1] using the c = 0.20 renormalization scheme with s = 2 step size. The gauge field gradient flow was driven by the same tree-level improved Symanzik gauge action which generated the gauge configurations. The lattice implementation of the gauge field operator E(t) = − 1 2 TrF µν F µν (t) used the clover construction which is known to reduce cut-off effects in the gradient flow [11]. SSC designates the setup with Symanzik gauge action driving both the gauge field gradient flow and the HMC evolution code, together with the clover operator implementation of E(t). The three target groups A, B, C of the precision tuned run sets tested the IRFP with negative conclusions, in disagreement with what was reported in [2]. In this work we target the step β-function in an extended range of the renormalized gauge coupling to cover the interval where the relocated IRFP was recently reported [3]. For consistency check, we added to the analysis the WSC scheme in section 3.2 where the Symanzik action is replaced by the simple Wilson plaquette action to drive the gradient flow on the gauge configurations. Only the Wilson plaquette action was used in [2,3] to drive the gradient flow. Importantly, we also test in the new work the influence of a 4 /L 4 cutoff effects in extrapolations to the continuum β-function and extend the analysis to the c = 0.25 renormalization scheme in section 3.3.
In Table 1 results are shown for gauge ensembles from the three new target groups D, E, F of precision tuned run sets using identical c = 0.2 renormalization scheme with s = 2 steps, as before.  Table 1: With previously tuned bare gauge couplings g 2 0 , the final 26 precision tuned runs are tabulated with 13 tuned runs and 13 paired s = 2 steps. The D, E, F run sets target g 2 approximately at 6.6, 6.8, 7.0 respectively.
The 26 runs were grouped into pairs for each step where the lower L/a value was precisely tuned to the target value of the renormalized gauge coupling. The higher L/a at the doubled physical size determined the step β-function at finite lattice spacing. Precision tuning for g 2 0 of the 13 steps of the three targets eliminated the systematic uncertainty in the step β-function from model-dependent interpolation in the bare gauge coupling. Figure 1 shows the remarkable accuracy of tuning for the three targets on the per mille accuracy level, with similar accuracy level of the renormalized g 2 entries in Table 1. The statistical analysis of the renormalized gauge coupling and step β-function of the precision tuned runs followed [17] and used similar software. For each run, extended in length between 5,000 and 20,000 time units of molecular dynamics time, autocorrelation times were measured in two independent ways, using estimates from the autocorrelation function and from the Jackknife blocking procedure. Errors on the renormalized couplings were consistent from the two procedures and the one from autocorrelation functions is listed in Table 1. Each run went through thermalization segments which were not included in the analysis. For detection of residual thermalization effects the replica method of [17] was used in the analysis. All 26 runs passed Q value tests when mean values and statistical errors of the replica segments were compared for thermal and other variations.
The leading cutoff effects in the step function at finite bare coupling g 2 0 appear at a 2 /L 2 order with a 4 /L 4 and higher order corrections. In our earlier work [1] the linear dependence on a 2 /L 2 was detected and fitted without including higher order corrections in reaching the continuum limit of the β-function. The high precision of the data allows for testing the a 4 /L 4 correction term leading to small increases of the continuum βfunction within one standard deviation in the SSC setup. This is illustrated in Figure 2 for targets D and F.  In addition to targeted precision tuning all the runs from the 6 targets A-F can be combined with additional trial runs from the tuning procedure into a new extended analysis which can project renormalized gauge couplings and step functions at any location in the bracketed range by using simple polynomial interpolation. Two samples of all the high quality polynomial interpolations are shown in Figure 3. This procedure allowed us to include the WSC analysis and the c = 0.25 renormalization scheme in the new work.  The results shown in Figure 3 are quite remarkable. The largest 28 → 56 step probes the smallest value of a 2 /L 2 with a positive step β-function over the whole bracketed g 2 range, incompatible with a zero in the β-function as first evidence against the existence of the IRFP.
The main results of the SSC analysis of the c = 0.20 renormalization scheme with step size s = 2 are shown in Figure 4.   Figure 5. The IRFP with red bar for statistical error and grey bar for systematic estimate is from [3]. The lower panel, with slightly shifted data in cyan and magenta colors for visibility, compares results from precision tuning and the combined method with an extrapolated point predicted slightly outside the bracketed interpolation range. Linear 3-point fits in a 2 /L 2 were used for precision tuned targets A and B (4-point fits before in [1] ). Both plots display the 4-loop and recent 5-loop results referenced in the main text. tion of the steps from finite bare couplings g 2 0 to the continuum β-function includes the a 4 /L 4 cutoff effects with typical fits shown in Figure 5.
In the 4-loop approximation the n f = 12 theory has an IRFP which disappears in the 5-loop approximation. The 5-loop βfunction predicts the lower edge of the conformal window between n f = 12 and n f = 13. It will require further investigation to understand the potential significance of the apparent consistency between the 5-loop β-function and our simulation results in different renormalization schemes. The authors of [2][3][4] prefer to show only the 4-loop β-function which exhibits a zero at the g 2 location close to where the original IRFP was published in [2] before being relocated to g 2 ≈ 7 in the c = 0.20 renormalization scheme [3,4]. The 4-loop zero in the β-function is inconsistent with the new 5-loop results. Independently, and non-perturbatively, the reported IRFP is ruled out in our new analysis of the extended data set with overwhelming statistical significance as illustrated in Figure 4.

WSC gradient flow with c=0.20 renormalization scheme
The WSC setup in our analysis designates the replacement of the tree-level improved Symanzik action by the simple Wilson plaquette action to drive the gradient flow on the gauge configurations which were generated by the Symanzik action in the HMC evolution code. The a 4 /L 4 cutoff contamination plays a critically important role in establishing consistency between the WSC scheme and the much less cutoff contaminated SSC scheme in reaching the correct continuum β-function. Since only the Wilson plaquette action was used for the gradient flow in [2,3], the WSC analysis is useful for completeness and consistency checks. Figure 5 shows the consistency of the two gradient flows converging to the same continuum β-function in the c = 0.20 renormalization scheme with step size s = 2. The analysis of extrapolations to the continuum β-function included both the a 2 /L 2 term and the a 4 /L 4 term in the fitting procedure. It is quite remarkable that consistency is clearly established although the cutoff effects in the WSC gradient flow are almost an order of magnitude larger than in the SSC gradient flow.
The origin of the large WSC cutoff effects: The lattice implementation of composite operators gets renormalized along the gradient flow as a function of flow time t. Figure 6 shows the renormalization of the clover lattice implementation of the composite operator E(t) = − 1 2 TrF µν F µν (t) along the gradient flow, proportional to g 2 (t) we target.   Table 1 was used at all 3 values of c for SSC and WSC.
The renormalization of the operator E(t) goes through transient effects at short flow times before the cutoff effects get sufficiently renormalized. These transient effects are almost an order of magnitude larger along the Wilson flow for the same clover lattice implementation as used in the Symanzik flow. In our WSC approach, much larger stepped volumes would be needed to match the smaller cutoff effects in extrapolation of the SSC approach to the continuum β-function. This is clearly demonstrated in Figure 5 and in Figure 6. We do not see how the stronger cutoff dependence of the Wilson flow allows one to ignore the a 4 /L 4 effects of small lattices when gauge and fermion actions different from ours are used in generating configurations in the simulations, like in [4][5][6].

The c=0.25 renormalization scheme
We also analyzed our data set in the c = 0.25 renormalization scheme which is further away from the c = 0 infinite volume scheme and less sensitive to cutoff effects. A subset of the full c = 0.25 analysis is shown in Figure 7. We combine again all  the runs from precision tuning of the 6 targets A-F with additional trial runs from the tuning procedure into a new extended analysis with SSC and WSC gradient flow and the c = 0.25 renormalization scheme at step size s = 2. In this new analysis we can predict again renormalized gauge couplings and step functions at any location in the bracketed range by using simple polynomial interpolation. Although the cutoff effect of the renormalization is much bigger when the Wilson action is driving the gradient flow, the two consistently converge to the same continuum β-function. The β-function in Figure 8 is consistent with the c = 0.20 renormalization scheme and without any trace of the IRFP reported in [2,3] for the c = 0.25 renormalization scheme.

Ten-flavor preview with conclusions
We explored SSC and WSC gradient flows in the c = 0.20 and c = 0.25 renormalization schemes for the determination of the continuum β-function with twelve flavors of massless fermions. Consistent results from our high precision simulations in large volumes do not show any infrared fixed point reported at several locations in [2][3][4]. This disagreement requires resolution and closure. Arguments were presented in [4] that the different results should not be viewed as disagreements in the simulations and in their analysis but as evidence for the violation of universality in the staggered formulation. Supporting this argument, three examples were invoked in [4]: (a) showing conformal fixed point structures in 3D statistical models with several fixed points, (b) disagreement between sextet β-functions using Wilson fermions [23] and rooted staggered fermions [24], (c) new results finding an IRFP with ten flavors of massless domain wall fermions implies conformal behavior for twelve flavors as well.
We note our disagreements in response: (a) Since staggered fermions at n f = 12 are built on a UV fixed point at zero gauge coupling, relevant or marginal operators, like in the examples of the 3D statistical models in [4], cannot be added to the staggered lattice fermion action which has correct locality and universality properties. The explicit construction is well-known in the literature. Besides, the controversy between [2][3][4] and our work exists for the staggered formulation itself.
(b) The theoretical framework for the rooting procedure, without universality violation when the sextet gauge coupling is targeted in fixed physical volume, was explained in [24]. It has not been challenged since in any follow up to the original criticism [23].
(c) We recently completed the comprehensive analysis of the theory with ten massless fermion flavors in the fundamental representation of the SU(3) color gauge group. An important result is shown in Figure 9 for the preview of the β-function with SSC analysis of the gradient flow in the c = 0.25 renormalization scheme using step size s = 2 [9]. Every detail of the ten-flavor analysis followed closely the procedure we implemented and used here for the analysis of the twelve-flavor model. Since rooting was used for ten flavors with staggered fermions, we tested the validity of the theoretical argument presented in [24]. The Dirac spectrum was closely analyzed in the runs and the expected behavior of the quartet structure was confirmed.  Figure 9: The IRFP as marked in the plot is taken from [5,6], in complete disagreement with our analysis.
In conclusion, we find unacceptable the resolution of the problem as conjectured in [4].