Infrared fixed point of SU(2) gauge theory with six flavors

We compute the running of the coupling in SU(2) gauge theory with six fermions in fundamental representation of the gauge group. We establish an infrared stable fixed point at strong coupling and measure also the anomalous dimension of the fermion mass operator at the fixed point. This theory therefore likely lies close to the boundary of the conformal window and will display novel infrared dynamics if coupled with the electroweak sector of the Standard Model.


I. INTRODUCTION
Determination of the vacuum phase of an SU(N ) gauge theory as a function of the number of massless flavors of Dirac fermions, N f , and their representations presents a challenge for our basic understanding of gauge theory dynamics at strong coupling. A lot of effort in the field of lattice gauge theory has been devoted to address the existence and properties of infrared fixed point (IRFP), which appears when N f is between a critical lower limit N crit f and the loss of asymptotic freedom. The bounds depend on N and the fermion representation. For recent reviews see [1][2][3]. A much studied benchmark case is SU(2) gauge theory with two Dirac fermions in the adjoint representation [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18], where the results indicate the existence of an IRFP.
In SU(2) gauge theory with fermions in the fundamental representation the precise dependence on N f remains uncertain, despite a large number of recent studies on the lattice [19][20][21][22][23][24]. The upper edge of the conformal window is robust: the asymptotic freedom is lost at N f = 11, where the 1-loop β-function coefficient changes sign. Just below the upper edge, at 10 flavors the theory has a perturbatively stable Banks-Zaks type infrared fixed point [25], which has also been observed on the lattice [21]. Recently, simulations of the 8 flavor theory have also shown the existence of a fixed point [24]. On the other hand, the theory with N f = 2 is well below the conformal window and breaks the chiral symmetry according to the expected pattern, and the theory with N f = 4 is expected to fall within this category as well [21]. However, for N f = 6 the previous results remain so far inconclusive [20][21][22][23].
In this article we give strong evidence that the six flavor theory indeed has an IRFP at strong coupling. * viljami.leino@helsinki.fi † kari.rummukainen@helsinki.fi ‡ joni.suorsa@helsinki.fi § kimmo.i.tuominen@helsinki.fi ¶ sara.tahtinen@helsinki.fi The result is based on a thorough state-of-the-art measurements of the running coupling and the anomalous dimension of the fermion mass operator. We use the HEX smeared Wilson-clover fermion lattice action and measure the coupling using the Yang-Mills gradient flow [26,27] in conjunction with the finite volume step scaling function with Dirichlet ("Schrödinger functional") boundary conditions [28]. The value of the coupling at IRFP, g 2 * , is scheme dependent and hence depends on the gradient flow time. With our benchmark scheme we find g 2 * = 14.5(4) +0.4 −1.2 with statistical and systematic errors.
We also measure two scheme independent quantities at the IRFP: the mass anomalous dimension γ * m and the leading irrelevant critical exponent γ * g , which gives the slope of the β-function at IRFP. The mass anomalous dimension is measured using two different methods: the mass step scaling method [29] and the Dirac operator spectral density method [30]. At the fixed point we observe γ * m = 0.283(2) +0.01 −0.01 . The slope of the β-function is directly measurable from the step scaling function of the coupling, obtaining γ * g = 0.648(97) +0. 16 −0.1 . In contrast to the value of the fixed point coupling, we observe that γ * g remains independent of the gradient flow time, in accord with the scheme independence of this quantity. This paper is structured as follows: In section II we define the model and outline the simulation methods. The numerical results are presented for running coupling, leading irrelevant exponent, and mass anomalous dimension, in sections III, IV, and V respectively. In section VI we present our conclusions.

II. LATTICE IMPLEMENTATION
The model we use and the algorithmic details we apply are described in detail in [17,24], and our discussion here will be brief so that we can then focus on the results we obtain in the case of N f = 6. The model is defined by the lattice action

arXiv:1707.04722v2 [hep-lat] 6 Apr 2018
The smeared gauge link V is defined by hypercubic truncated stout smearing (HEX smearing) [31], and we mix smeared, S G (V ), and unsmeared, S G (U ), Wilson gauge actions with mixing parameter c g = 0.5. This partial smearing allows us to reach significantly larger couplings by avoiding the unphysical bulk phase transition in the region of interest of the parameter space [32]. We use clover Wilson fermion action with the Sheikholeslami-Wohlert coefficient set to tree level value of unity, c SW = 1, which is the standard choice for smeared clover fermions. We have verified that this value is very close to the true non-perturbatively fixed c SW coefficient, canceling most of the O(a) errors.
On a lattice of size L 4 we use Dirichlet boundary conditions at the temporal boundaries x 0 = 0, L by setting the gauge link matrices U = V = 1 and the fermion fields to zero. The spatial boundaries are periodic. These boundary conditions enable simulations at vanishing fermion mass, and allow the mass anomalous dimension to be measured using the same configurations as for the running coupling.
We run our simulations using the hybrid Monte Carlo algorithm with 2nd order Omelyan integrator [33,34] and chronological initial values for the fermion matrix inversions [35]. We tune the step length to have an acceptance rate larger than 85%. We run the simulations with bare couplings varying within the range and tune the hopping parameter κ c (β L ) so that the absolute value of the PCAC fermion mass [36] is less than 10 −5 at lattices of size 24 4 . The same critical hopping parameter values are used for all the lattice sizes, and for each β L (and corresponding κ c (β L )). The critical hopping parameters are available in the Table III. We use lattices of size L =8, 10,12,16,18,20,24, and 30 chosen to allow step scaling with either s = 2 or s = 3/2. For our analysis we choose s = 3/2 as it includes more pairs within the larger lattices. We generate (5 − 100) · 10 3 trajectories for each combination of β L and L. For the exact number of trajectories used, see Table IV.
To define the running coupling, we apply the Yang-Mills gradient flow method [26][27][28]. This method defines a flow that smooths the gauge fields and removes UV divergences and automatically renormalizes gauge invariant objects [37]. The method is set up by introducing a fictitious flow time t and studying the evolution of the flow gauge field B µ (x, t) according to the flow equation where G µν (x; t) is the field strength of the flow field B µ and D µ = ∂ µ + [B µ , · ]. The initial condition is B µ (x; t = 0) = A µ (x), where A µ is the original continuum gauge field. In the lattice formulation the (unsmeared) lattice link variable U replaces the continuum flow field, which we then evolve using either the tree-level improved Lüscher-Weisz pure gauge action (LW) [38] or the Wilson plaquette gauge action (W). In the continuum limit these should yield identical results, providing a check of the reliability of the limit. The coupling at scale µ = 1/ √ 8t [39] is defined via energy measurement as where a is the lattice spacing. The shift parameter τ 0 is introduced to reduce the O(a 2 ) discretization effects caused by the flow [40] and can be numerically estimated during the analysis. The normalization factor N for the chosen boundary conditions has been calculated in [41] to match the MS coupling in the tree level. As the translation symmetry is broken by the chosen boundary conditions, the coupling g 2 GF is measured only on the central time slice x 0 = L/2. To quantify the effects of different discretizations, we measure the energy density E(t) using both symmetric clover and simple plaquette discretizations of the flow.
In order to limit the scale into a regime where both lattice artifacts and finite volume effects are minimized, we relate the lattice and the renormalization scales by defining a dimensionless parameter c t such that µ −1 = c t L = √ 8t as described in [41,42]. The chosen boundary conditions have a reasonably small cutoff effects and statistical variance within the range of c t = 0.3 − 0.5 [41]. Value of this parameter defines the renormalization scheme.

III. EVOLUTION OF THE COUPLING
Our "benchmark" set of results presented here are obtained with gradient flow Eq. (2) evolved with Lüscher-Weisz action (LW), clover definition of energy density Eq. (3) and c t = 0.3. In order to estimate systematic errors, we vary discretizations of the flow and the observable and the parameter values of the flow. The raw data is available in tables V to XII.
The measured couplings with the aforementioned parameters are shown in the top panel of Fig. 1. It is clear from the figure that the finite volume effects become substantial on smaller lattices as the coupling grows larger.
To quantify the running of the coupling we use the finite volume step scaling function [43]: which describes the change of the measured coupling when the linear size of the system is increased from L to sL. The step scaling obtained from the data in Fig. 1 is shown in Fig. 2, using a scaling factor s = 3/2, that is pairs of lattices with L/a = 8 − 12, 12 − 18, 16 − 24 and 20 − 30. At the smallest volume pair L/a = 8 − 12 the  step scaling deviates significantly from the others, and will be excluded from the continuum analysis.
For the continuum limit σ(g 2 GF ) of the step scaling function we use the extrapolating function At weak coupling the cutoff effects are regulated by the proximity of the ultraviolet fixed point, and the lowest order discretization effects of the Wilson-clover action are expected to be of order O(a 2 ), motivating the use of Eq. (5). At large coupling, so long as the coupling remains below the possible IRFP, the continuum limit is ultimately also reached at the UV fixed point. However, due to the smallness of the β-function this would require astronomically large scale hierarchy between the lattice size L and lattice spacing a and hence is impossible to observe in simulations. Nevertheless, if the anomalous exponents of the fields remain small near the infrared fixed point, one can assume that the power counting of operators is applicable and the cutoff effects (dominated by dimension 6 operators) decrease with a power of the lattice spacing a. The naive a 2 behaviour may be modified by anomalous exponents, though. In section V we observe that the mass anomalous dimension at the IRFP remains relatively small, γ m ≈ 0.28, suggesting that the a 2 behaviour in Eq. (5) may also receive only minor corrections. The available range in our data does not allow us to numerically determine differences from Eq. (5) and hence we use it at all couplings. Near the IRFP we indeed observe that the continuum limit becomes somewhat less robust, which is taken into account in our systematic error estimation.
In order to determine the continuum limit of the step scaling function we need the measurements of the step scaling at constant value of the coupling. However, in practice the simulations are carried out at a fixed set of bare lattice couplings which do not correspond to same value of g 2 GF when a/L is varied. Hence, we adopt the customary interpolation procedure of the measured couplings to intermediate couplings. We do this using a polynomial fit 1 where we use m = 10 for lattices smaller than L = 16 and m = 9 for the larger lattices. With this choice we obtain the χ 2 /d.o.f's reported in Table I for each used c t and in Table XVI for each used discretization. We study the robustness of the fits by repeating the analysis with m decreased by one. While this choice increases the χ 2 /d.o.f., the results stay compatible with those obtained with larger m. In Fig. 3 we show the continuum limit extrapolation of the step scaling function when g 2 GF is varied from weak to strong coupling, obtained using Luscher-Weisz or Wilson flow actions and clover or plaquette field strength observables. At small couplings the continuum limit is very well  under control: different discretizations extrapolate very close to the same value. At couplings g 2 GF > ∼ 10 the continuum limits start to show a few per cent scatter. This is taken into account in the systematic uncertainties of the final results.
The τ 0 -correction parameter in Eq. (3) can be tuned to reduce most of the O(a 2 ) errors from the continuum limit extrapolation of the step scaling function Eq. (5). The parameter τ 0 should have a small effect in the continuum extrapolation, so long as it is not too large [44].
In practice, we have observed that τ 0 which depends logarithmically on g 2 GF works well at small coupling [24]. With c t = 0.3, Luscher-Weisz flow action and clover field strength observable we use which makes the interpolation errors almost vanish at g 2 GF < ∼ 10, as can be observed in Fig. 3. At larger couplings τ 0 correction cannot remove O(a 2 ) significantly without ruining the continuum limit. We note that in order to have consistent O(a 2 ) shift in the step scaling analysis, the τ 0 correction should be a function of g 2 GF instead of the bare coupling g 2 0 [28]. Because adjusting τ 0 changes the value of the measured g 2 GF , the final value of τ 0 is found by iterating equations (3) and (7), starting from the initial value g 2 GF = g 2 0 . In Fig. 4 we show the continuum limit of the step scaling at c t = 0.35, 0.4 and 0.45, evaluated at the IRFP of each c t . The couplings g 2 GF = g 2 * at the fixed point are shown in Table II. Because different c t -values correspond to different coupling constant scheme, the values of g 2 * vary significantly. It is evident that as c t is increased the difference between the Lüscher-Weisz and Wilson flow actions grows, contributing to increasing systematic errors.  6). It is evident that as c t increases the reliability of the continuum limit extrapolation decreases, and already at c t = 0.4 the final result has clearly unphysical strong "wavy" structure. This is caused by the use of the polynomial interpolation functions in Eq (6). We note that if we would use polynomials of smaller degree (for example, m = 8 in (6)) the wavy structure would be strongly reduced and error bands would be much narrower; however, the χ 2 /d.o.fvalues would not be acceptable.
Each value of c t corresponds to different coupling constant scheme, and the value of the fixed point coupling strongly depends on the value of c t . The location of fixed point for each c t is reported in Table II and individually for all discretizations of the flow in Table XIV in the appendix. The first error is the statistical uncertainty and the second error estimates systematic effects by including the full range of different discretization choices that were present in the Fig. 3. For our benchmark value c t = 0.3 we find that our continuum extrapolated results are compatible within 1σ level with respect to all these effects except in the interval g 2 GF ∈ [8,12] where LW and

IV. LEADING IRRELEVANT CRITICAL EXPONENT
We can also obtain the leading irrelevant exponent γ * g at the fixed point, defined by the slope of the β-function at the IRFP. This quantity is scheme independent, and thus should not depend on c t .
In the proximity of the fixed point we can approximate the β-function as Measuring the slope of the step scaling function σ(g 2 ) around the fixed point gives the exponent γ * g = 0.648(97) +0. 16 −0.1 at c t = 0.3; the results with other c t are shown in Table II. While there is noticeable variance between different discretizations, as indicated by the second set of errors, the result is compatible with the recent scheme independent estimate of γ * g = 0.6515 in Refs. [45,46]. The results obtained with different discretizations are shown individually in Table XV. When c t is varied, the value of γ * g remains constant within errors, in accord with the scheme independence of this quantity.
The results obtained above rely on the accurate continuum limit of the step scaling function σ(g 2 ). However, as discussed in section III, the continuum limit may be in effect somewhat different at the close proximity of the IRFP. We can verify the consistency of the results by using a finite size scaling method developed in [47][48][49][50] to get an alternative measurement of γ * g . In the close proximity of the IRFP, by integrating (8) we obtain a finite size scaling relation between lattices of size L ref and L [49]: This equation relies on the evolution of the coupling towards the fixed point as the lattice size is increased from L ref to L. Hence, it cannot be used exactly at the fixed point where there is no evolution, but only in some environment around it. We note that this also assumes vanishing discretization artifacts, and thus it can be used only if the lattices are already close enough to the continuum (L large). In Fig. 6 we show the fit to Eq. (9) to individual measurements of g 2 GF at β L ≤ 0.8, corresponding to measurements which are close to the fixed point. A good fit to Eq. (9) is obtained if we choose L ref /a ≥ 18, allowing us to extract an estimate for γ * g . Instead of using individual measurements we use the interpolated values of g 2 GF (β L , L), because this allows us to freely tune the value of β L . In Fig 7 we show the resulting γ * g from fits to Eq. (9), plotted as functions of g 2 ref ≡ g 2 GF (β L , L ref ). The red lines correspond to the values given in Table XV, measured from the slope of the step scaling function. The shaded error bands correspond to statistical errors while keeping the values of g 2 * fixed to central values in Table XIV, and the dashed lines show the variation of the result if we allow g 2 * to vary within the statistical error range.
The resulting γ * g is expected to be close to the true γ * g only in close proximity of the IRFP. However, too close to  the IRFP Eq. (9) becomes unstable, which is indicated by a sudden drop in the γ * g measurements. Indeed, at g 2 ref ≈ 12 we observe the c t = 0.3 case to give γ * g which is in agreement with the one obtained from the slope of the β-function. At small g 2 ref the measurement of γ * g using Eq. 9 approaches zero.

V. ANOMALOUS DIMENSION OF THE MASS OPERATOR
In order to measure the anomalous dimension of the fermion mass operator γ * m we use two different methods, the mass step scaling method and the spectral density method. In the step scaling method we measure γ m from the running of the pseudoscalar density renormalization constant [29,51] where f P and f 1 are pseudoscalar current densities defined explicitly in e.g. [24,36]. The mass step scaling function is defined as [29]: Σ P (u, s, L/a) = Z P (g 0 , sL/a) Z P (g 0 , L/a) g 2 GF (g0,L/a)=u (11) As in the case of the coupling, we choose s = 3/2. The continuum limit σ P (u, s) is obtained by interpolating the measured Z p by 8th order polynomials and assuming O(a 2 ) errors. The mass anomalous dimension is then obtained as [51] γ The results are shown in Fig. 8 and the raw data is given in Table XIII. The method gives results comparable to one loop perturbation theory predictions at small gauge coupling g 2 GF . However, the method becomes unstable at large coupling, which implies that at the fixed point g 2 * ≈ 14.5 the continuum limit cannot be trusted. The second way to measure γ m is based on the fact that the it also determines the scaling of the spectral density of the massless Dirac operator. The explicit calculation of the eigenvalue distribution is prohibitively costly, but the stochastic methods [52] have made it possible to determine the mass anomalous dimension from the scaling of the mode number of the Dirac operator [30]. The mode number is known to follow a scaling behavior in some energy range between the infrared and the ultraviolet in the vicinity of a fixed point. Here γ * m is the mass anomalous dimension γ m at the fixed point.
We calculate the mode number per unit volume of Eq. (13) by using where the operator P(Λ) projects from the full eigenspace of M = m 2 − / D 2 to the eigenspace of eigenvalues smaller than Λ 2 . The trace is evaluated stochastically [52], and fitted to the power law behavior of Eq. (13). However, the energy range where this power law behavior holds is not known beforehand, and needs to be determined by observing the quality of the fit in a given range. We use L/a = 24 lattices from the step scaling analysis, and take 12 to 20 well separated configurations for each value of the gauge coupling. We calculate the mode number for 90 values of Λ 2 ranging from 10 −3 to 0.3. The results are then fitted to Eq. (13). The fit range is determined by varying its lower and the upper limits and observing the stability and the quality of the fit. As a cross reference at weak coupling, the fitted value of γ * m and the value obtained with the step scaling method are compared.
In Fig. 9 we plot the mode number divided by the fourth power of the eigenvalue scale, where the the fit range and the fit function of Eq. (13) are shown overlaid in red. According to Eq. (13) in the proximity of the fixed point the infrared behavior should be a power law in the absence of lattice artifacts. We observe this at strongest couplings, however, at small couplings the low eigenvalues appear in discrete energies, which manifests as the bumps in the mode number curve, making the power law less evident. To illustrate the evolution of the mass anomalous dimension we use the same fit range for both weak and strong couplings.
The final result of the spectral density method is shown in Fig 10, where the mass anomalous dimension γ * m , obtained by fitting the data with Eq. (13), is shown as a function of the gauge coupling g 2 GF . The shaded band illustrates the uncertainty resulting from varying the upper and lower limits of the fit range by ∼ 50%. The  Fig. 9 is shown with black points and the one loop perturbative result with a red line. The shaded regions are estimates for reasonable ranges of values obtainable using the method, and were obtained by varying the fit range shown in Fig. 9 slightly.
largest uncertainty arises at small gauge couplings, where the bumps in the data cause the changes in the fit range to change the fit dramatically. The error band of Fig. 10 becomes narrower towards the larger couplings as the ensembles near the IRFP are less sensitive to variations of the fit range.
At the fixed point g 2 * = 14.5 we obtain γ * m = 0.283(2) +0.01 −0.01 . However, this result is obtained at fixed lattice size L/a = 24. A proper continuum limit requires extrapolation to infinite L, but at smaller L/a the finite size effects make the usable range for the power law fit too narrow.
Interestingly, the mass step scaling method and the spectral density method complement each other: while the mass step scaling is stable and accurate at weak couplings, where the spectral density method fails, at strong coupling the roles are reversed.

VI. CONCLUSIONS
We have studied the running coupling in the SU(2) lattice gauge theory with 6 fermions in the fundamental representation. Gradient flow algorithm with Dirichlet boundaries was shown to provide robust results on the large coupling behavior of this theory allowing to conclude the existence of IRFP at g 2 * = 14.5(4) +0.4 −1.2 in our benchmark scheme. The scheme-independent slope of β-function at IRFP was measured to be γ * g = 0.648(97) +0. 16 −0.1 . We also determined the mass anomalous dimension γ m in this theory using the spectral density method and the mass step scaling method. The step scaling method gives results compatible with perturbation theory at weak coupling. At intermediate couplings the step scaling method can be matched on the results from the spectral density method which remain stable in the vicinity of the fixed point where the step scaling computation breaks down. With the spectral density method we estimated the mass anomalous dimension at the fixed point as γ * m = 0.283(2) +0.01 −0.01 , albeit a proper continuum limit is still lacking.
Our results establish that the SU(2) gauge theory with six fermion flavors in the fundamental representation is within the conformal window. The theory likely is near the lower boundary of the conformal window, which makes it an interesting candidate for beyond the standard model theories: when coupled with the electroweak gauge currents the chiral symmetries are explicitly broken, the theory is pulled outside the conformal window and may constitute a concrete example of a walking technicolor theory.    (17) (4) 11.86 (7) 11.39 (2) 10.94 (2) 10.98 (2) 11.09 (6) 11.37 (6) 11.77 (10)