The running coupling of the minimal sextet composite Higgs model

We compute the renormalized running coupling of SU(3) gauge theory coupled to N_f = 2 flavors of massless Dirac fermions in the 2-index-symmetric (sextet) representation. This model is of particular interest as a minimal realization of the strongly interacting composite Higgs scenario. A recently proposed finite volume gradient flow scheme is used. The calculations are performed at several lattice spacings with two different implementations of the gradient flow allowing for a controlled continuum extrapolation and particular attention is paid to estimating the systematic uncertainties. For small values of the renormalized coupling our results for the beta-function agree with perturbation theory. For moderate couplings we observe a downward deviation relative to the 2-loop beta-function but in the coupling range where the continuum extrapolation is fully under control we do not observe an infrared fixed point. The explored range includes the locations of the zero of the 3-loop and the 4-loop beta-functions in the MSbar scheme. The absence of a non-trivial zero in the beta-function in the explored range of the coupling is consistent with our earlier findings based on hadronic observables, the chiral condensate and the GMOR relation. The present work is the first to report continuum non-perturbative results for the sextet model.


Introduction
We study SU (3) gauge theory coupled to N f = 2 flavors of massless Dirac fermions in the 2-index-symmetric (sextet) representation. The model may be a minimal realization of the strongly interacting composite Higgs scenario [1][2][3][4][5][6][7][8] which is our primary motivation. This possibility has attracted active experimental interest recently [9]. Since the dynamics is intrinsically strongly coupled only a non-perturbative approach like the lattice can determine the ultimate viability of the model. Several simulation results have been reported, but none which fully controls all systematic errors and provides a consistent picture of the various aspects of the model such as the running coupling, the anomalous mass dimension, finite temperature, hadronic observables, etc. Recent reviews are given in [8,10].
Even though the currently available simulations are performed at finite lattice spacing and may be affected by other uncontrolled systematics, the evidence suggested by these results is that the infrared behavior of the model is QCD-like. Chiral symmetry is broken spontaneously [5,8,[11][12][13] leading to a thermal phase transition at finite temperature [14][15][16][17][18][19]. Results on rough lattices concerning the running coupling in the Schroedingerfunctional scheme can be found in [20,21].
The goal of the present work is to calculate the β-function of the renormalized running coupling over a wide range of couplings in such a way that all systematics are fully controlled. We are interested in the massless theory and our simulations are set up in such a way that the bare mass can be set to zero directly. The running is implemented via changing the physical volume of the system, introducing a scale µ = 1/L, hence finite volume dependence is translated into the renormalization group flow of the coupling. This leaves us with the continuum extrapolation as the only source of systematic uncertainties.
By performing the simulations at several pairs of lattice volumes in the framework of a step scaling analysis, we gain control over the finite lattice spacing effects. In fact at each value of the renormalized coupling, results at five values of the lattice spacing are used, with two different discretizations for the observables in question and a reliable continuum extrapolation is hence possible. Our work is the first to report on fully controlled continuum results for the N f = 2 sextet model.

The gradient flow running coupling scheme
The gradient flow [22][23][24][25][26][27][28] is a particularly useful tool for studying the running coupling of a non-abelian gauge theory. There are various finite volume setups which mainly differ in the choice of boundary conditions for the gauge field [29,30,[40][41][42][43][44][45]. In the present work we follow [29,30] where the gauge field is taken as periodic in all four directions. For the fermions on the other hand we impose anti-periodic boundary conditions again in all four directions. Other applications of the gradient flow can be found in [46][47][48].
More precisely, in our scheme a 1-parameter family of couplings is defined in finite 4-volume L 4 by , where t is the flow parameter, N corresponds to the gauge group SU (N ), c = √ 8t/L is a constant, E(t) is the field strength squared at t > 0 and the numerical factor δ(c) = − c 4 π 2 3 + ϑ 4 e −1/c 2 − 1 (2.2) is chosen such that at leading order g 2 c agrees with the coupling in MS for all c; ϑ is the 3rd Jacobi elliptic function. Hence the coupling g c (µ) runs via the scale µ = 1/L, for more details see [29,30].
There is a peculiarity of our running coupling scheme related to the fact that we impose periodicity on the gauge fields, leading to zero modes [49][50][51][52][53][54][55][56][57]. These gauge zero modes cause the perturbative expansion of g 2 c in g MS to contain both even and odd powers and potentially logarithms too. However for N > 2 logarithms do not appear in the first two non-trivial orders, only polynomials [30]. The first unusual odd power results in only the 1-loop β-function coefficient being the same as in MS.
The constant 0 < c ≤ 1/2 specifying the scheme can in principle be chosen at will. However, as discussed in [29,30], a small c leads to small statistical errors but large cut-off effects and a larger c results in larger statistical errors and smaller cut-off effects. In the present work we set c = 7/20, slightly higher than the value c = 3/10 in [29,58], in order to reduce cut-off effects.

Rooted staggered formulation
The fermion doublet in the staggered fermion implementation requires the square root of the fermion determinant, also known as the rooting procedure. With the mass of the fermion doublet set to zero, the continuum step β-function as determined from a scale-dependent renormalized coupling g 2 R (L) shows no sign of turning zero in the explored range, as we will see. Our results are consistent with chiral symmetry breaking when probed with finite fermion mass deformations in the p-regime [8].
Some preliminary work, with goals similar to ours but in the Wilson fermion formulation, reports consistency with a zero in the β-function in the renormalized running coupling in the same range where our β-function is positive and monotonically growing [59][60][61]. If confirmed and further supported with conformal scaling laws, the new work in Wilson formulation might suggest a conformal infrared fixed point at vanishing fermion mass in the sextet model, inconsistent with our results and spontaneously broken chiral symmetry.
Motivated by this controversy, doubts were raised about our results questioning the application of the rooting procedure with the mass of the staggered fermion doublet set to zero in the simulations at finite lattice spacing. This lead to the speculation that in the staggered rooting procedure setting the fermion mass m to zero at fixed finite lattice spacing a might be incorrect because of the non-locality of the rooted staggered action we appear to deploy by interchanging the so-called required limit of a → 0 first, while holding the fermion mass non-zero before taking the m → 0 limit in the continuum theory as the last step, after the cutoff is removed. As a consequence of this issue of non-locality concerns were raised whether a rooted staggered theory is in the correct universality class of the continuum theory and whether rooting can identify a conformal theory [59][60][61].
To alleviate the concerns, we will show that the rooting procedure is correct when the fermion mass is set to zero at finite lattice spacing while the finite physical volume of the continuum limit is held fixed. Consequently, the conformality of a model would not be missed and the rooted staggered formulation in finite physical volume and in the infinite volume limit are expected to remain in the correct universality class.

Review of rooting in infinite volume
The method to address the rooting procedure properly has been developed in a series of papers by Bernard, Golterman, Shamir, and Sharpe [31][32][33][34][35][36][37] when the renormalized fermion mass is kept finite before the continuum limit is taken. We adapt their analysis to our model in finite physical volumes to demonstrate that the rooting procedure we apply at vanishing fermion mass and finite lattice spacing a should remain valid on the level of their reasoning.
The main results of the analysis at finite fermion mass and infinite volume are summarized first from two succinct exposures of the rooting issues [35,37] that we closely follow here. Accordingly, the rooted staggered action is defined on a fine-grained lattice with lattice spacing a f and connected with infrared physics using n renormalization group steps to a blocked physically equivalent lattice action on a coarse lattice with lattice spacing a c which is held fixed on some physical scale [35,37]. At fixed a c Λ IR , where Λ IR designates some non-perturbative infrared scale, the continuum limit a f → 0 is investigated when n → ∞.
In the technical implementation of the RG procedure, the blocked and unrooted staggered Dirac operator D stag,n is split into the taste invariant part D inv,n = D n ⊗ 1 4 with exact taste symmetry of four degenerate fermions and the taste breaking part ∆ n after each blocking step, where Tr denotes the trace in taste space. The trace of ∆ n vanishes in taste space, and 1 4 designates the taste identity matrix. The local taste invariant theory represented by D n ⊗ 1 4 has four degenerate fermions in taste space and the fourth root of the fermion determinant is trivially given by After taking the fourth root of D stag,n , an estimate is needed to show the convergence of Det 1/4 (D stag,n ) to the local single taste determinant Det(D n ) in the n → ∞ limit. As shown in Eq. (3.5), the convergence of this expansion is controlled by D −1 inv,n ∆ n . We need to estimate a c ∆ n and (a c D inv,n ) −1 separately on the scale of the coarse lattice.
Following [35] we can safely assume the bound a c ∆ n < ∼ a f /a c to hold on the coarse lattice and scaling with a f . This is a basic feature of unrooted fermions simply stating that taste-breaking disappears in the continuum limit. By exploiting the proximity of the local re-weighted theory defined with D inv,n after a large number n of blocking steps, it was argued that the bound a c ∆ n < ∼ a f /a c and its scaling with a f /a c is also valid in rooted theories [35]. We will adapt this argument. The estimate for the upper bound on the inverse of the taste invariant operator is given by (a c D inv,n ) −1 ≤ 1/(a c m R (a c )) with the renormalized fermion mass set at the physical scale a c . The important combined estimate follows with D −1 inv,n ∆ n ≤ a f /(a 2 c m R (a c )) (3.3) and the small expansion parameter where in the first step the lattice spacing is doubled by the change from staggered fermion basis to Dirac basis followed by n blocking steps in the Dirac basis. This small expansion parameter implies the convergence of the rooted staggered theory to a local action of a single taste in the n → ∞ limit, .

(3.5)
It is important to note that in estimating a lower bound on the norm of (a c D inv,n ) −1 the finite renormalized fermion mass m R (a c ) provides the infrared cutoff of the Dirac spectrum when the volume is infinite. Since renormalization is multiplicative in the staggered formulation, it is implemented on the physical scale a c requiring the adjustment of the bare mass in the n → ∞ limit which is equivalent to the a f → 0 continuum limit. The choice m R (a c ) is arbitrary once a physical scale a c is set in the theory. This is expected because the RG invariant fermion mass related to m R (a c ) is arbitrary but as long as it is kept finite in the estimate of the expansion in Eq. (3.5) the convergence of the rooted theory to the local single taste action is assured. This line of reasoning, based on [35,37], explains why the sextet model of the rooted fermion doublet with finite renormalized mass is expected to be in the correct universality class when the continuum limit is taken. In our simulations of the volume dependent renormalized coupling g 2 R the renormalized fermion mass m R (a c ) is set to zero on any physical scale a c since the bare mass m itself is set to zero and the mass renormalization is multiplicative from the chiral symmetry of the staggered formulation. Since the estimate of the expansion for the convergence of the rooted determinant to the local single taste determinant as given in Eq. (3.5) is not applicable in the work presented here, the rooting procedure at m R (a c ) = 0 requires separate discussion.

Rooting in finite physical volume at zero bare mass
It is important to precisely define our rooting procedure where the required limit suggested by Eq. (3.5) is not followed. In all of the simulations reported in this paper the bare fermion mass is set to zero at finite lattice spacing a. This also sets the renormalized mass to zero on any choice of physical scale a c on the coarse lattice. Although the estimate on the bound and its scaling with the RG steps in Eq. (3.4) is lost, there is no problem with the rooting procedure since anti-periodic boundary conditions are imposed on the fermions in all four directions.
As we will now show, the choice of anti-periodic boundary conditions for the fermions restores the validity of the rooting procedure. The simulations always target some chosen values of the scale-dependent renormalized coupling. Each choice selects the corresponding linear size L of the physical volume in the continuum. The renormalized coupling g 2 R (a c ) in the finite volume L also depends on the ratio a c /L from a 1-parameter family of schemes. As the number of RG steps keeps increasing toward the continuum limit, the coupling on the scale a f (bare coupling on the cutoff scale) has to be adjusted while g 2 R (a c ) is held fixed, and similarly the lattice size measured in a f units is adjusted to keep the physical size L fixed together with a c /L. The scheme we introduced in section 2 defines a different but related finite volume scheme without affecting the reasoning. The former is built on the RG procedure and the other is defined on the gradient flow.
In the finite volume scheme a finite gap λ gap (a c ) is created in the Dirac spectrum which depends on g 2 R (a c ). Weak couplings correspond to small physical scales and the gap is approximately determined by the minimum momentum π/L in each direction with O(g 2 R (a c )) corrections. As the renormalized coupling become stronger with increasing volume and the interacting energy levels increasingly repel, they settle into a gradually decreasing but finite gap λ gap (a c ) set by the physical scale of the volume. The estimate for the bound on the taste breaking operator remains unchanged with a c ∆ n < ∼ a f /a c . The new estimate for the upper bound on the inverse of the taste invariant operator is given by (a c D inv,n ) −1 ≤ 1/(a c λ gap (a c )) with the gap of the spectrum set at the physical scale a c . The important combined estimate follows with D −1 inv,n ∆ n ≤ a f /(a 2 c λ gap (a c )) and the small expansion parameter is changed to .
The finite gap in the Dirac spectrum implies the convergence of the rooted staggered theory to a local action of a single taste in the n → ∞ limit, . (3.7) The conclusion from this simple analysis is that in the calculation of the volume dependent running coupling with a rooted and massless fermion doublet the role of the renormalized mass m R (a c ) on the physical scale is replaced by the λ gap (a c ) gap of the Dirac operator in the finite physical volume set by the targeted renormalized coupling g 2 R (a c ) for fixed a c /L. It is easy to see how this works at weak coupling and is sustained with growing volume if the gap does not collapse to zero at some critical volume size. At any targeted value of g 2 R (L) while holding the physical size L fixed, the eigenvalues of the infrared Dirac spectrum will collapse into degenerate quartets in the a f → 0 limit, consistent with the locality of the rooted action in the continuum limit.
We know, however, that although the simulations become increasingly difficult with increasing volume, the ensemble-averaged gap cannot disappear in finite physical volumes not even after some rapid crossover into the phase which either has chiral symmetry breaking or is instead conformal. With chiral symmetry breaking, the low end of the spectrum is expected to scale as λ ∼ 1/V which protects the gap. In the conformal theory the spectral density is expected to scale as ρ(λ) ∼ λ α with some critical exponent α and λ ∼ (1/L) 4/(1+α) for the low infrared part of the spectrum which also protects the gap from complete collapse. The finite gap cannot disappear in finite physical volumes even if the rooted model is conformal. In the conformal case the beta function is expected to turn zero at some critical coupling g 2 crit which can only be reached asymptotically at infinite volume. Our method with rooted staggered fermions can clearly distinguish a conformal model from one with chiral symmetry breaking.
The simulations, as reported in this paper, reach a limited range in the renormalized coupling without any sign of the β function turning zero. Beyond our reach, the results do not rule out conformality in large volumes, although they remain consistent with chiral symmetry breaking. Although it is tempting to pursue further simulations in large volumes at vanishing fermion mass, it is not practical at very small values of the gap in the Dirac spectrum.

The bridge to large volume physics and simulations at finite cutoff a f
Studying small fermion mass deformations in large volumes at finite cutoff a f can clearly differentiate between phases with chiral symmetry breaking, or conformality. Taste breaking at finite a f is described by operators in the Symanzik effective theory (SET) as calculated in [38,39]. As pointed out in [37], when the goal is to match the rooted theory to the Symanzik effective theory the a c m R (a c ) term can be dropped from the denominator in Eq. (3.3) since matching to the taste breaking operators is done at some finite momentum p ≫ Λ IR which serves in the matching loop diagrams as an IR cutoff.
The bound in Eq. (3.3) is much weaker than needed in the derivation of the SET, and it implies for infinite volume that the chiral m → 0 limit can only be taken after the continuum (a f → 0) limit. In Eq. (3.3) the role of m R (a c ) was to establish the existence of the correct continuum limit of the full rooted theory on any scale including the far infrared when the volume is infinite [35]. It follows from [37] that the Symanzik effective theory is well-defined in the chiral limit, together with the chiral effective theory that can be derived from the SET. The requirement that the zero mass limit for staggered fermions should be taken only after the continuum limit is then reproduced by calculations within staggered ChPT [31] for certain operators.
To bridge the current work with inherently non-perturbative large volume analysis we follow the procedure just outlined with mass deformed analysis at finite cutoff. What we observe is consistent with chiral symmetry breaking of the non-perturbative phase in large volumes. To build the bridge to the results in this paper we are interested in a scale-dependent and volume independent renormalized coupling in the symmetry breaking phase matching the scale dependent coupling g 2 R (a c ) presented here. This would leave no room for the β function turning zero on any scale. This strategy is outlined in more detail in [8] with results of a preliminary implementation.

Numerical simulation
The details of the simulations are similar to [29,58]. In particular we use the staggered fermion action with 4 steps of stout improvement [22] and stout parameter ̺ = 0.12. The bare fermion mass is set to zero, 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 [62,63]. For integration along the gradient flow we use both the Wilson plaquette and the tree-level improved Symanzik discretizations. The observable E(t) is discretized as in [25]. Hence, in the terminology of [64], we consider the discretizations W SC and SSC for Wilson-flow and tree-level improved Symanzik-flow, respectively.
As detailed in section 3 a gap in the Dirac spectrum is needed for the validity of rooting hence the available physical volume is limited. This translates into the limitation that the renormalized coupling cannot be explored above a certain value with a given set of lattice volumes. This limitation is however not unique to our running coupling scheme and not even unique to staggered fermions. All running coupling studies that are directly at the massless limit (by either setting the mass to zero using staggered or chiral fermions,  (1)  or tuning κ to the massless point κ c using Wilson fermions) will be limited to a certain renormalized coupling range with a given set of lattice volumes. This is because on a given set of lattice volumes a quite large renormalized coupling can only be achieved by increasing the bare gauge coupling which in turn will produce small Dirac eigenvalues which in turn will cause the (R)HMC algorithm to break down because the condition number of the Dirac operator might be very large on some configurations. We will see that in our scheme we are able to explore the range 0 < g 2 R < 6.5 which is however quite large and includes the location of the 3-loop and 4-loop fixed point in the MS scheme [65,66].
There is also a practical issue related to the rooting procedure. Rooting is implemented by the RHMC algorithm which relies on the Remez algorithm. The latter is used for the computation of the coefficients in the partial fraction expansion of the fourth root. A  necessary input for the Remez algorithm is an upper and lower bound on the spectrum of the Dirac operator squared D † D. For m > 0 a strict lower bound with staggered fermions is m 2 . However we set m = 0 and use the anti-periodic boundary conditions to produce a gap in the spectrum and no strict lower bound is available in this case. Hence we first need to measure the lowest and highest Dirac eigenvalues in all runs and then set the lower and upper bounds accordingly for the subsequent production runs. We found that this procedure is robust and a carefully chosen lower and upper bound on the spectrum is not violated in the production runs. Histories of the lowest eigenvalue for various parameters are shown for illustration in figure 1. As expected, increasing β leads to a larger lowest eigenvalue and similarly decreasing the lattice volume also leads to larger lowest eigenvalues. In a lattice setting a convenient and practical method of calculating the running coupling or its β-function is via step scaling [67,68]. In this context the finite volume L is increased by a factor s and the change of the coupling, (g 2 (sL) − g 2 (L))/ log(s 2 ), is defined as the discrete β-function. Note that in this convention asymptotic freedom corresponds to a positive discrete β-function for small values of the renormalized coupling. If the ordinary infinitesimal β-function of the theory possesses a fixed point, the discrete β-function will have a zero as well. Note that as s → 1 the discrete β-function turns into the infinitesimal variant. 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 volume pairs 8 4 → 12 4 , 12 4 → 18 4 , 16 4 → 24 4 , 20 4 → 30 4 and 24 4 → 36 4 . The continuum limit corresponds to L/a → ∞. Hence our data set has 5 pairs of lattice volumes over a range of lattice spacings to cover a desired range of renormalized couplings.
The collected number of thermalized unit length trajectories at each bare coupling and volume was between 2000 and 20000 depending on the parameters and every 10 th was used for measurements. The acceptance rates were between 65% and 95%. The measured renormalized coupling values are listed in tables 1 and 2 and the resulting discrete βfunctions are shown in figure 2 for the two discretizations we considered, SSC and W SC.
Clearly, at finite lattice spacing, or equivalently at finite lattice volume, the qualitative features of the two discretizations are quite different. While the discrete β-function is  It is important to point out that the observed zeros of the discrete β-functions of the W SC setup for the roughest four lattice spacings are however such that as the lattice spacing decreases, the location of the zero increases.
Let us emphasize that the behavior of the discrete β-function at finite lattice volume, whether it crosses zero or not, is entirely irrelevant as far as the continuum model is concerned. The measured data at finite lattice volume need to be continuum extrapolated and zeros of the discrete β-function may or may not survive the continuum limit. It will turn out in the next section that in fact the zeros of the W SC setup do disappear in the continuum limit while there aren't any zeros to begin with in the SSC setup, and the continuum results for the W SC and SSC setups agree, as they should, and show no sign of a fixed point in the explored coupling range.

Continuum extrapolation
The simplest way to perform the continuum extrapolation of our data is to interpolate the renormalized coupling, g 2 (β), as a function of the bare coupling β at each lattice volume. We choose the interpolating functions as similarly to [69]. The order of the above polynomial is allowed to be n = 3, 4 or 5 for the volumes L/a = 8, 12, 16, 18, 24 and n = 3 or 4 for the volumes L/a = 20, 30, 36. The corresponding degrees of freedom of the fits are 1, 2 or 3 for the first set and 1 or 2 for the second set. Once the parametrized curves g 2 (β) are obtained for all volumes the discrete β-function (g 2 (sL) − g 2 (L))/ log(s 2 ) can be computed for arbitrary g 2 (L) for fixed L/a and s = 3/2. Then assuming that corrections are linear in a 2 /L 2 the continuum extrapolation can be performed for each g 2 (L).
In [64,70] we calculated the tree-level improvement of our observables in order to have smaller slopes in the continuum extrapolations. For the SU (3) fundamental model with N f = 4 tree-level improvement did indeed decrease the slopes over the full considered coupling range, however with N f = 8 we observed in [58] that tree-level improvement only decreased the slopes for small couplings but in fact increased it for larger couplings. In the current work we observe the same. More precisely for approximately g 2 (L) 3.0 tree-level improvement decreased the absolute value of the slope of the continuum extrapolation but for g 2 (L) 3.0 it increased it. The reason most probably is the same as for N f = 8, namely that the large fermion content enhances the fermion loops which are completely absent from the tree-level calculation and these fermion loops are bound to increase with increasing coupling. For this reason we do not include tree-level improvement in the current work because the phenomenologically interesting region is in the larger coupling range, g 2 (L) ∼ 6.
At small values of the renormalized coupling the continuum discrete β-function can be reliably calculated in continuum perturbation theory. For the SU (3) sextet model with N f = 2 we have, where b 1 = 13/3 and b 2 = −194/3. Due to the small volume gauge dynamics mentioned in section 2 only the first coefficient is the same in our finite volume scheme as above, but nevertheless for comparison we show both the 1-loop and 2-loop expressions. The numerical results, after continuum extrapolation, should agree with the perturbative result for small renormalized coupling and this test is an important cross-check of our procedures. Following the above procedure with a fixed polynomial order for each volume for the interpolations one obtains a continuum result for both the SSC and W SC setups. The interpolations (5.1) are linear in the free parameters hence the statistical errors are easy to propagate to the final result. Of course one needs to make sure that the continuum extrapolations are acceptable from a statistical point of view, for example the χ 2 /dof values are not very large, and one needs to test whether all 5, or perhaps only 4, or perhaps only 3 lattice spacings are in fact in the scaling region. The more lattice spacings that are useable in the continuum extrapolation, the more reliable the result is.

Systematic error estimate
Apart from the statistical errors we would like to estimate the systematic errors too as precisely as possible. The only source of systematic error is the continuum extrapolation. However two distinct types of systematic errors are present in our procedures. One, various polynomial orders can be used for the interpolation (5.1) for each lattice volume and two, one may perform the continuum extrapolation using 5 or 4 lattice spacings (assuming of course that all 5 lattice spacings are actually in the scaling region), i.e. dropping the roughest lattice spacing. As we discussed in section 4 the rooting trick of the staggered formulation itself does not introduce unwanted systematic effects.
We will apply the histogram method [71] in order to estimate the systematic uncertainties. The polynomial order n for the interpolation (5.1) is allowed to be n = 3, 4, 5 for L/a = 8, 12, 16, 18, 24 and n = 3, 4 for L/a = 20, 30, 36. All together this leads to 3 5 · 2 3 = 1944 interpolations and correspondingly to 1944 continuum results for a given discretization. Following [72] a Kolmogorov-Smirnov test is applied to the 1944 interpolations and only those are deemed acceptable to which the Kolmogorov-Smirnov test assigns at least a 30% probability. This requirement results in 240 and 306 acceptable interpolations for the SSC and W SC cases, respectively. These all correspond to continuum extrapolations using 5 lattice spacings.
In order to include the systematic effect coming from performing continuum extrapolations using 4 lattice spacings only, i.e. dropping the roughest, 8 4 → 12 4 , we include such extrapolations too. Using the volumes L/a = 12, 16,18,20,24,30,36 only with the polynomial orders as above, we have a total number of 3 4 · 2 3 = 648 interpolations. Out of these the Kolmogorov-Smirnov test allows 240 and 249 for the SSC and W SC cases, respectively 1 .
Summarizing the above, we have 240 + 240 = 480 continuum results for the SSC case and 306 + 249 = 555 continuum results for the W SC case. In both cases these are binned into a weighted histogram where the weights are given by the Akaike Information Criterion (AIC) [73][74][75]. We take 68% of the full distribution around the average to estimate the systematic error. Further details are given in [58] where it is explained how to perform the  Kolmogorov-Smirnov test in a running coupling setup and also for the precise definition of the AIC weights.

Final results
In show that the final continuum result using the W SC data actually agrees within errors with the one obtained using the SSC data. This agreement between two different discretizations is a reassuring consistency check of our procedures, especially because we have seen in the W SC case in figure 2 that at the smaller lattice volumes the β-functions did cross zero. A remnant of the small lattice volume β-functions crossing zero is that for approximately g 2 (L) 5.0 some of the β-function values that are used in the extrapolation are negative. But it is clear from figure 4 that the continuum value is positive for both g 2 (L) = 5.0 and 6.0 and in fact over the entire range. The zeros of the small volume β-functions hence did not survive the continuum limit and the W SC and SSC final results agree within errors.
It is worth emphasizing again: in a given discretization the finite (perhaps small) volume discrete β-functions can perfectly well cross zero while in another discretization the same thing may not happen. This in itself however is in no way indicative of the behavior in the continuum as these small volume zeros may disappear in the continuum. The present model, SU (3) gauge theory with N f = 2 flavors of sextet massless fermions using the W SC and SSC discretizations serves as an example.
Another cautionary note is in order regarding small volumes. It is clear from figures 3 and 4 that for approximately g 2 (L) 5.5 using only the 3 roughest lattice spacings, 8 4 → 12 4 , 12 4 → 18 4 and 16 4 → 24 4 would in fact give a continuum result compatible with the one including all 5 lattice spacings in the SSC setup. Only at around g 2 (L) ∼ 6.0 the 3 roughest lattice spacings alone are not usable in a continuum extrapolation. The same, however, is not the case for the W SC discretization. Already for approximately g 2 (L) 2. essential, without these one may obtain a much smaller β-function which actually would be totally unreliable. This is all the more important in the phenomenologically important larger coupling region g 2 (L) ∼ 6.
In [59][60][61] preliminary results were reported on the sextet model favoring an infrared fixed point in the continuum. The scheme used is the same as ours except that the fermions were anti-periodic in one direction only. The lattice discretization was different, Wilson fermions were used, however the largest lattice volumes used in the step function were 16 4 → 24 4 . We speculate the reported fixed point was a lattice artefact due to the large cut-off effects inherent in using small lattice volumes. This scenario would be analogous to some extent with our W SC setup and using only our roughest 3 lattice spacings. Similarly, the inconclusive findings in [20,21] we speculate are also the result of using small lattice volumes without reaching the scaling regime.
If one were to work with a fixed discretization one of course would not know a priori how large volumes are needed for a reliable continuum extrapolation. That is why it is extremely important to consider several discretizations, check that the lattice spacings are in fact in the scaling region, estimate the systematic uncertainty coming from the continuum extrapolation reliably, and only trust results in the continuum if they agree for the considered discretizations. In our work we have performed this analysis in a fully controlled fashion.
We show the final continuum result in figure 5. Clearly the β-function stays positive over the entire range and is monotonically increasing, and agreement is found for g 2 (L) < 2.5 between our result and the 2-loop perturbative result within 1.3σ.

Conclusion and outlook
We have studied SU (3) gauge theory coupled to N f = 2 flavors of massless Dirac fermions in the sextet representation. Our primary motivation was the fact that this model may be a possible realization of a composite Higgs scenario. The current work is an integral part of our program to understand the infrared dynamics of this model and its ultimate viability as a building block of Beyond Standard Model physics. Using lattice simulations we have studied the hadron spectrum and low energy behavior of the model in large volumes in previous work and concluded that the model is consistent with spontaneous chiral symmetry breaking. The current work supports this picture in that the running coupling, when carefully continuum extrapolated, shows no sign of a fixed point in the range 0 < g 2 < 6.5. This range in fact includes the zero of the 3-loop and 4-loop β-function in the MS scheme which is g 2 ≃ 6.28 and g 2 ≃ 5.73, respectively [65,66]. Even though our scheme is of course different from MS and we were only able to explore a limited renormalized coupling range we are tempted to speculate that the zero of the 3-loop and 4-loop β-function is a perturbative artefact just as the zero of the 2-loop β-function [76,77] is a perturbative artefact at g 2 ≃ 10.58. The absence of a fixed point and hence chiral symmetry breaking in the infrared would be consistent with the ladder resummation in the Schwinger-Dyson approach predicting non-conformal behavior [1,78]. What we can firmly conclude from our work is that in the scheme we use, a fixed point is ruled out in the range 0 < g 2 < 6.5 with high confidence. In future work we would like to connect the running coupling at some large renormalized coupling value to another scheme we define in nearly infinite volumes using the flow time as running scale µ = 1/ √ 8t. Since our large volume simulations are in the chirally broken phase this connection would follow the coupling from the perturbative regime all the way to the chirally broken regime, ruling out conformal behavior completely.
It is important to emphasize that only non-perturbative lattice calculations are able to reliably answer questions about the infrared dynamics of non-abelian gauge theories such as our sextet model. These lattice calculations are nevertheless plagued by systematic uncertainties and reliable results are only possible to obtain if these are fully controlled. In our work we were able to fully control all systematics leading to reliable continuum results. As far as the sextet model is concerned, our present work is the first to do so.
The lack of a reliable and fully controlled continuum extrapolation prior to our work made it possible that seemingly contradictory claims appeared in the literature. These were not contradictory in a sense that at finite lattice spacing different discretizations may indeed lead to different conclusions. A contradiction only arises if results at finite lattice spacing are interpreted to reflect the properties of the continuum model. In the continuum of course all correct discretizations should agree and all conclusions should converge.