The running coupling of 8 flavors and 3 colors

We compute the renormalized running coupling of SU(3) gauge theory coupled to N_f = 8 flavors of massless fundamental Dirac fermions. The recently proposed finite volume gradient flow scheme is used. The calculations are performed at several lattice spacings allowing for a controlled continuum extrapolation. The results for the discrete beta-function show that it is monotonic without any sign of a fixed point in the range of couplings we cover. As a cross check the continuum results are compared with the well-known perturbative continuum beta-function for small values of the renormalized coupling and perfect agreement is found.


Introduction and summary
It is well-known that SU (3) gauge theory with N f flavors of fundamental Dirac fermions has a so-called conformal window [1,2]. The window refers to the range of flavor numbers where the theory possesses an infrared fixed point, i.e. the long distance behavior is conformal but the short distance behavior is still asymptotically free. Clearly, N f = 2 is outside the window and N f = 16 is inside because the perturbative β-function has a zero at a coupling where perturbation theory is trustworthy. Unambiguously pinning down the lower end of the conformal window has been nevertheless difficult so far because only ab initio non-perturbative studies have a chance to settle the question but these same studies are plagued by systematic uncertainties. It is important to note that these systematic uncertainties are of a practical nature only and the associated error bars can in principle be reduced to arbitrary small values. In currently available results with large fermion content the systematic uncertainties were rarely, if at all, quantitatively estimated. These are however important. The statistical uncertainties can be reduced by simply increasing the statistics but after a certain point systematic uncertainties will dominate and further increasing the statistics will be pointless.
The most important of the systematic uncertainties is related to the continuum limit. In previous work we have considered N f = 4 and performed a careful continuum extrapolation; the N f = 4 theory is well-known to be outside the conformal window and the results obtained were consistent with this expectation. As the next step in approaching the lower end of the conformal window we continue our investigation with N f = 8 in the present work and pay special attention to estimating the systematic uncertainties.
The question of finding the lower end of the conformal window is an interesting field theory problem on its own, however there are reasons to be interested in the N f = 8 model for phenomenological purposes. A large class of Beyond Standard Model extensions involve a composite Higgs particle. A natural framework for a composite Higgs state is strong dynamics which solves the hierarchy problem of the Standard Model and gives rise to dynamical electro-weak symmetry breaking. One challenge (among others) that all these models should overcome is that they should contain a relatively light scalar particle, which when coupled to the Standard Model is to be identified by the 125 GeV particle found at the LHC in 2012. A light scalar might arise in strongly coupled non-abelian gauge theories with many fermions if it is not far from the conformal window. Hence the N f = 8 model might play a useful role in studying the properties of the hypothetical composite Higgs.
Past studies of the N f = 8 model include [3] using the Schroedinger functional, [4] for thermodynamics aspects and [5,6] for hadron spectroscopy.
The organization of the paper is as follows. In section 2 we briefly summarize the finite volume gradient flow running coupling scheme that we use; for more details see [7,8].
In section 3 the details of our numerical simulations are given and the collected data is presented, while in section 4 the continuum extrapolation of the discrete β-function is performed. Finally in section 5 we end with a conclusion and provide avenues for future studies.
As we were preparing our manuscript a closely related study appeared [9] using our running coupling scheme. Although the conclusions are similar to ours there are some puzzling aspects of the results which we will comment on in section 5 emphasizing the importance of controlling all the systematic uncertainties.

The gradient flow running coupling scheme
In order to investigate the infrared behavior of a model the running coupling is a natural choice. There are many well-defined schemes and one is free to choose any one of them. If the β-function in one scheme has a non-trivial zero indicating conformal behavior in the infrared, its existence is universal in every other well-defined scheme. In the current work the recently proposed finite volume gradient flow scheme [7,8] is used, which is based on Luscher's Wilson flow [10][11][12][13] related to earlier constructions by Morningstar and Peardon [14] as well as Lohmayer and Neuberger [15]. In this scheme a 1-parameter family of couplings is defined in finite 4-volume by , 1) where N corresponds to the gauge group SU (N ), t is the flow parameter, c = √ 8t/L is a constant, E(t) is the field strength squared at t and where ϑ is the 3rd Jacobi elliptic function. The numerical factors are chosen such that at leading order g 2 c = g 2 MS for all c. The gauge field is chosen to be periodic and the massless fermions are anti-periodic in all 4 directions. The coupling g c (µ) runs via the scale µ = 1/L; for more details on the gradient flow in general see [10][11][12][13] while more details on the finite volume gradient flow scheme can be found in [7,8].
A peculiar but well-known property of the femtoworld in a periodic box [16][17][18][19][20][21][22][23][24] or small volume dynamics is that the perturbative expansion of g 2 c in g MS contains both even and odd powers, g 2 c = g 2 MS (1 + O(g MS )) for N > 2. This property results in only the 1-loop β-function coefficient being the same as the well-known coefficient in MS. The case of N = 2 has been described in [8] and in the present work we focus on N = 3.
There are two considerations affecting the choice of the constant c of the 1-parameter family. If chosen too small cut-off effects will be large, if chosen too large the statistical errors will be large. In [7,8] we have found for N f = 4 that the value c = 3/10 was optimal in this sense and will use c = 3/10 also in the current work. We will however show some results at c = 1/5 in order to illustrate the c-dependence of our procedure. Henceforth the index c will nevertheless be dropped and it will be clear from the context when c = 3/10.
The expression (2.1) for the coupling was so far considered in the continuum. In a previous work [25] we have determined the lattice spacing dependence of the tree-level correction factor δ(c). The tree-level, finite volume and finite lattice spacing perturbative calculation led to the expression where p µ = 2πn µ /L with an integer non-zero 4-vector n µ and S f,g,e (p) are the tree-level expressions for the action along the flow, dynamical gauge action and the observable in momentum space, respectively and G(p) is a gauge fixing term. The finite lattice momentum sums can easily be evaluated numerically to arbitrary precision. See [25] for more details. This expression can be used to tree-level improve the coupling by simply introducing .
However, we have found that even though tree-level improvement worked very well in reducing the slope of the continuum extrapolations with N f = 4 in [25] it did not work as well in the current work with N f = 8. The reason probably lies in the fact that the larger N f is, the larger the fermion loop contributions are. And of course tree-level improvement is not sensitive to fermionic radiative corrections. We will illustrate these issues in some select cases in section 4. It should be noted that Schroedinger functional boundary conditions can also be implemented together with the gradient flow leading to a closely related scheme [26][27][28][29][30]. An advantage of the periodic boundary conditions used in the present work is that translational symmetry is not broken in the time direction. A third option is using twisted boundary conditions which was explored for SU (2) pure gauge theory to high precision in [31].

Numerical simulation
The technical details of the simulations closely follow our work on N f = 4 in [7,8]. In particular we use the staggered fermion action with 4 steps of stout improvement with ̺ = 0.12 [14]. The bare fermion mass is set to zero and anti-periodic boundary conditions in all four directions are imposed on the fermions and the gauge field is periodic. The gauge action is the tree-level improved Symanzik action [44,45]. The observable E(t) is discretized by the clover-type construction as in [11].
Along the gradient flow we use two discretizations, the Wilson plaquette action and the tree-level improved Symanzik gauge action. These setups correspond to the W SC and SSC cases in the terminology of [25]: the notation is Flow-Action-Observable and W stands for Wilson plaquette action, S for tree-level improved Symanzik action and C for the clover discretization. Both setups lead to the same continuum limit, only the size of cut-off effects is different. This fact allows for the introduction of yet another coupling definition at finite lattice spacing, which however again leads to the same continuum limit [46], Here the parameter X is arbitrary, the choice of the two coefficients, X and 1 − X, guarantees that the continuum limit of g 2 X is the same as that of g 2 SSC or g 2 W SC , i.e. the correct one. It is important to note that X is a constant and does not depend on the bare gauge coupling β or the lattice volume L/a. In practice we have found that the choice X = 1.75 is most useful. Note that in principle X could depend on the renormalized coupling but in the present work we do not explore this possibility.
Just as in [7,8] where N f = 4 was considered we do not need to take the root of the fermion determinant. Hence the results do not depend on the validity of the fourth-roottrick commonly used for QCD. The evolution along a trajectory of the hybrid Monte Carlo algorithm [47] is implemented with multiple time scales [48] and Omelyan integrator [49].
In a lattice setting the most practical method of calculating the running coupling or its β-function is via step scaling [50,51]. In this context the linear size L is increased by a factor s and the difference of couplings is defined as the discrete β-function. If the ordinary infinitesimal β-function of the theory possesses an infrared fixed point, the discrete β-function will have a zero as well. On the lattice the linear size L is easily increased to sL by simply increasing the volume in lattice units, L/a → sL/a at fixed bare gauge coupling. In the current work we set s = 3/2 and use volumes 12 4  The collected number of thermalized trajectories at each bare coupling and volume was in the range between 5000 and 10000 and every 10 th was used for measurement. The measured renormalized couplings at each β and lattice volume are shown in table 1 for the definition (3.1) using X = 1.75. By taking the difference of renormalized couplings for lattice volumes scaled by a factor s = 3/2 and at the same bare β one obtains the discrete β-function at finite lattice spacings; see figure 1. Clearly, there is no sign of a fixed point, the running is monotonically increasing, at least at finite lattice spacing, i.e. finite lattice volumes. We are of course interested in the continuum results where we turn next.

Continuum extrapolation
In order to perform a continuum extrapolation we parametrize the renormalized coupling as a function of the bare coupling, g 2 (β) at each fixed lattice volume L/a by similarly as in [52]. The order n of the polynomial may be chosen such that acceptable fits are obtained, however in this work we would like to estimate the systematic errors that come from various choices for n; see section 4.1.
Using the parametrized curves the discrete β-function (3.2) can be obtained for arbitrary g 2 (L) for fixed L/a and s = 3/2. Estimating the error on the interpolated values is straightforward because the interpolation is linear in the fit parameters C m . Then assuming that corrections are linear in a 2 /L 2 the continuum extrapolation can be performed.

Systematic error
In our previous work [7,8] the polynomial order for the interpolation (4.1) was fixed at each lattice volume. However different choices lead to similarly acceptable interpolating fits and these in turn lead to slightly different continuum results. Even though the final continuum result varies only a bit and generally within 1-σ of the statistical error in the current work we would like to estimate the systematic error as precisely as possible. In order to achieve this the histogram method introduced in [54] is used. There are two sources of systematic uncertainties. First, it is a priori unknown what interpolation function to use for the renormalized coupling as a function of β at fixed lattice volumes, and second, one may perform continuum extrapolations using 3 or 4 lattice spacings. We interpolate using (4.1) for each lattice volume, 12 4 , 16 4 , 18 4 , 20 4 , 24 4 , 30 4 , 36 4 , with three choices of polynomial orders, n = 4, 5 and 6. All together these produce 3 7 = 2187 combination of interpolations and correspondingly lead to 2187 different continuum results. Since the data on different volumes at different β are all independent we perform a Kolmogorov-Smirnov test on the 2187 interpolations and demand that only those assignments of polynomial orders are allowed to which the Kolmogorov-Smirnov test assigns at ings. Some of these are of course the same, but needs to be counted in order to have the proper weight in the final histogram. The 1233 + 813 = 2046 continuum results at each g 2 (L) can be binned in a weighted histogram and the weight can be the goodness of the fit, a weight provided by the Akaike information criterion (AIC) 1 or no weight at all. Examples of the AIC-weighted histograms are shown in figures 2-3.
Our continuum central values at each g 2 (L) are the medians of the histograms and the systematic uncertainty can then be determined by counting 68% of the total starting symmetrically from the central value. The three types of weights lead to compatible results and for our final results we use the AIC-weighted histograms.
The systematic and statistical errors are of the same order, there is never a larger factor between them than two.

Final results
At 6 chosen values of g 2 (L) the histograms of the discrete β-function for all continuum extrapolations are shown in the right panels of figures 2-3. On the left we show typical continuum extrapolations from within a 1 − σ systematic uncertainty around the median of the histograms. Clearly, all 4 lattice spacings are in the scaling region and nicely fit on a straight line with good χ 2 /dof . In fact, the choice X = 1.75 was motivated by exactly the requirement that all 4 lattice spacings should be in the scaling region. This is not a sharp requirement, one may choose any value in the approximate range 1.6 < X < 1.9.
It is quite instructive to look at the details of these figures and discuss the source of the most important systematic error, the continuum extrapolation. Our theory is a confining one in which large bare couplings (small βs) correspond to large lattice spacings. As table 1 shows large renormalized couplings are obtained with large lattice volumes and small β values. Thus, for a given renormalized coupling one reaches the continuum limit by increasing both β and the lattice volume. Since the largest volume, independently of β, was 36 4 , large renormalized couplings correspond within our parameter set to large lattice spacings and obviously large cutoff effects.
It is of obvious interest to turn this qualitative statement to a quantitative one and to determine the size of the systematic uncertainty related to this question. Most importantly, we want to know where to stop with the present lattice sizes because no controlled continuum extrapolation can be carried out any further. As our g 2 = 6 case illustrates for this large value of the renormalized coupling one has a two peak structure for the histogram. The two peaks are the result of the significant difference between using only the finer lattices with 3 points or taking 4 points (including also the coarsest lattices) for the continuum extrapolations. This phenomenon clearly indicates that the results from the coarsest lattices are starting to deviate from the a 2 scaling showed by the finer lattices. The difference between the peaks still quantifies the systematic uncertainty for g 2 = 6 and tells us that for even larger g 2 values the control over this systematic effect could be lost and finer lattices with larger lattice volumes are needed.
The discrete β-function may reliably be calculated in (continuum) perturbation theory for small values of the renormalized coupling. In terms of the well-known infinitesimal 1 and 2 loop β-function coefficients, b 1 and b 2 the discrete variant is given by As noted already in our finite volume gradient flow scheme only b 1 is the same as in every other well-defined scheme. The reason is a well-understood feature of the finite 4-volume or femtoworld [7,8]. Nevertheless we include not only the 1-loop continuum β-function but also the 2-loop approximation in our comparisons, even though strictly speaking agreement is only expected with the 1-loop result.  As mentioned in section 3 tree-level improvement [25] did not reduce the slope of the continuum extrapolations as dramatically as for N f = 4 in our previous study. The reason is presumably that the larger fermion content results in larger fermionic contributions which are, of course, completely absent from the tree-level expressions. We illustrate both points, the smaller scaling region without employing the linear combination (3.1) and the less effective tree-level improvement in figure 4. Clearly, the continuum results are always consistent, as they should be, the various choices (improvement vs. non-improvement, linear combination vs. no linear combination) only affect the slopes of the extrapolations and the size of the scaling region.
In figure 5 we illustrate another aspect mentioned in section 3, namely c-dependence. Different choices of c define different schemes, i.e. the β-function will be c-dependent. For small coupling the difference should be very small since regardless of what c is, agreement is expected with the perturbative 1-loop result. Furthermore, the expectation is that a smaller c leads to smaller statistical errors because of smaller autocorrelations and also to larger cut-off effects because of the smaller flow time t. This is illustrated in figure 5 where the continuum extrapolation is shown for g 2 = 3 and both for c = 3/10 and c = 1/5.
Finally, in figure 6 we show the continuum extrapolated β-function over the entire 0.9 < g 2 < 6.3 range accessible to our simulations together with the 1-loop and 2-loop results. The linear combination method (3.1) was used with X = 1.75 and c = 3/10 was chosen. Our non-perturbative continuum result is in nice agreement with the perturbative results for small renormalized coupling and deviates from it for larger values. Most importantly, the deviation from the perturbative 1-loop result is downward. This could have been  expected because at some higher N f value we do expect a fixed point and by continuity one might argue that this is only possible if the running is slower than the monotonically increasing 1-loop result, at least for some N f value which is not far below the conformal window. At N f = 8 we do not see a sign of a fixed point in any case, at least in the explored range 0.9 < g 2 < 6.3.

Conclusion and outlook
In this work we have continued our study of SU (3) gauge theory with many fermions. The representation was fundamental and after having examined N f = 4 in our previous work the β-function of the N f = 8 model was computed in the present work, in the continuum. The β-function does not appear to "bend back" in the coupling range we have studied hence does not support the idea that the N f = 8 model is already inside the conformal window. This result is consistent with our study of the mass spectrum which indicated spontaneous breaking of chiral symmetry at zero fermion mass [5]. The running coupling does deviate from the perturbative β-function downwards though. While preparing our manuscript the work [9] appeared. The method used there is similar to ours and the conclusions were also similar, i.e. the behavior was compatible with a monotonically increasing β-function. The explored coupling range was larger, 2.0 < g 2 < 14, but the continuum results were not compatible with perturbation theory even at g 2 = 2. This is a puzzling feature since one would expect perturbation theory to be reliable at such a small coupling. In our work we in fact show consistency with perturbation theory up to approximately g 2 = 5 and only detect deviations for larger couplings which is more in line with expectations. The reason for the discrepancy in [9] might be due to the fact that the systematic uncertainties were not adequately addressed. In our work we controlled both types, one from the a priori unknown interpolation as a function of the bare coupling at finite lattice volume and also the one coming from the continuum extrapolation.