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

We present a new lattice study of the discrete beta function for SU(3) gauge theory with Nf=8 massless flavors of fermions in the fundamental representation. Using the gradient flow running coupling, and comparing two different nHYP-smeared staggered lattice actions, we calculate the 8-flavor step-scaling function at significantly stronger couplings than were previously accessible. Our continuum-extrapolated results for the discrete beta function show no sign of an IR fixed point up to couplings of g^2~14. At the same time, we find that the gradient flow coupling runs much more slowly than predicted by two-loop perturbation theory, reinforcing previous indications that the 8-flavor system possesses nontrivial strongly coupled IR dynamics with relevance to BSM phenomenology.

1 Introduction and review of previous work SU(3) gauge theory with N f = 8 flavors of massless fermions in the fundamental representation is interesting both theoretically and in the context of phenomenology for new physics beyond the standard model (BSM). Theoretical interest comes from the possibility that N f = 8 may be close to the lower boundary of the conformal window, the range of N (c) f ≤ N f < 16.5 for which the theory flows to a chirally symmetric conformal fixed point in the infrared (IRFP) [1,2]. The connection to BSM phenomenology stems from expectations that mass-deformed models with N f near N (c) f will possess strongly coupled approximately conformal dynamics, producing a large mass anomalous dimension and slowly running ("walking") gauge coupling across a wide range of energy scales [3][4][5]. In models of new strong dynamics, these features are invoked to evade constraints from flavor-changing neutral currents, to produce a phenomenologically viable electroweak S parameter, and to justify a relatively light and SM-like composite Higgs boson with mass M H = 125 GeV. See refs. [6,7] for brief reviews of these issues.
The onset of IR conformality is an inherently nonperturbative phenomenon: at the two-loop perturbative level, the conformal window opens with the appearance of an IR fixed point in the infinite-coupling limit. This occurs at a non-integer N (c) f ≈ 8.05 very close to N f = 8. Although both three-and four-loop perturbative calculations of the renormalization group β function in the MS scheme predict an 8-flavor IRFP, the resulting fixed-point coupling is still quite strong, g 2 MS ≈ 18.4 and 19.5, respectively [8]. There is no reason to trust perturbation theory at such strong couplings. We also do not wish to rely on arguments that spontaneous chiral symmetry breaking should be induced for g 2 MS ∼ 10, which combine perturbation theory with an approximate analysis of Schwinger-Dyson equations [9]. The resulting estimates of N f ≈ 8 in ref. [10] to N (c) f ≈ 12 in ref. [11], while a bound N (c) f 12 follows from a conjectured thermal inequality [12].
Since interest in 8-flavor SU(3) gauge theory revolves around its strongly coupled IR dynamics, lattice gauge theory is an indispensable approach to study the system nonperturbatively, from first principles. A wide variety of methods have been employed by existing lattice studies. These include investigation of the running coupling and its discrete β function [13,14]; exploration of the phase diagram through calculations at finite temperature [15][16][17][18][19]; analysis of hadron masses and decay constants [20][21][22][23][24]; and study of the eigenmodes of the Dirac operator [20,25]. These various analyses are complementary, and in combination offer the most reliable information about the IR dynamics of the system. Let us summarize the strengths of each of these approaches, and review the current state of knowledge for the 8-flavor system, to motivate the new work that we will present.
In this paper we will report on a new step-scaling study of the 8-flavor discrete β function, exploiting several recent improvements to this method. Generically, running coupling studies are carried out in the am = 0 chiral limit, and connect the perturbative (asymptotically free) UV regime to the strongly coupled IR. The IR limit of a massless theory is characterized by either spontaneous chiral symmetry breaking or renormalization group flow to an IR fixed point. Lattice running coupling studies, after extrapolation to the continuum, directly search for an IRFP within the range of renormalized couplings probed by the study. At the same time, the use of massless fermions prevents these studies from exploring chirally broken dynamics, which finite-temperature or spectral techniques are better suited to investigate.
For example, the (pseudo)critical couplings g cr of chiral transitions at finite temperature T and nonzero fermion mass am depend on the lattice spacing a, or equivalently on the temporal extent of the lattice N t = 1/(aT ). In a chirally broken system, these transitions must move to the asymptotically free UV fixed point g cr → 0 as the UV cutoff a −1 → ∞. At the same time the fermion mass must be extrapolated to the am → 0 chiral limit to ensure that the observed chiral symmetry breaking is truly spontaneous. In an IR-conformal system, in contrast, the finite-temperature transitions in the chiral limit must accumulate at a finite coupling as N t → ∞, and remain separated from the weak-coupling conformal phase by a bulk transition. Spectral studies can proceed more directly by attempting to fit nonzero-mass lattice data to chiral perturbation theory. Since the chiral regime is inaccessible to existing studies, these investigations typically search for simpler signs that the pseudoscalar mesons behave as Goldstone bosons in the chiral limit, for instance by considering whether the ratio of vector and pseudoscalar meson masses M V /M P → ∞ as am → 0. In a similar vein, eigenmode studies can investigate chiral symmetry breaking by comparing the low-lying Dirac spectrum with random matrix theory, or by considering the scale dependence of the effective mass anomalous dimension predicted by the eigenmode number.
Spectral and eigenmode studies have further applications beyond simply searching for spontaneous chiral symmetry breaking. The hadron masses themselves are phenomenologically interesting. In addition to exploring whether the system may possess a sufficiently light Higgs particle, these calculations predict the properties of further resonances that may be observed at the Large Hadron Collider or future experiments. The low-energy constants of the effective chiral Lagrangian are also experimentally accessible, for example in the form of the electroweak S parameter and WW scattering lengths [26,27]. Finally, in approximately conformal systems, finite-size scaling of spectral data can probe the effective mass anomalous dimension γ eff (µ), the scale dependence of which can be extracted from the Dirac eigenmodes [25,28].
In the context of the 8-flavor system, a pioneering lattice investigation performed a running coupling study based on the Schrödinger functional [13,14]. This work could access the continuum-extrapolated discrete β function up to g 2 SF 6.6, in which range reasonable agreement with two-loop perturbation theory was found. In part, computational expense limited the strength of the renormalized coupling that could be considered. In addition, the study had to avoid a bulk phase transition at stronger bare couplings, a typical restriction that prevents lattice calculations from probing arbitrarily strong couplings. Since refs. [13,14] used unimproved staggered fermions, one may expect to reach stronger couplings and to reduce computational costs by improving the lattice action, which is one of the steps we take in the present work.
Given the evidence from refs. [13,14] for rough consistency with perturbation theory up to g 2 SF ≈ 6.6, we can turn to finite-temperature and spectral studies to explore whether chiral symmetry is spontaneously broken at these couplings. The pioneering 8-flavor finitetemperature study of ref. [15], later extended by ref. [16], investigated N t = 6, 8 and 12 with fixed am = 0.02, 1 for which mass the chiral transitions move to weaker coupling for larger N t in agreement with two-loop perturbation theory. In order to explore the approach to the chiral limit, in recent work we carried out finite-temperature investigations for a range of fermion masses am ≤ 0.02 with N t = 12, 16 and 20 [18,19]. For sufficiently large masses am ≥ 0.01 we also observed two-loop scaling, but this did not persist at smaller am ≤ 0.005, where the finite-temperature transitions merged with a bulk transition into a lattice phase. (We will review our lattice phase diagram in section 3.) Even ongoing studies using a rather large 48 3 ×24 lattice volume, part of a joint project with the Lattice Strong Dynamics Collaboration, have not yet established spontaneous chiral symmetry breaking, as we will report in a future publication [29].
Similarly, studies of the 8-flavor spectrum and Dirac eigenmodes have not clearly demonstrated spontaneous chiral symmetry breaking. In ref. [21] the LatKMI Collaboration argued that at lighter fermion masses 0.015 ≤ am ≤ 0.04 the spectrum of the theory may be described by chiral perturbation theory, while data at heavier 0.05 ≤ am ≤ 0.16 appear to exhibit some remnant of IR conformality despite chiral symmetry breaking. At smaller masses 0.004 ≤ am ≤ 0.01 and larger lattice volumes up to 48 3 ×96, however, a US-BSM project could not confirm spontaneous chiral symmetry breaking [23]. Recent work by the Lattice Strong Dynamics Collaboration using the domain wall fermion formulation (as opposed to the staggered fermions used by all other studies discussed above) observed a slight but steady increase in the ratio M V /M P for smaller fermion masses in the range 0.0127 ≤ am ≤ 0.0327, even though their data were not within the radius of convergence of chiral perturbation theory.
In summary, although existing lattice studies are all consistent with 8-flavor SU(3) gauge theory being chirally broken, with no evidence for IR conformality, spontaneous chiral symmetry breaking has not yet been conclusively established. The implications of this situation extend well beyond a simple categorization of the system. In particular, the lattice results provide indications that the N f = 8 model exhibits the desirable phenomenological features expected for N f ≈ N (c) f . When analyzed within the framework of mass-deformed IR conformality, the spectral studies mentioned above prefer a large effective mass anomalous dimension γ eff 1. Our investigations of Dirac eigenmode scaling find that this large γ eff (µ) persists across a wide range of energy scales [25]. Arguably the most exciting recent development is the observation of a light flavor-singlet scalar Higgs particle by the LatKMI Collaboration [22]. 2 From these considerations, we conclude that further lattice studies of the 8-flavor system are well motivated. In this paper we present a new study of the discrete β function, taking two novel steps in order to access stronger couplings than were previously probed for N f = 8. First, instead of using the traditional Schrödinger functional running coupling discussed above, we employ a recently introduced alternative based on the gradient flow, which offers improved statistical precision for lower computational costs. We review gradient flow step scaling in the next section, also summarizing several recent improvements that make this method more robust against systematic errors. In addition, we make use of highly improved lattice actions, comparing two staggered-fermion actions with either one or two nHYP smearing steps. (The once-smeared action is also being used in separate finite-temperature [29], spectral [23] and eigenmode [25] studies, which offer complementary insight into additional aspects of this system.) In section 3 we describe our numerical setup and lattice ensembles, focusing on the issue of how to reach strong renormalized couplings without encountering a bulk transition into a lattice phase.
Our step-scaling analyses and results are presented in section 4. Our nonperturbative study predicts the continuum-extrapolated discrete β function of 8-flavor SU(3) gauge theory up to renormalized couplings g 2 c ≈ 14. For much of this range we find that the coupling runs much more slowly than in two-loop perturbation theory, and also more slowly than the (IR-conformal) four-loop MS prediction. We conclude in section 5 with discussion of important directions to pursue in further future studies of eight flavors on the lattice.

Gradient flow step scaling and its improvement
The gradient flow is a continuous transformation that smooths lattice gauge fields to systematically remove short-distance lattice cutoff effects [33]. Following the demonstration that the gradient flow is mathematically well defined and invertible [34], it has been used in a wide variety of applications (recently reviewed by ref. [35]). We are interested in step-scaling studies of a renormalized coupling defined through the gradient flow. This coupling is based on the energy density E(t) = − 1 2 ReTr [G µν (t)G µν (t)] after flow time t, which defines [36] The normalization N is set by requiring that g 2 GF (µ) agrees with the continuum MS coupling at tree level. To use the gradient flow coupling in step-scaling analyses, we tie the energy scale to the lattice volume L 4 by fixing the ratio c = √ 8t/L, as proposed by Refs. [37][38][39]. Each choice of c defines a different renormalization scheme, producing a different renormalized coupling g 2 c (L) and predicting a different discrete β function in the continuum limit. If periodic boundary conditions (BCs) are used for the gauge fields, these β functions are only one-loop (and not two-loop) universal [37].
At nonzero bare coupling g 2 0 , the gradient flow renormalized couplings g 2 c have cutoff effects that must be removed by extrapolating to the (a/L) → 0 continuum limit. The cutoff effects depend on the lattice action used to generate the configurations, on the gauge action used in the gradient flow transformation, and on the lattice operator used to define the energy density E(t). While it is possible to systematically remove lattice artifacts by improving all three quantities simultaneously, this approach is not always reasonable in practice. Another option proposed by ref. [40] is to modify the definition of the renormalized coupling to perturbatively correct for cutoff effects, Here the function C(L, c) is a four-dimensional finite-volume sum in lattice perturbation theory, which depends on the action, flow and operator. It is computed at tree level by ref. [40], and we use that result to include this correction in our definition of g 2 c . Since we use periodic BCs for the gauge fields, the correction C(L, c) also includes a term that accounts for the zero-mode contributions. Even with this tree-level improvement, the gradient flow step scaling can show significant cutoff effects. These can be reduced to some extent by working with relatively large c 0.3, at the price of increased statistical uncertainties [39]. In ref. [32] we introduced a different modification of the renormalized coupling that replaces the energy density E(t) with the value resulting from a small shift in the flow time, with |τ 0 | t/a 2 . This t-shift τ 0 can be either positive or negative. In the continuum limit τ 0 a 2 → 0 and g 2 GF (µ) = g 2 GF (µ). For O(a)-improved actions like those we use, a simple calculation shows that it is possible to choose an optimal τ 0 value τ opt such that the t-shift removes the O(a 2 ) corrections of the coupling g 2 GF (µ; a) defined in eq. 2.3. In our previous studies of both the 4-and 12-flavor SU(3) systems [32], this τ opt depended only weakly on g 2 GF (µ), and simply setting it to a constant value sufficed to remove most observable lattice artifacts throughout the ranges of couplings we explored in each case.
Since the gradient flow is evaluated through numerical integration, replacing g 2 c → g 2 c by shifting t → t+τ 0 does not require any additional computation. The t-shift also does not interfere with the perturbative correction in eq. 2.2, and in the following we will combine both improvements, searching for the optimal τ opt after applying the tree-level perturbative corrections. Using the resulting g 2 c gradient flow running coupling, we will investigate the 8-flavor discrete β function corresponding to scale change s, .
This quantity is sometimes called the step-scaling function σ s (u, L) with u ≡ g 2 c (L; a), and we will use these terms interchangeably. Our final results for the continuum discrete β function β s ( g 2 c ) = lim (a/L)→0 β s ( g 2 c , L) are then obtained by extrapolating (a/L) → 0. We emphasize that different values of τ 0 should all produce the same β s ( g 2 c ) in the continuum limit [32]. In section 4 we will see that this is not actually the case for one of the lattice actions we consider. With two nHYP smearing steps the continuum extrapolations with different t-shifts disagree by statistically significant amounts. We will account for these discrepancies as one source of systematic uncertainty.

Numerical setup and lattice ensembles
We carry out numerical calculations using nHYP-smeared staggered fermions with smearing parameters α = (0.5, 0.5, 0.4) and either one or two smearing steps. The gauge action includes fundamental and adjoint plaquette terms with couplings related by β A /β F = −0.25. We keep the fermions exactly massless, which freezes the topological charge at Q = 0. We impose anti-periodic BCs for the fermions in all four directions, but the gauge fields are periodic.
Previous studies of this lattice action with one nHYP smearing step observed an " S 4 " lattice phase in which the single-site shift symmetry (S 4 ) of the staggered action is spontaneously broken [18,19,41]. In the massless limit, a first-order transition into the S 4 phase occurs at β (c) F ≈ 4.6. The twice-smeared action also has an S 4 phase that is separated from the weak-coupling phase around β (c) F ≈ 3.6. In this work we consider only weaker couplings safely distant from the S 4 lattice phase. Although the bare couplings β F for these two different lattice actions are not directly comparable, we find that two smearing steps do allow us to access stronger renormalized couplings before encountering the S 4 phase (cf. figure 1). This is consistent with our expectations; the possibility of probing stronger couplings was our main motivation for investigating the twice-smeared action in addition to the once-smeared case. Another benefit of considering two lattice actions is that we obtain two independent sets of results. In the continuum limit both analyses should predict the same discrete β function, so by comparing our final results from the two different actions we can check for systematic errors.
Using each action, we generate ensembles of gauge configurations with six different  smearing step we study twelve couplings in the range 5 ≤ β F ≤ 11; with two smearing steps we study nine couplings in the range 4.75 ≤ β F ≤ 7. The resulting ensembles (72 with one smearing step and 54 with two) are summarized in Tables 1 and 2 in the appendix, respectively. In figure 1 we show the gradient flow renormalized coupling g 2 c (L) measured on each ensemble for c = 0.25. These data use the optimal t-shift values τ opt determined in the next section, and also include the tree-level perturbative correction factor C(L, c) in eq. 2.2. The perturbative corrections are fairly mild for the plaquette gauge action we use for both lattice generation and gradient flow, and the clover operator we use to define the energy density. Although we do reach stronger renormalized couplings with two smearing steps, the gain is fairly modest, only ∼15% with τ 0 = 0 and less after t-shift improvement. As shown in figure 1, however, a good deal of freedom remains to extend the twice-smeared runs to stronger couplings before encountering the S 4 phase, which is located at the left edge of each plot. The computational cost of such runs prevents us from including them in the present work. As tabulated in Tables 1 and 2, the gauge fields generated in the twice-smeared runs are already quite rough, with average plaquettes approaching 1/3. In addition, as we will show below our twice-smeared results already exhibit cutoff effects significantly larger than those we observe with one smearing step, suggesting that pushing this action to stronger couplings may not be worth the computational expense.

Step-scaling analyses and results
Following the standard procedure for lattice step-scaling analyses, we will first fit our input data to some interpolating function to determine the finite-volume discrete β functions β s ( g 2 c , L) with fixed L (eq. 2.4), and then extrapolate these to the (a/L) 2 → 0 continuum limit. Because we consider the same input bare couplings β F (giving the same lattice spacings a) on every lattice volume, we can either interpolate the renormalized couplings g 2 c (L) as functions of β F , or at each input β F we can compute β s ( g 2 c , L) directly from eq. 2.4 and interpolate these as functions of g 2 c (L). We will carry out analyses using both approaches, and interpret any disagreement between them as a systematic error. A similar procedure was used by ref. [37]. In this work, we find that our results from the two approaches always agree within statistical uncertainties.
Fitting the renormalized coupling g 2 c (L) on each lattice volume to some interpolating function in the bare coupling g 2 0 ≡ 12/β F is the more traditional approach. While the choice of interpolating function is essentially arbitrary, typically some functional form motivated by lattice perturbation theory is used. For example, refs. [14,37] both fit 1 g 2 − 1 g 2 0 to polynomials in g 2 0 . Inspired by ref. [42], we instead consider rational function interpolations. Specifically, the interpolating curves shown in figure 1 use the "(2, 2)" rational function which reduces to the expected g 2 c ∝ g 2 0 at weak coupling. Most of the fits shown are of good quality, with 0.2 χ 2 /d.o.f. 1.6, corresponding to confidence levels 0.94 CL 0.16. The main outlier is the twice-smeared L = 12 interpolation, which has χ 2 /d.o.f. ≈ 3.8 and CL ≈ 0.004. While the quality of fits can be improved by adjusting the number of terms in the rational function, the final results are unchanged within statistical uncertainties.
When we interpolate the finite-volume β s (L) from eq. 2.4 as functions of g 2 c (L), it is reasonable to use the same sort of polynomial interpolating function that perturbation theory predicts for the continuum β function, In every case the unshifted τ 0 = 0 results show significant dependence on (a/L) 2 , despite the tree-level perturbative correction discussed in section 2. We wish to optimize τ 0 by finding the value τ opt for which these cutoff effects are minimized. As discussed in section 2, we will use constant τ opt for all g 2 c , which will only reduce and not completely remove O(a 2 ) effects. With one nHYP smearing step, our choice τ opt = 0.07 is satisfactory for all couplings we consider. Figure 3 shows the resulting removal of cutoff effects for one particular u = 10, while the left panel of figure 2 considers the full range of g 2 c with c = 0.25. The three curves in figure 2 are the finite-volume discrete β functions that we extrapolate to the continuum, and with τ opt = 0.07 they nearly overlap for all couplings.
The twice-smeared action is quite different. In this case we choose τ opt = 0.18, more than 2.5 times larger than the once-smeared τ opt = 0.07, which already indicates more severe cutoff effects. While this τ opt = 0.18 produces the desired nearly constant continuum   extrapolations for u = 10 in figure 4, from the right panel of figure 2 we can see that cutoff effects remain for both smaller and larger couplings. Specifically, a smaller t-shift τ 0 ≈ 0.12 produces better improvement for u 8, while a larger τ 0 ≈ 0.24 is more effective for u 12.
Our choice of constant τ opt = 0.18 is a compromise that "over-improves" at small g 2 c and "under-improves" at large g 2 c . While it is possible to use a u-dependent τ opt , we prefer to keep this improvement as simple as possible, to avoid the risk of losing predictivity by introducing too many optimization parameters [32]. Although we end up with non-trivial continuum extrapolations as indicated by the right panel of figure 2, they remain reliably linear in (a/L) 2 .
More problematically, the larger lattice artifacts of the twice-smeared action affect even the continuum-extrapolated discrete β function results. Reliable continuum extrapolations should behave as shown for the once-smeared action in figure 3, where different values of τ 0 predict the same (a/L) 2 → 0 limit β s ( g 2 c ) well within statistical uncertainties. In this case the t-shift improvement simply stabilizes the extrapolations by removing cutoff effects, without changing the continuum results. The contrast with figure 4 for two nHYP smearing In each plot we include once-and twice-smeared results with τ 0 = 0 as well as with the optimal τ opt = 0.07 and 0.18, respectively. While τ 0 optimization can change the twicesmeared continuum limit, removing a source of systematic error, the once-smeared results always agree within uncertainties. steps is dramatic, especially for smaller c. In this case different t-shifts produce continuumextrapolated β s ( g 2 c ) that disagree by statistically significant amounts. These discrepancies must be included among our systematic uncertainties, as we will now discuss.
We account for three potential sources of systematic errors: Optimization: To determine how we should account for any sensitivity to the t-shift improvement parameter τ 0 , consider figure 5. Each panel in this figure compares onceand twice-smeared continuum-extrapolated results with both τ 0 = 0 and the optimal τ opt , which should all predict the same β s ( g 2 c ). While the once-smeared results always agree within uncertainties, optimizing τ 0 produces a statistically significant change with two nHYP smearing steps, just as in figure 4. In fact, the t-shift brings the twice-smeared results into better agreement with the once-smeared action, removing systematic errors that would be present for an unimproved analysis with τ 0 = 0. The only remaining systematic uncertainties from optimization therefore result from our restriction to constant τ opt . As discussed above, τ opt = 0.07 is satisfactory for all g 2 c , so these systematic uncertainties vanish for the once-smeared action. With two nHYP smearing steps, however, τ 0 = 0.12 (0.24) is preferred for small (large) g 2 c . We conservatively define as systematic errors any discrepancies between results for either of these two τ 0 compared to those for τ opt = 0.18. These systematic errors tend to be quite mild, at least 3.5 times smaller than the statistical uncertainties. Interpolation: As discussed at the start of this section, we analyze our data both by interpolating g 2 c (L) as functions of β F and by interpolating β s ( g 2 c , L) as functions of g 2 c (L). We take our final results from the latter analysis. Any discrepancies between the two approaches we include as a systematic error. For the 8-flavor analyses we carry out in this work, these systematic errors always vanish. Extrapolation: Even after accounting for tree-level perturbative corrections and t-shift improvement, our continuum extrapolations are not always perfectly linear in (a/L) 2 .
To determine the resulting systematic effects, we repeat all analyses without including the smallest-volume L = 12 → 18 data, considering only 16 → 24 and 20 → 30 points in linear (a/L) 2 → 0 extrapolations. Any discrepancies between the two-and threepoint continuum extrapolations defines our third systematic uncertainty. Although this source of systematic error also often vanishes, for some u it can be up to four times larger than the statistical uncertainty. In all three cases, we take the systematic errors to vanish when the results being compared agree within 1σ statistical uncertainties. This ensures that statistical fluctuations are not double-counted as both systematic and statistical errors. Note that to determine the systematic uncertainties from τ 0 optimization, it was important to compare multiple lattice actions. We will return to this point in section 5.
We are now ready to present our final results for the 8-flavor system. Figure 6 shows the continuum-extrapolated s = 3/2 discrete β function for two different renormalization schemes, c = 0.25 and 0.3. In both panels we include our nonperturbative results for the once-and twice-smeared actions. The darker error bands show the statistical uncertainties, while the lighter error bands indicate the total uncertainties, with statistical and systematic errors added in quadrature.
We compare our numerical results with perturbation theory, where for N f fermions transforming in representation R of the gauge group. For the fundamental representation of SU(3) gauge theory, 4) so that N f = 8 gives b 0 = 17 3 and b 1 = 2 3 . Higher-order coefficients b i are renormalization scheme dependent. In the MS scheme, ref. [8] reports numerical values b 2 ≈ −423 and b 3 ≈ 374 for 8-flavor SU(3) gauge theory. Both the three-and four-loop β functions predict an IR fixed point, but only at strong couplings g 2 MS ≈ 18.4 and 19.5 where perturbation theory is not reliable.
Along with our numerical results we include the two-and four-loop perturbative predictions for the s = 3/2 discrete β function in figure 6. The once-and twice-smeared actions predict consistent continuum results, which are significantly smaller than the twoloop perturbative curve, by more than a factor of three for g 2 c = 12. At the weakest coupling that we probe, g 2 c ≈ 2, our results are still approaching the perturbative predictions from below. (Although we mentioned in section 2 that the gradient flow discrete β function is only one-loop universal, the one-and two-loop perturbative results are almost indistinguishable across the range shown in figure 6.) Due to the large negative b 2 coefficient in the MS scheme, the four-loop discrete β function also becomes much smaller than the two-loop prediction. Even at the strongest g 2 c = 13.5 that we are able to reach with two nHYP smearing steps in the c = 0.25 scheme, our numerical results remain even smaller than four-loop perturbation theory. In the c = 0.3 scheme the twice-smeared β function becomes comparable to the maximum of the four-loop curve at the largest accessible g 2 c = 14.3. Of course, since the discrete β function is scheme dependent the c = 0.25 and 0.3 results do not have to agree. The perturbative four-loop MS β function not only corresponds to another different scheme, but is also of questionable validity at such strong couplings. Our comparisons with perturbation theory are for illustration only.

Discussion and conclusions
Before we attempt to interpret our nonperturbative results for the discrete β function in figure 6, let us review the motivations for and goals of this work. We are attracted to 8flavor SU(3) gauge theory primarily by the possibility that it may possess strongly coupled near-conformal IR dynamics, leading to desirable BSM phenomenology including a light Higgs particle [22] and large effective mass anomalous dimension across a wide range of energy scales [25]. While a variety of existing lattice studies have not yet been able to establish chiral symmetry breaking in the m → 0 limit for N f = 8, their results are all consistent with such dynamics [13][14][15][16][17][18][19][20][21][22][23][24][25]. To address this situation, we have carried out a new step-scaling study of the discrete β function, exploiting two different improved lattice actions and the recently introduced gradient flow running coupling that enabled us to investigate significantly stronger couplings than were previously accessible.
Our results in figure 6 indicate a coupling that runs much more slowly than predicted by two-loop perturbation theory, even more slowly than the four-loop MS prediction, which possesses a strongly coupled IR fixed point. Despite considering a second lattice action with two nHYP smearing steps, in addition to our usual once-smeared action, we could not reach strong enough couplings either to see a similar IRFP in our numerical results, or to obtain a clear deviation from the IR-conformal four-loop result. 3 We see no sign of spontaneous chiral symmetry breaking for running couplings as large as g 2 c ≈ 19 with τ 0 = 0 on 30 4 lattice volumes. In part, our ability to access stronger renormalized couplings is limited by the more severe lattice artifacts in our twice-smeared results. In this case, significant t-shift improvement is required to obtain agreement with the once-smeared results, illustrating the importance of cross-checking continuum predictions by comparing different lattice actions. As shown by figure 5, the large τ opt = 0.18 needed with two smearing steps significantly reduces the range of g 2 c that we can reach on lattice volumes from 12 4 to 30 4 . Even though we could still push twice-smeared computations to stronger bare couplings before encountering the S 4 lattice phase (figure 1), the severe cutoff effects we already observe suggest that doing so may not be worth the computational expense.
Future investigations of the 8-flavor system will benefit from several studies currently being carried out with the once-smeared action we considered in this work. As discussed in section 1, the Lattice Strong Dynamics Collaboration is studying the finite-temperature phase diagram with N t = 24, which still seems to be too small to establish chiral symmetry breaking in the massless limit [29]. At the same time, the lattice ensembles generated by USBSM [23] are being analyzed in search of a light scalar Higgs particle, and we have improved our techniques to extract the effective mass anomalous dimension from the Dirac eigenmode spectrum [25,28]. Although the combination of these complementary studies will shed further light on N f = 8 and its phenomenological viability as the basis of new BSM physics, our results in this work also highlight the importance of comparing studies using different lattice actions, preferably including different fermion formulations, when exploring such unfamiliar and nontrivial systems.  Table 1. Lattice ensembles with one nHYP smearing step. For each ensemble specified by the volume L 4 and gauge coupling β F , we report the total molecular dynamics time units (MDTU), the thermalization cut, and the resulting number of 100-MDTU jackknife blocks used in analyses. We also list the average plaquette (normalized to 3), to illustrate the roughness of the gauge fields.