Improved study of the $\beta$-function of $SU(3)$ gauge theory with $N_f = 10 $ massless domain-wall fermions

I perform an improved study of the $\beta$-function of $ SU(3) $ lattice gauge theory with $N_f=10$ massless optimal domain-wall fermions in the fundamental representation, which serves as a check to what extent the scenario in the previous work [arXiv:1603.08854; Proc. Sci. LATTICE2016 (2017) 228] is valid. In the finite-volume gradient flow scheme with $ c = \sqrt{8t}/L = 0.3 $, the renormalized couplings $g^2 (L,a) $ of four primary lattices ($ L/a = \{ 8, 10, 12, 16 \}$) are tuned (in $ 6/g_0^2 $) to the same $ g_c^2 $ with statistical error less than $0.5 \% $, in contrast to the previous work where $ g^2(L,a) $ were obtained by the cubic-spline interpolation. Then the renormalized couplings $ g^2(sL, a) $ of the scaled lattices ($ sL/a = \{16, 20, 24, 32\} $ with $s=2$) are computed at the same $ 6/g_0^2 $ of the corresponding primary lattices. Using the renormalized couplings of four lattice pairs $ (sL,L)/a = \{ (16,8), (20,10), (24,12), (32,16) \} $, the step-scaling $\beta$-function $ [g^2(sL,a) - g^2(L,a)]/\ln (s^2) $ is computed and extrapolated to the continuum limit $ \beta(s,g_c^2) $, as summarized in Table III. Based on the four data points of $ \beta(s,g_c^2) $ at $ g_c^2 = \{ 6.86(2), \ 6.92(3), \ 7.03(2), \ 7.16(2) \} $, I infer that the theory is infrared near-conformal, or conformal with the fixed-point $ g_*^2 = 7.55(36) $. This corrects the scenario in the previous work with $ g_*^2 \sim 7.0 $, and also suggests that the interpolation method cannot give a reliable determination of the $\beta$-function, especially in the regime close to the infrared fixed-point.

is computed and extrapolated to the continuum limit β(s, g 2 c ), as summarized in Table III. Based on the four data points of β(s, g 2 c ) at g 2 c = {6.86(2), 6.92(3), 7.03(2), 7.16(2)}, I infer that the theory is infrared near-conformal, or conformal with the fixed-point g 2 * = 7.55(36). This corrects the scenario in the previous work with g 2 * ∼ 7.0, and also suggests that the interpolation method cannot give a reliable determination of the β-function, especially in the regime close to the infrared fixed-point.

I. INTRODUCTION
In Refs. [1,2], I investigated the β-function of the SU (3) gauge theory with N f = 10 massless optimal domain-wall fermions in the fundamental representation. The motivation was to see whether this theory possesses a non-trivial infrared fixed point, which is not only a fundamental problem in quantum field theory, but also relevant to beyond the Standard Model scenarios with a composite Higgs boson. (For recent reviews, see, e.g., Refs. [3][4][5].) The results in Refs. [1,2] suggest that the theory might possess an infrared fixed point (IRFP) around g 2 c ∼ 7. However, the major systematic uncertainty in Refs. [1,2] was that interpolation was used to obtain the renormalized couplings g 2 (L, a) and g 2 (sL, a). This could lead to a large systematic error in the strong-coupling regime where the renormalized coupling varies rapidly with respect to the bare coupling g 0 (or 6/g 2 0 ), which in turn may give incorrect results for the step-scaling β-function β(s, a/L, g 2 ) = g 2 (sL, a) − g 2 (L, a) ln(s 2 ) , as well as its extrapolated value in the continuum limit (a/L → 0). The purpose of the present study is to eliminate this systematic uncertainty, by tuning 6/g 2 0 such that the renormalized couplings g 2 (L, a) of all primary lattices (L/a = 8, 10,12,16) have the same value with statistical error less than 0.5%. The tuning process implies that many simulations on the primary lattices have to be performed, which are rather challenging in terms of computing resources, time, and effort. After the value of 6/g 2 0 is determined for a chosen g 2 c = g 2 (L, a), the simulation on the scaled (s = 2) lattice is performed at the same 6/g 2 0 to obtain the renormalized coupling g 2 (sL, a). Since the results in Refs. [1,2] suggest that the theory may possess an infrared fixed point around g 2 c ∼ 7, four targeted values of g 2 c = {6.86(2), 6.92(3), 7.03(2), 7.16(2)} around g 2 c = 7.0 are chosen. Also, a point at g 2 c = 3.51(2) is picked to check whether the interpolation used in Refs. [1,2] works well in the regime where the renormalized coupling varies slowly with respect to the bare coupling. Moreover, in view of a recent study of the SU (3) gauge theory with N f = 10 massless domain-wall fermions [6] (which reported 2-3 standard deviations compared with the results of Ref. [2] for 4.5 < g 2 c < 6.0), a point at g 2 c = 5.25(2) is chosen to check whether the discrepancy is due to the systematic error of the interpolation used in Ref. [2]. All together, the targeted values of g 2 c in this study are g 2 c = {3.51(2), 5.25(2), 6.86(2), 6.92(3), 7.03(2), 7.16(2)}. The outline of this paper is as follows. In Section II, we describe our hybrid Monte Carlo (HMC) simulation of SU (3) gauge theory with N f = 10 massless optimal domainwall fermions, and summarize the residual masses of all gauge ensembles in Table   I. In Section III, we present our results for the renormalized couplings in the finitevolume gradient flow scheme with c = √ 8t/L = 0.3, for all gauge ensembles, as summarized in Table II. In Section IV, we perform the extrapolation of the step-scaling βfunction β(s, a/L, g 2 c ) to the continuum limit (a/L → 0) with the linear fit [A + B(a/L) 2 ] and quadratic fit [A + B(a/L) 2 + C(a/L) 4 ]. The results are summarized in Table III. In Section V, we perform the extrapolation of β(s, g 2 c ) with the linear fit, using four data points at g 2 c = {6.86(2), 6.92(3), 7.03(2), 7.16(2)}, and determine the IRFP g 2 * and the slope of β(s, g 2 ) at the IRFP. In Section VI, we determine the universal scaling exponent γ * g of the conventional β-function β(g 2 (µ)) in the continuum, with the input of the slope of β(s, g 2 ) at the IRFP. In Section VII, we summarize the results of this paper, and discuss the discrepancies between the results in this paper and those obtained with N f = 10 massless staggered fermions in a recent study [7].

II. GENERATION OF THE GAUGE ENSEMBLES
Since we are dealing with massless fermions, it is vital to use lattice fermions with exact chiral symmetry at finite lattice spacing (i.e., domain-wall [8] /overlap [9] fermions) with exactly the same flavor symmetry as their counterpart in the continuum. Theoretically, the effective four-dimensional lattice Dirac operator of the domain-wall fermion with infinite extent in the fifth dimension (N s = ∞) is exactly equal to the overlap Dirac where m q is the bare fermion mass, , m 0 ∈ (0, 2), The chiral symmetry can be maximally preserved on a lattice with finite N s by optimal domain-wall fermions [13], with the effective four-dimensional lattice Dirac operator exactly equal to the Zolotarev optimal rational approximation of the overlap Dirac operator. In this paper, we use optimal DWFs with the R 5 symmetry [14], whose effective four-dimensional lattice Dirac operator exactly equal to the "shifted" Zolotarev optimal rational approximation of the overlap operator, with the approximate sign function S(H) where d Z is the maximum deviation |1− √ xR Z (x)| max of the Zolotarev optimal rational polynomial R Z (x) of 1/ √ x for x ∈ [1, λ 2 max /λ 2 min ], with degrees (n − 1, n) for N s = 2n. The action of one-flavor optimal DWFs can be written as where the indices x and x ′ denote the sites on the four-dimensional spacetime lattice, s and s ′ are the indices in the fifth dimension, and the lattice spacing a and the Dirac and color indices have been suppressed. Here D w is the standard Wilson-Dirac operator minus the parameter m 0 ∈ (0, 2). The operator L is independent of the gauge field, and it can be written as where m q is the bare fermion mass, m 0 ∈ (0, 2), and N s is the number of sites in the fifth dimension, For massless DWFs, m q is set to zero. Besides Eq. (2), the action for the Pauli-Villars fields with m q = 2m 0 has to be included for the cancellation of the bulk modes, which is exactly the same as Eq.
(2) except for m q = 2m 0 in L ± [Eq. 3]. Thus the action for SU (3) lattice gauge theory with N f = 10 massless optimal DWFs can be written as where S g (U ) is the gauge action. In this paper, we use the Wilson plaquette gauge action where g 0 is the bare coupling. For the fermion action, we set m 0 = 1.8, and N s = 16. The optimal weights ω s [14] are computed with λ max /λ min = 6.2/0.05.
Simulating N f = 10 DWFs amounts to simulating five pairs of N f = 2 DWFs. Starting from the action (2) and following the procedures of even-odd preconditioning and the Schur decomposition given in Ref. [15], the partition function for the SU (3) gauge theory with N f = 10 massless optimal DWFs can be written as where φ i and φ † i are pseudofermion fields, and However, HMC simulations with Eq. (4) turn out to be rather time consuming for large lattices at strong couplings, e.g., 32 4 at 6/g 2 0 = 6.45. To resolve this difficulty, we use a novel N f = 2 pseudofermion action based on the exact pseudofermion action for one-flavor DWFs [16], which turns out to be more efficient than (4). This novel N f = 2 pseudofermion action for the optimal DWFs can be written as Here Note that K(m) is defined on the four-dimensional lattice, while H T (m) is a Hermitian operator defined on the five-dimensional lattice. The general form of the novel twoflavors pseudofermion action for domain-wall fermions with H = cγ 5 D w (1 + dD w ) −1 will be presented in a forthcoming paper [17].
Then, the partition function for the SU (3) gauge theory with N f = 10 massless optimal DWFs can be written as We perform HMC simulations of all gauge ensembles with Eq. (5)  The chiral symmetry breaking due to finite N s = 16 can be measured in terms of the residual mass of the massless fermion [21], where D −1 c denotes the massless fermion propagator, "tr" denotes the trace running over the color and Dirac indices, and the brackets · · · U denote averaging over all configura-tions of the gauge ensemble. The residual masses of all gauge ensembles in this work are summarized in Table I.
We observe that the variation of the residual mass is quite mild, ranging from ∼ 4.4 × 10 −5 to ∼ 8.5 × 10 −5 , i.e., less than a factor of 2. Moreover, the residual mass of any lattice size L 4 is much smaller than the energy scale µ ≃ (cL) −1 of the finite-volume gradient flow Even for the smallest µ of the largest lattice 32 4 in this work, the residual mass of any gauge ensemble satisfies Thus the effect of the residual masses on the renormalized couplings should be negligible for our analysis.

III. RENORMALIZED COUPLING OF THE FINITE-VOLUME GRADIENT FLOW SCHEME
To obtain the renormalized coupling of gauge theory on a finite lattice with volume L 4 , we use the finite-volume gradient flow scheme [22], which is based on the idea of continuous smearing [23] or equivalently the gradient flow [24] to evaluate the expectation value t 2 E , where E is the energy density of the gauge field and t is the flow time. This amounts to solving the discretized form of the following equation As shown in Ref. [24], the gradient flow is a process of averaging gauge field over a spherical region of root-mean-square radius R rms = √ 8t. Moreover, since t 2 E is proportional to the renormalized coupling, one can use c = √ 8t/L as a constant to define a renormalization scheme on a finite lattice, and obtain where a is the lattice spacing depending on the bare coupling g 0 , E is the energy density, and the numerical factor on the rhs of Eq. (6) is fixed such that g 2 c (L, a) = g 2 MS to  finite-lattice-spacing corrections [25]. In this paper, we use the Wilson flow, the Wilson action, and the clover observable, the so called WWC scheme, which is known to have very small tree-level cutoff effects [25]. Moreover, we fix c = √ 8t/L = 0.30. For each targeted value of g 2 c , the renormalized couplings g 2 (L, a) of four primary lattices (L/a = {8, 10, 12, 16}) are tuned (in 6/g 2 0 ) to the same g 2 c with statistical error less than 0.5%. Here the statistical error is estimated using the jackknife method with a bin size of 10-15, of which the statistical error saturates. In Fig. 1 After the value of 6/g 2 0 is determined for a chosen g 2 c = g 2 (L, a), the simulation on the scaled (s = 2) lattice is performed at the same 6/g 2 0 to obtain the renormalized coupling g 2 (sL, a). All renormalized couplings of g 2 (L, a) and g 2 (sL, a) are summarized in Table II. Each row gives the values of g 2 (L, a) and g 2 (sL, a) at the same 6/g 2 0 . Every four rows are grouped for the same targeted value of g 2 c .
Moreover, if β(s, g 2 c ) is determined for several values of s, then it can be extrapolated to where β(g 2 (µ)) is equal to the continuum β-function in the momentum space. To fix our notation, we recall the β-function to two-loop order in the SU (3) gauge theory with N f massless fermions in the fundamental representation, If β(g 2 ) has an IRFP, then β(s, g 2 ) also has a corresponding IRFP, and vice versa. In this paper, we determine β(2, g 2 ) of the SU (3) Table III for both linear and quadratic fits. In the following, we compare the results in the second column of Table III with those of Ref. [2], and Ref. [6].
First, we check the value of β(s, g 2 c ) = 0.234(30) at g 2 c = 3.51(2), which is in good agreement with the value 0.23(1) obtained in Ref. [2]. This suggests that cubic-spline interpolation can work well in the regime where the renormalized coupling varies slowly with respect to the bare coupling 6/g 2 0 . In other words, the β-function β(s, g 2 c ) reported in Ref. [2] should be valid for 0 ≤ g 2 c ≤ 3.51. Next, we check the value of β(s, g 2 c ) = 0.318(45) at g 2 c = 5.25 (2), which is quite smaller than the value 0.43(2) reported in Ref. [2]. This implies that cubic-spline interpolation fails in the regime where the renormalized coupling varies rapidly with respect to the bare coupling 6/g 2 0 . Now the value of β(s, g 2 c ) at g 2 c = 5.25 (2) is compatible with the result of a recent study of the SU (3) gauge theory with N f = 10 massless domain-wall fermions [6]. This suggests that the discrepancy between the results of Ref. [2] and Ref. [6] [2]. This confirms that using interpolation would give unreliable results for g 2 (L, a) and g 2 (sL, a), especially in the regime where they vary rapidly with respect to the bare coupling 6/g 2 0 , and consequently yield an incorrect β(s, a/L, g 2 c ) as well as the extrapolated β(s, g 2 ) in the continuum limit. Nevertheless, the resulting β(s, g 2 c ) seems to be able to capture some salient features of the β-function, e.g., the increasing/decreasing trend of β(s, g 2 c ) with respect to g 2 c , even though it cannot give the precise shape of the entire β-function in the (g 2 c , β) plane. Finally, we note that as g 2 c is increased from 5.25(2) to 6.86(2), β(s, g 2 c ) increases from 0.318(45) to 0.371(32). This implies that the slope of β(s, g 2 c ) is positive for , where β(s, g 2 c ) reaches the local maximum at g 2 max . Then for g 2 c > g 2 max , the slope of β(s, g 2 c ) becomes negative, and β(s, g 2 c ) decreases to 0.371(32) at g 2 c = 6.86 (2). To determine the exact location of g 2 max as well as the precise shape of β(s, g 2 c ) in the vicinity g 2 max is very challenging, since it requires many targeted values of g 2 c , and also we cannot rely on the renormalized couplings from interpolation, especially in this regime where the slope of β(s, g 2 c ) changes sign (from positive to negative).  Table III).
In Fig. 3, we plot β(s, g 2 c ) versus g 2 c , for g 2 c = {6.86(2), 6.92(3), 7.03(2), 7.16(2)}. In Fig. 3 (a), the data points are obtained by continuum extrapolation with the linear fit, as listed in the second column of Table III, while in Fig. 3 (b), the data points are obtained by continuum extrapolation with the quadratic fit, as listed in the fourth column of Table   III. In both cases, the data points are well fitted by the linear approximation of β(s, g 2 ), β(s, g 2 ) = dβ(s, g 2 ) dg 2 In Fig. 3 (a), the linear fit gives with χ 2 /d.o.f. = 0.06, while in Fig. 3 (b), the linear fit gives with χ 2 /d.o.f. = 0.16. Note that our convention for β(s, g 2 ) in Eq. (7) is the negative of the conventional β-function in the continuum (8) and thus gives a negative slope β (1) s at the IRFP, as shown in Fig. 3. We omit the negative sign in Eqs. (11) and (13) to conform with the conventional β-function in the continuum.
These two sets of results (10)- (13) are consistent with each other within error bars, which seems to imply the existence of an IRFP at g 2 * ∈ [7.20, 8.03]. However, we have not measured β(s, g 2 c ) for g 2 c > 7.16. Thus it is uncertain whether β(s, g 2 c ) would behave like Eq. (9) all the way from g 2 c = 7.16 to g 2 * , or if it would start to bounce back at some point g 2 min > 7. 16 and become an increasing function of g 2 c for g 2 c > g 2 min . The former scenario implies that the theory is infrared conformal with the fixed point at g 2 * ∈ [7.20, 8.03], while the latter suggests that the theory is near-conformal, depending on how closely β(s, g 2 c ) approaches zero.
Note that in the limit s → 1, Eq. (15) gives The significance of Eqs. (15) and (16) is that the slope of β(s, g 2 ) at the IRFP (with s 1) can be used to determine that at s = 1, i.e., the slope of the β-function β(g 2 ) at the IRFP, which is equal to γ * g /2. Substituting Eq. (11) into Eq. (15) gives while putting Eq. (13) into Eq. (15) gives These two results are consistent with each other within error bars. They are also compatible with the results in the weak-coupling perturbative theory, 0.473 (the schemeindependent value to the fifth order) and 0.853 (to four-loop order in the MS scheme), as given in Ref. [27].

VII. DISCUSSION AND CONCLUSION
In this paper, I performed an improved study of the β-function of SU (3) gauge theory with N f = 10 massless optimal domain-wall fermions in the fundamental repre- In the first scenario, β(s, g 2 c ) would behave like Eq. (9) all the way from g 2 c = 7.16 to g 2 * , and the theory is infrared conformal. Combining the fitting results from Figs. 3(a) and 3(b) gives g 2 * = 7.55 ± 0.36, and the universal scaling exponent of β(g 2 ), γ * g = 0.88 ± 0.40. In the second scenario, β(s, g 2 c ) would behave like a decreasing function of g 2 c for g 2 c > 7.16 until it reaches the local minimum at g 2 min , when it bounces back and becomes an increasing function of g 2 c for g 2 c > g 2 min . The question is how closely the minimum β(s, g 2 min ) approaches zero. To investigate whether the theory is near-conformal or conformal for g 2 c > 7.16 requires much more computing resources than that was available to this study. Note that the HMC simulations become more expensive as g 2 c becomes larger (or, equivalently, 6/g 2 0 becomes smaller).
Recently, a study of the β-function in the SU (3) lattice gauge theory with N f = 10 massless staggered fermions in the fundamental representation was presented in Ref. [7], with a preview in Ref. [28]. The continuum β-function β(s, g 2 c ) in Ref. [7] is a monotonic increasing function of g 2 c ∈ [5.0, 7.7], in complete disagreement with the four data points of β(s, g 2 c ) in Fig. 3. Such a dramatic discrepancy looks rather striking. In the following, I compare the results of Ref. [7] at g 2 c = 7.0 with those in this study at g 2 c = 7.03 (2). In Ref. [7], the step-scaling β-function β(s, a/L, g 2 c ) was obtained with five lattice pairs (sL, L)/a = { (24,12), (32, 16), (36,18), (40,20), (48, 24)}, which is a monotonic decreasing function of (a/L) 2 , for g 2 c = 7.0. This is completely different from the β(s, a/L, g 2 c ) in this paper, which is a monotonic increasing function of (a/L) 2 , as shown in the top-right panel of Fig. 2 for g 2 c = 7.03 (2). Consequently, the continuum β-function in Ref. [7] became very large, β(s, g 2 c ) = 0.75(4) at g 2 c = 7.0, which is completely different from the β(s, g 2 c ) = 0.299(35) at g 2 c = 7.03 (2) in this paper (see Table III). What would cause such a dramatic discrepancy between these two studies of the β-function of SU (3) lattice gauge theory with N f = 10 massless lattice fermions ?
First, could it be due to the residual mass at finite N s = 16 in this study ? As shown in Table I, the residual masses are all very tiny and quite uniform across all lattice sizes and couplings. Even if the residual mass has some additive correction to the renormalized coupling, say, g 2 (L, a) → g 2 (L, a) + δ(m res a), it would be canceled in the step-scaling βfunction [g 2 (sL, a) − g 2 (L, a)]/ ln(s 2 ), since (m res a) sL ∼ (m res a) L . Thus the residual mass has almost no effect on either the step-scaling function β(s, a/L, g 2 c ) or its value in the continuum limit. So we rule out the possibility that the residual mass could change the slope of the step-scaling function β(s, a/L, g 2 c ) at g 2 c = 7.03(2) (see the top-right panel of Fig. 2) from positive to negative. Next, does the residual mass have any effect on the shape/location of the continuum β-function in the (g 2 c , β) plane ? Now, for g 2 c itself, g 2 c → g 2 c + δ(m res a) without cancellation. Since (m res a) L is almost constant (with fluctuations less than 20%), for all g 2 c on the primary lattices (L/a = 8, 10, 12, 16) it gives δ(m res a) ∼ δ for all g 2 c , and the curve of β(s, g 2 c ) in the (g 2 c , β) plane is shifted to β(s, g 2 c + δ) with almost no change in its shape. If the theory is infrared conformal, the IRFP is shifted from g 2 * to g 2 * + δ, while the slope β (1) s of β(s, g 2 c ) at the IRFP and γ * g are not affected. In view of the tiny residual masses in Table I, I suspect that δ is already much smaller than the error of g 2 c resulting from tuning g 2 (L, a) = g 2 c for all primary lattices. In general, for any study with DWFs, if δ(m res a) is a monotonically increasing function of g 2 (L, a), then the shape of the curve β(s, g 2 c ) would be a little bit stretched along the positive direction of the g 2 c axis, due to the non-uniformity of δ(m res a). If the theory is infrared conformal, the measured location of the IRFP would be a little larger than the exact g 2 * (at zero residual mass), and also the measured slope of the β-function at the IRFP would be smaller than its exact β (1) s . Consequently, the measured universal scaling exponent would be a little smaller than the exact γ * g (at zero residual mass). Likewise, if the theory is infrared nearconformal, the measured g 2 min would be a little larger than the exact g 2 min (at zero residual mass). From the above discussions, the effect of the residual mass in this study should be very small. Thus it is impossible to change the slope/curvature of β(s, g 2 c ) in Fig. 3 from negative to positive. So we rule out the possibility that the residual mass could produce such a dramatic discrepancy in β(s, g 2 c ) at g 2 c ∼ 7.0, namely, 0.299(35) in Table III versus 0.75(4) in Ref. [7].
Next, could this be due to the volumes being too small in this study ? Would it be possible to make a dramatic change in the continuum extrapolation if we include a larger volume, say, 48 4 in our analysis ? From the data for β(s, a/L, g 2 c ) at g 2 c = 7.03(2), as shown in the top-right panel of Fig. 2, the rate of change of β(s, a/L, g 2 c ) with respect to (a/L) 2 is rather small at any (a/L) 2 . Even if we add a larger volume, say 48 4 , with an additional data point of β(s, a/L, g 2 c ) at (a/L) 2 = (1/24) 2 ∼ 0.00174, it is very unlikely that the slope of β(s, a/L, g 2 c ) would undergo a dramatic change from a small positive slope to a large negative slope in the limit (a/L) → 0. Note that the (a/L) 4 correction gets smaller for larger L as (a/L) → 0. Consequently, the deviation of the step-scaling β-function β(s, a/L, g 2 c ) from the linear function of (a/L) 2 gets smaller as L gets larger. In other words, in this study, adding an extra data point with a larger volume for the step-scaling function β(s, a/L, g 2 c ) at g 2 c = 7.03(2) is very unlikely to make a dramatic change to its value in the continuum limit (a/L → 0). Note that in this study only one data point of β(s, a/L, g 2 c ) at the largest g 2 c = 7.16 (2) and at the smallest volume with (a/L) 2 = (1/8) 2 ∼ 0.016 has a noticeable correction from the (a/L) 4 Table III. Thus we rule out the possibility that adding data points of β(s, a/L, g 2 c ) with larger volumes in this study could produce such a dramatic difference in β(s, g 2 c ) at g 2 c ∼ 7.0, namely, ∼ 0.3 in Table III versus ∼ 0.75 in Ref. [7].
Finally, we compare the actions in this study with those in Ref. [7]. The gauge action in Ref. [7] is the tree-level improved Symanzik gauge action, which is different from the Wilson plaquette action in this study. However, we do not expect that different gauge actions would cause such dramatic differences in any observables. Then we come to the possibility that the dramatic discrepancies are due to two different lattice fermion actions. If both lattice fermion Dirac operators belong to the same universality class of the continuum Dirac operator, then they should produce consistent results in the continuum limit. Could the staggered fermion operator violate fermion universality in the vicinity of the IRFP ? This conjecture has been addressed by the authors of Ref. [6]; however, it was refuted by the authors of Ref. [28]. A nonperturbative analytic proof seems to be required to settle the issue of whether the (rooted) staggered fermions belong to the same universality class of the continuum Dirac operator, especially in the vicinity of the IRFP. At the moment, the results of this study could not rule out those in Ref. [7], and vice versa. Moreover, I do not see any other (systematic/statistical) possibilities that can reconcile the dramatic discrepancies between these two studies of the β-function of the SU (3) lattice gauge theory with N f = 10 massless lattice fermions.
To conclude, based on the four data points of β(s, g 2 c ) as shown in Fig. 3, I infer that the theory is infrared near-conformal, or conformal with the fixed-point g 2 * = 7.55(36). This also implies that the SU (3) gauge theory with N f = 12 massless fermions in the fundamental representation is most likely infrared conformal with IRFP g 2 * < 7.2. This prediction is consistent with a recent study with N f = 12 domain-wall fermions [29], which suggests that the theory is infrared conformal with an IRFP g 2 * ∼ 6.