Safety versus triviality on the lattice

We present the first numerical study of the ultraviolet dynamics of nonasymptotically free gauge-fermion theories at large number of matter fields. As test bed theories, we consider non-Abelian SU(2) gauge theories with 24 and 48 Dirac fermions on the lattice. For these numbers of flavors, asymptotic freedom is lost, and the theories are governed by a Gaussian fixed point at low energies. In the ultraviolet, they can develop a physical cutoff and therefore be trivial, or achieve an interacting safe fixed point and therefore be fundamental at all energy scales. We demonstrate that the gradient flow method can be successfully implemented and applied to determine the renormalized running coupling when asymptotic freedom is lost. Additionally, we prove that our analysis is connected to the Gaussian fixed point as our results nicely match with the perturbative beta function. Intriguingly, we observe that it is hard to achieve large values of the renormalized coupling on the lattice. This might be an early sign of the existence of a physical cutoff and imply that a larger number of flavors is needed to achieve the safe fixed point. A more conservative interpretation of the results is that the current lattice action is unable to explore the deep ultraviolet region where safety might emerge. Our work constitutes an essential step toward determining the ultraviolet fate of nonasymptotically free gauge theories.


I. INTRODUCTION
Asymptotically free theories [1,2] are fundamental according to Wilson [3,4] since they are well defined from low to arbitrary high energies. This remarkable property stems simply from the fact that asymptotically free theories are obviously conformal (and therefore scale invariant) at short distances, given that all interaction strengths vanish in that limit. This is one of the reasons why asymptotic freedom has played such an important role when building extensions of the Standard Model (SM). Non-Abelian gauge-fermion theories, like QCD, with a sufficiently low number of matter fields are time-honored examples of asymptotically free gauge theories. These theories feature a single four-dimensional marginal coupling induced by the gauge dynamics, and no further interactions are needed to render the theories asymptotically free.
On the other hand, purely scalar and scalar-fermion theories are not asymptotically free. Adding elementary scalars, and upgrading gauge-fermion theories to gauge-Yukawa systems, one discovers that scalars render the existence of asymptotic freedom less guaranteed. In particular, complete asymptotic freedom in all marginal couplings is no longer automatically ensured by requiring a sufficiently low number of scalar and fermion matter fields.
To determine the asymptotically free conditions on the low energy values of the accidental couplings that may lead to complete asymptotic freedom, a one-loop analysis in all marginal couplings is sufficient. One discovers that the gauge interactions are essential to tame the unruly behavior of the accidental couplings provided the latter start running within a specific region in coupling space at low energies.
Nonasymptotically free theories can belong either to the trivial or the safe category. Triviality occurs when the theories develop a physical cutoff and can therefore be viewed as low energy effective descriptions of a more fundamental but typically unknown quantum field theory. Triviality literally means that the only way to make sense of these theories as fundamental theories (when trying to remove the cutoff) is by turning off the interactions. In truth, it is rather difficult to demonstrate that a theory is trivial beyond perturbation theory since the couplings become large in the UV and therefore one can imagine a nonperturbative UV fixed point to emerge leading to nonperturbative asymptotic safety. Nevertheless, we can be confident that certain theories are indeed trivial. Perhaps the best known example is the four-dimensional λϕ 4 theory, which was studied on the lattice in a series of papers by Lüscher and Weisz [5][6][7]. The triviality here was established using large order hopping parameter expansion with perturbative renormalization group evolution, and corroborated with several lattice Monte Carlo simulations (see, for example, Refs. [8,9]). Related to our study here, the triviality versus conformality has also been investigated in large-N f QCD with staggered fermions [10].
Analytically, using a-maximization and violation of the a-theorem, one can also demonstrate that certain nonasymptotically free supersymmetric gauge-Yukawa theories such as super QCD with(out) a meson and their generalizations are trivial once asymptotic freedom is lost by adding enough super matter fields [11].
We move now to safe theories. These achieve UV conformality while remaining interacting, meaning that the interaction strengths freeze in the UV without vanishing. The first four-dimensional safe field theories were discovered in Ref. [12] within a perturbative study of gauge-Yukawa theories in the Veneziano limit of a large number of flavors and colors. This discovery opened the door to new ways to generalize the Standard Model as envisioned first in Ref. [13] and then investigated in Refs. [14][15][16][17][18][19][20] with impact in dark matter physics [21,22] and cosmology [18]. Additionally, it allowed us to use the large charge [23,24] method to unveil new controllable conformal field theory (CFT) properties for four-dimensional nonsupersymmetric quantum field theories [25].
Although both safe and free theories share the common feature of having no cutoff, the respective mechanisms and dynamics for becoming fundamental field theories are dramatically different [12]. For example, within perturbation theory, it is impossible to achieve safety with gauge-fermion theories [26]. Yukawa interactions and the consequent need for elementary scalars is an essential ingredient to tame the UV behavior of these theories. Of course, this is a welcome discovery given that it provides a pleasing theoretical justification for the existence of elementary scalars, such as the Higgs, and Yukawa interactions without the need to introduce baroque symmetries such as supersymmetry. Beyond perturbation theory, however, little is known, and it is therefore worth asking whether scalars are needed to achieve asymptotic safety.
Interesting hints come from the knowledge of the beta function for Abelian and non-Abelian gauge-fermion theories at leading order in 1=N f [27][28][29][30]. This result suggests, as we shall review below, the potential existence of short distance conformality [30,31]. However, due to the fact that the UV zero in the beta function stems from a logarithmic singularity and more generally that the beta function is not a physical quantity, careful consistency checks of these results are crucial [32].
For these reasons, and because it is important to uncover the phase diagram of four-dimensional quantum field theories, we initiate here a consistent lattice investigation of the ultraviolet fate of non-Abelian gauge-fermion theories at a small number of colors but at a large number of flavors where asymptotic freedom is lost. These parameters are currently inaccessible with other methods. Specifically, we will investigate the SU(2) gauge theory with 24 and 48 massless Dirac flavors. These two numerical values are substantially larger than the value where asymptotic freedom is lost, which is 11, and they are also very roughly estimated to be close to the region where the 1=N f squared corrections become relevant [29,31]. This was taken to be a sign of where one would expect the UV fixed point to disappear [31] and the theory develop a cutoff.
To numerically investigate the UV properties of these two theories, in our lattice analysis, we focus on the determination of the renormalized running gauge coupling. This is achieved by implementing the Yang-Mills gradient flow method at finite volume and with Dirichlet boundary conditions [33], enabling us to perform lattice simulations at vanishing fermion mass. We are able to demonstrate that our simulations are performed in the physical region connected to the Gaussian fixed point. This is corroborated by the fact that our results match onto the schemeindependent two-loop perturbative beta function in the infrared. We also discover that it is difficult to achieve large values of the renormalized coupling on the lattice. This might be interpreted as an early sign of the existence of a physical cutoff both for 24 and 48 flavors and that a larger number of flavors would be needed to achieve the safe fixed point. However, a closer look at the large-N f predicted UV behavior of these theories, if applicable, would suggest a safe fixed point to occur at a much stronger value of the gauge coupling 1 than achievable with the current simulations. This suggests a more conservative interpretation, that we are not yet able to explore the deep ultraviolet region where safety might emerge. Nevertheless, we feel that our work constitutes a necessary stepping stone toward unveiling the ultraviolet fate of nonasymptotically free gauge theories.
The paper is organized as follows. In Sec. II, we briefly review the analytical results for (non-)Abelian gauge 1 We observe that the leading N f beta function is scheme independent. theories at leading order in 1=N f . We compare the large-N f beta function with two-loop perturbation theory and discuss the conformal window 2.0 [31]. It is straightforward to specialize the results of this section to the case of SU (2) gauge theory investigated on the lattice. Section III is constituted by several subsections with the goal to make it easier for the reader to focus on the relevant aspects of the lattice setup and results. We begin with an introduction to the lattice action its features. We then discuss how the gauge coupling and its running is defined through the Yang-Mills gradient flow method. Then, we summarize the lattice results for N f ¼ 24 and N f ¼ 48. In Sec. IV, we offer our conclusions and directions for further studies.

II. CONFORMAL WINDOW 2.0: REVIEW OF THE ANALYTIC RESULTS
Consider an SUðN c Þ gauge theory with N f fermions transforming according to a given representation of the gauge group. Assume that asymptotic freedom is lost, meaning that the number of flavors is larger than N AF f ¼ 11C G =ð4T R Þ, where the first coefficient of the beta function changes sign. In the fundamental representation, the relevant group theory coefficients are At the one-loop order, the theory is free in the infrared, i.e., non-Abelian QED, and simultaneously trivial. As discussed in the Introduction, this means that the theory has a sensible continuum limit, by sending the Landau pole induced cutoff to infinity, only if the theory becomes noninteracting at all energy scales.
At two loops, in a pioneering work, Caswell [26] demonstrated that the sign of the second coefficient of the gauge beta function is such that an UV interacting fixed point, which would imply asymptotic safety, cannot arise when the number of flavors is just above the value for which asymptotic freedom is lost. This observation immediately implies that for gauge-fermion theories triviality can be replaced by safety only above a new critical number of flavors. In order to investigate this possibility, consider the large-N f limit at fixed number of colors. The leading order large-N f beta functions for QED and non-Abelian gauge theories were constructed in Refs. [27,28], while a summary of the results and possible investigation for the existence of UV fixed points appeared first in Refs. [29,30] with the scaling exponents computed first in Ref. [12]. Although in this work we will concentrate on non-Abelian gauge theories, we now briefly comment on the status of QED. Even though the large-N f beta function develops a nontrivial zero, it was demonstrated in Ref. [31] that at the alleged UV fixed point the fermion mass anomalous dimension violates the unitarity bound and hence the UV fixed point is unphysical. At this order in 1=N f , we conclude that QED is trivial.
For the non-Abelian case, using the conventions of Refs. [27,29], the standard beta function reads with g the gauge coupling. At large N f , it is convenient to work in terms of the normalized coupling A ≡ N f T R α=π.
Expanding in 1=N f , we can write 3 2 where the identity term corresponds to the one-loop result and constitutes the zeroth order term in the 1=N f expansion. If the functions jH i ðAÞj were finite, then in the large-N f limit, the zeroth order term would prevail, and the Landau pole would be inevitable. This, however, is not the case due to the occurrence of a divergences in the H i ðAÞ functions. According to the large-N f limit, each function H i ðAÞ resums an infinite set of Feynman diagrams at the same order in N f with A kept fixed. To make this point explicit, consider the leading H 1 ðAÞ term. The nth-loop beta function coefficients b n for n ≥ 2 are polynomials of order The coefficient with the highest power of T R N f will be b n;n−1 , and this is the coefficient contributing to H 1 ðAÞ at the nth-loop order. Moreover, it was shown in Ref. [34] that the b n;n−1 terms are invariant under the scheme transformations that are independent of N f (as appropriate for the large-N f limit). Now, the nth-loop beta function will have an interacting UV fixed point (UVFP) when the following equation has a physical zero [30]: This expression simplifies at large N f . Truncating at a given perturbative order n max , one finds that the highestloop beta function coefficient b n max contains just the highest power of ðT R N f Þ n max −1 multiplied by the coefficient b n max ;n max −1 , as can be seen from Eq. (3). Since this highest power of ðT R N f Þ n max −1 dominates in the N f → ∞ limit and since in this limit b 1 < 0, the criterion for the existence of an UV zero in the n max -loop beta function becomes [30] for N f → ∞; βðαÞ has an UVFP only if b n max ;n max −1 > 0: In perturbation theory, only the first few coefficients b n;n−1 are known, but, remarkably, it is possible to resum the perturbative infinite sum to obtain H 1 ðAÞ. From the results in Refs. [27,28], By inspecting I 1 ðxÞ and I 2 ðxÞ, one notices that the C G term in I 2 has a pole in the integrand at x ¼ 1 (A ¼ 3). This corresponds to a logarithmic singularity in H 1 ðAÞ, and will cause the beta function to have an UV zero already at this order in the 1=N f expansion, and by solving 1 þ H 1 ðAÞ=N f ¼ 0, this nontrivial UV fixed point occurs at [12] where k ¼ 16T R and l ¼ 18 Performing a Taylor expansion of the integrand in Eq. (5) and integrating term by term, we can obtain the nth-loop coefficients b n;n−1 and check our criteria above for the existence of the safe fixed point [34,35]. A preliminary investigation was performed in Ref. [34] up to 18th-loop order, where it was also checked that the first four loops agree with the known perturbative results. It was found that, even though up to the 12th-loop order the resulting coefficients are scattered between the positive and negative values, starting from the 13th-loop order all b n;n−1 are positive for the fundamental representation, two-index representations, and symmetric or antisymmetric rank-3 tensors. Unfortunately, the positivity of the coefficients is insufficient to prove the stability of the series and determine its radius of convergence. The first complete study of the analytic properties of the leading nontrivial large-N f expansion appeared recently in Ref. [35]. Here, it was demonstrated that an analysis of the expansion coefficients to roughly 30 orders is required to establish the radius of convergence accurately, and to characterize the (logarithmic) nature of the first beta function singularity.
These studies agree with the existence of a singular structure at leading order in 1=N f leading to a zero in the beta function. Although not a proof, see e.g., Ref. [36], it can be viewed as lending support for the possible existence of an UV fixed point in these theories. These results have been confirmed when extended to theories with Yukawa interactions [37][38][39] and employed to build realistic asymptotically safe extensions of the SM [15,17,18,40,41].
Using the results above, we can sketch a complete phase diagram, as a function of N c and N f , for an SUðN c Þ gauge theory with fermionic matter in a given representation. A robust feature of this phase diagram is the line where asymptotic freedom is lost, i.e., N AF f ¼ 11C G =ð4T R Þ. As is well known, decreasing N f slightly below this value, one achieves the perturbative Banks-Zaks infrared fixed point (IRFP), that at two loops yields α Ã ¼ −b 1 =b 2 . This analysis has been extended to the maximum known order in perturbation theory in Refs. [30,42,43].
Just above the loss of asymptotic freedom, as already mentioned, Caswell [26] demonstrated that no perturbative UVFP can emerge. By continuity, there should be a region in color-flavor space where the resulting theory is non-Abelian QED with an unavoidable Landau pole. This is the unsafe QCD region. The theories in this region are low energy effective field theories featuring a trivial IRFP. This means that one can expect the existence of a critical value of number of flavors N Safe f above which safety emerges. This region extends to infinite values of N f , i.e., the safe QCD region [31].
For the fundamental representation, the leading 1=N f expansion is applicable only for N c ≲ N f =10, while for the adjoint representation, we find N f ≳ 7 for any N c . Following Ref. [31], it is sensible to use these values as a first rough estimate of the lower boundary of the safe QCD region. Altogether, these constraints allowed the authors to draw the corresponding phase diagrams in Ref. [31]. For the reader's convenience, we draw in Fig. 1 again the phase diagrams presented first in Ref. [31] both for the fundamental (panel a) and adjoint representations (panel b).
Before specializing to the theories that we will investigate on the lattice, let us comment also on the safe status of supersymmetric gauge theories. An UV safe fixed point can, in principle, flow to either a Gaussian IR fixed point (noninteracting) or to an interacting IR fixed point. So far, for the nonsupersymmetric case, we discussed the first class of theories because it is theoretically and phenomenologically important to assess whether nonasymptotically free theories can be UV complete, up to gravity. We provided an affirmative answer for gauge-Yukawa theories that are remarkably similar in structure to the Standard Model in Ref. [12]. The situation for gauge fermions is more involved, and this is the reason we further investigate it here via first-principle lattice simulations. Given the above, the general conditions that must be (nonperturbatively) abided by nonasymptotically free supersymmetric theories to achieve safety were put forward in Ref. [11] generalizing and correcting previous results of Ref. [63]. To make the story short, at least one chiral superfield must achieve a large-R charge at the safe fixed point to ensure that the variation of the a-function between the safe and Gaussian fixed points is positive as better elucidated in Ref. [64]. Models of this type were shown to exist in Refs. [63][64][65].
Another possible way to elude the constraints [11,64] is to consider UV fixed points flowing to IR interacting fixed points. Within perturbation theory, nonsupersymmetric theories of this type were discovered in Ref. [66], and supersymmetric theories were discovered in Refs. [65,67]. Summarizing the status for supersymmetric safety, we can say that theories abiding the constraints of Ref. [11] exist. Nevertheless, specializing Ref. [11] to the case in which asymptotic freedom is lost (N f ≥ 3N c ) for super QCD with(out) a meson [11,63], one can show that these theories are unsafe for any number of matter fields. This is in agreement with the 1=N f studies of the supersymmetric beta functions [63,68] as explained in Ref. [35]. The unsafety of super QCD with(out) a meson should be contrasted with QCD at large number of flavors for which, as we argued above, safety is possible [31] and QCD with a meson for which safety is a fact [12].
After the supersymmetric parentheses, it is time to resume the investigation of the nonsupersymmetric conformal window 2.0. What has been reviewed so far clearly motivates a nonperturbative study of the large number of flavor dynamics via lattice simulations. Recalling that the rough estimate for the lower boundary of the safe side of the conformal window 2.0 requires N f > 10N c , it is clear that to minimize the number of flavors we should use N c ¼ 2, with N f > 20. Therefore, in Fig. 2, we present the leading N f beta function for an SU(2) gauge theory with either 24 or 48 flavors. The large-N f beta function is shown by solid curves, while the dotted curves show the five-loop perturbative beta function. The upper panel is for the N f ¼ 24 case, and the lower for the N f ¼ 48 one. At leading order in 1=N f , the beta functions support the presence of an UVFP for both flavors. It is instructive to show also the four-and five-loop perturbative results. These are the two theories we choose to investigate on the lattice, i.e., an SU(2) gauge theory with either 24 or 48 flavors transforming according to the fundamental representation of the underlying gauge group.

A. Lattice formulation
In this section, we will define the model we consider and discuss the methods we use. Our treatment of the general features will be brief here as more detailed description can be found in Refs. [49,69]. The model is defined by the lattice action where U is the standard SU(2) gauge link matrix, V is smeared gauge link defined by hypercubic truncated stout smearing (HEX smearing) [70], S G is the single plaquette gauge action, and S F is the Wilson fermion action with the clover term S SW . We set the Sheikholeslami-Wohlert coefficient c SW ¼ 1, which is a good approximation with HEX smeared fermions [69]. The coupling constant is measured using the gradient flow method with Dirichlet boundary conditions [33]. This method has been used successfully to measure the evolution of the coupling constant in SU(2) gauge theory with N f < 10, motivating its use also here [49,69]. On a lattice of size L 4 , we use periodic boundary conditions at the spatial boundaries. At the temporal boundaries x 0 ¼ 0, L, we use Dirichlet boundary conditions by setting the gauge link matrices to U ¼ V ¼ 1 and the fermion fields to zero. These boundary conditions enable simulations at vanishing fermion mass. We run the simulations using the hybrid Monte Carlo algorithm with the fifth order Omelyan integrator [71,72] and chronological initial values for the fermion matrix inversions [73]. The step length is tuned to have an acceptance rate of the order of 80% or higher. We run the simulations with bare couplings varying within the range β L ≡ 4=g 2 0 ∈ ½−1; 6, and for each value of β L , we tune the hopping parameter κ to its critical value κ c ðβ L Þ, for which the absolute value of the PCAC fermion mass [74] is of the order of 10 −4 on a lattice of size 24 4 . The obtained critical hopping parameter values, κ c ðβ L Þ, are then used for all used lattice sizes L=a ¼ 12, 18, 24, and 30. A summary of the simulation parameters and corresponding PCAC quark masses, as well as the acceptance rates and accumulated statistics for each simulation, is given in the Tables I and II. Because we need to use strong coupling in the UV, we are forced to use small values for β L . We include even negative values for β L in our study. This is to compensate for the effective positive β L shift induced by the Wilson fermions [75,76]. To leading order, this shift is proportional to the number of flavors and can therefore be substantial at large N f . As a consequence, even at the smallest β L ¼ −1 at N f ¼ 48, the lattice gauge field observables (for example, the plaquette) behave as if the effective gauge coupling remains positive. A qualitatively similar effect has been observed with staggered fermions [10]. We note that SU(2) pure gauge lattice theory has a smooth bulk crossover at g 2 0 ∼ 1.7-1.8, which we do not observe in the large-N f case within the range of couplings used.
To define the running coupling, we apply the Yang-Mills gradient flow method [33,77,78]. This method defines a flow that smooths the gauge fields, removes UV divergences, and automatically renormalizes gauge invariant objects [79]. 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 where A μ is the original continuum gauge field. In the lattice formulation, the lattice link variable U replaces the continuum flow field, which is then evolved using the tree-level improved Lüscher-Weisz pure gauge action [80].

B. What to expect
The gradient flow equation has the form of a heat equation, which, in the case of a free gluon field in the Landau gauge, is for example solved by the heat kernel [81]  Hðy − x; tÞ ¼ 1 Hence, by evolving a gauge configuration with the gradient flow for some time t, one performs a Gaussian smearing with smearing radius σ ¼ ffiffiffiffi 2t p , which means that physical processes at length scales below 2σ ¼ ffiffiffiffi 8t p ≕ λðtÞ get effectively integrated out. The scale μðtÞ ¼ 1=λðtÞ can therefore be interpreted as renormalization scale of gauge observables measured at flow time t.     The nonperturbative renormalized coupling at scale μ [81] is then defined via energy measurement as where N is a normalization factor that has been calculated in Ref. [82] to match the MS coupling at tree level and the gauge field energy density E is measured only on the central time slice of the lattice: x 0 ¼ L=2.
On the lattice, the finite lattice spacing a (which is a function of the bare lattice parameters β L and κ) and the finite system size L ¼ Na (with N being the number of lattice sites in each direction) restrict the renormalization scales, accessible with the gradient flow, to the range 1=L < μðtÞ < 1=a. Per construction, the flow always starts on the UV side of this interval and evolves the gauge field toward the IR. At the lattice scale 1=a, the lattice theory deviates strongly from its continuum counterpart, and the gradient flow smearing scale λðtÞ ¼ ffiffiffiffi 8t p should therefore reach at least two to three lattice spacings before the renormalized coupling (or any other observable of the lattice gauge field) at flow time t can be expected to behave like the corresponding continuum quantity.
One can now ask how the lattice gradient flow coupling g 2 GF ðλðtÞ; β L Þ should behave, as function of the flow scale λðtÞ, if the theory possesses an interacting UV fixed point. Before addressing this question for the lattice, let us recall how the running continuum coupling g 2 ðλ=λ 0 ; g 2 0 Þ, i.e., the solution to the differential equation (1) with μ ¼ 1=λ, behaves as function of increasing λ=λ 0 > 0 for different choices of the initial condition or reference coupling g 2 0 ¼ gð1; g 2 0 Þ at reference scale λ 0 , if there is an interacting UV fixed point at coupling g 2 ¼ g 2 cr . The behavior is illustrated schematically in Fig. 3: the fixed point is unstable, and therefore, for increasing λ=λ 0 > 1, the running coupling gðλ=λ 0 ; g 2 0 Þ will: (a) decrease, if g 2 0 < g 2 cr ; (b) remain constant, if g 2 0 ¼ g 2 cr ; and (c) increase, if g 2 0 > g 2 cr . The qualitatively distinct behavior of the running coupling on the two sides of the UV fixed point implies a drastic change in the long-distance behavior of the system when the reference coupling g 2 0 crosses the critical value g 2 cr . We note that if g 2 0 is much smaller than g 2 cr the behavior of the running coupling for g 2 ðλ=λ 0 ; g 2 0 Þ < g 2 cr can be almost indistinguishable from a Landau pole type behavior. The coupling constant evolution curves demonstrating an UV fixed point shown in Fig. 3 have been obtained by integrating the perturbative four-loop beta function at N f ¼ 48, which indeed has a zero at g 2 cr ≈ 6.8 (dashed curve in lower panels of Fig. 2). The Landau pole curve has been obtained from the corresponding one-loop beta function. While we naturally cannot expect purely perturbative results to be reliable, these curves give us qualitatively plausible scenarios for the coupling constant evolution. We can conclude that it can be difficult to distinguish between UV fixed point and Landau pole behaviors by looking at the behavior of the running coupling for g 2 0 ≪ g 2 cr . Switching now to the lattice, we can think of the lattice scale a as defining a reference length scale λ 0 for the corresponding continuum theory. The gradient flow scale λðtÞ, given in lattice units, then corresponds to the ratio of scales λ=λ 0 , and we can for λ=λ 0 ≫ 2 identify the running continuum coupling g 2 ðλ=λ 0 ; g 2 0 Þ with the lattice gradient flow coupling g 2 GF ðλðtÞ; β L Þ. As already mentioned above, for λ=λ 0 ≤ 2, discretization effects are strong. Therefore, the lattice gradient flow coupling cannot be expected to behave like the running coupling of the continuum theory in this regime, 2 which means that the relation between the bare inverse lattice coupling β L and the corresponding continuum coupling g 2 0 ¼ gðλ=λ 0 ¼ 1; g 2 0 Þ cannot be read off from the value of g 2 GF ðλðtÞ; β L Þ at flow scale λðtÞ ¼ a. In order to verify the existence of an interacting UV fixed point via lattice simulations, one has to find a value of β L , for which a ¼ λ 0 is sufficiently small, so that the corresponding g 2 0 is larger than g 2 cr . One should then observe, that the gradient flow coupling increases with increasing flow  3. A cartoon of the running coupling g 2 ðλ=λ 0 ; g 2 0 Þ as a function of increasing (relative) renormalization length scale λ=λ 0 in the presence of an interacting UV fixed point, located at g 2 ¼ g 2 cr , for different choices of the initial coupling g 2 0 ¼ g 2 ð1; g 2 0 Þ and corresponding reference scale λ 0 . Given the fixed point at g 2 cr , the running coupling gðλ=λ 0 ; g 2 0 Þ as a function of increasing λ=λ 0 ≥ 1 will decrease if g 2 0 < g 2 cr , remain constant if g 2 0 ¼ g 2 cr , and increase if g 2 0 > g 2 cr . The fixed point is repelling, so that one can only flow away from it. If g 2 0 is much smaller than g 2 cr , the UVFP and Landau pole running couplings are almost indistinguishable for g 2 ðλ=λ 0 ; g 2 0 Þ < g 2 cr . 2 In fact, as we use in our simulations the clover Wilson fermion action with HEX smeared gauge links (meaning that the elementary cells on which the action is defined have a linear size of at least 2a instead of just a), g 2 GF ðλðtÞ; β L Þ and g 2 ðλ=λ 0 ; g 2 0 Þ will even start to agree only for λðtÞ ¼ λ=λ 0 ≳ 4. scale λðtÞ, whereas for g 2 0 < g 2 cr , one will always flow toward the trivial IR fixed point. In addition, one should also observe a phase transition associated to the aforementioned change of the long-distance behavior of the system when g 2 0 crosses g 2 cr . For a theory which is IRinstead of UV-free, the dependency of the lattice spacing a on the bare lattice parameter β L is unfortunately not necessarily such that a can become arbitrarily small, and it might therefore not be possible for a ¼ λ 0 to reach values for which g 2 0 ≥ g 2 cr . The discovery of an interacting UV fixed point on the lattice in this way is possible only if the true, nonperturbative beta function is such that the discussion in connection with Fig. 3 applies. As discussed above, the four-loop beta function from Fig. 2 was used to model the UV fixed point behavior in Fig. 3. This four-loop beta function is smooth everywhere, and its slope in the neighborhood of the UV fixed point is moderate. However, if the nonperturbative beta function behaves more like the leading order large-N f beta function (red curve in lower panels of Fig. 2), the situation is significantly more complicated, as the UV fixed point arises due to a singularity as β → −∞. This abrupt change makes the running coupling for g 2 0 < g 2 cr look even more as if there is a Landau pole at λ=λ 0 ≤ 1 (cf. Fig. 4), and because of the singularity right after the UV fixed point, the behavior of the running coupling for g 2 0 > g 2 cr is unknown. The slope of the beta function at the fixed point was computed in Ref. [12] and linked to the anomalous dimension of glueballs in Ref. [32]. The divergent slope can either lead to violation of unitarity or be interpreted as the physical decoupling of glueballs at the fixed point if the fixed point persists. A glueball decoupling is in line with the physical mechanism for the existence of a fixed point which is enforced by the fermion dynamics. A more conservative viewpoint is that the consistent inclusion of subleading 1=N f terms will smooth out the singular structure. It is worth noting that the comparison with the large-N f beta function is relevant also away from the fixed point since it contains all the leading N f orders in α. As we explained earlier the series in αN f is convergent, and the coefficients are scheme independent. One should therefore observe a better agreement of the large-N f beta function with 48 flavors than with 24 away from the singularity. These figures show that the measurements of g 2 GF are dominated by lattice artifacts both at too small and at too large flow scales, λ ≲ 3a and λ ≳ 0.3L, respectively. When the flow scale is of the order of lattice spacing or smaller, large artifacts can be expected, and indeed, from the definition of the gradient flow coupling, Eq. (10), we see that g 2 GF → 0 as λ → 0 (t → 0). This causes the characteristic "peak" at small λ in Figs. 5 and 6, as g 2 GF develops toward more continuumlike behavior. The use of the HEX smeared clover fermions also increases the range of interaction terms of the lattice action to approximately 3-4a, which may also affect the flow at small λ.
Apart from these UV cutoff effects, the gradient flow is also affected by IR cutoff effects due to the finite lattice size. In our case, these seem to be particularly strong, as can be seen by comparing the curves for different system sizes in the individual panels of Figs. 5 and 6: the flow scales at which the g 2 GF for different system sizes start do deviate marks the scale at which IR cutoff effects start to dominate in the smaller of the two systems. This seems here to happen already at scales around λ ≳ 0.3L almost independently of the bare inverse lattice coupling β L .
Due to the strong finite size and finite volume effects, the range of flow scales λ, for which the lattice gradient flow FIG. 4. The running coupling g 2 ðλ=λ 0 ; g 2 0 Þ for the leading order large-N f beta function from Eq. (2) with N f ¼ 48 (cf. lower panels of Fig. 2) as a function of increasing (relative) renormalizaion length scale λ=λ 0 for different choices of initial coupling g 2 0 ¼ g 2 ð1; g 2 0 Þ and corresponding reference scale λ 0 . The red, horizontal dotted line indicates the location of the UV fixed point, In contrast to the four-loop N f ¼ 48 beta function used for the discussion in Fig. 3, which had a moderate slope at the fixed point, the slope of the leading order large-N f beta function at the fixed point g 2 cr is very steep. This makes it almost impossible to draw the line of constant running coupling, that would start at g 2 0 ¼ g 2 cr and remain there. Consequently, the running coupling for g 2 0 < g 2 cr looks even more as if there is a Landau pole at λ=λ 0 ≤ 1, as can be seen by comparing for example the g 2 0 ¼ 1 curve with the corresponding one-loop result (dashed line), which has a Landau pole at λ=λ 0 ¼ s lp ≈ 0.041. The leading order large-N f running coupling matches the oneloop curve almost perfectly up to g 2 cr , where it suddenly stops. As this beta function ends right after the UV fixed point in a singularity at −∞, we do not have any information about the behavior of the running coupling above g 2 cr .
coupling g 2 GF ðλ; β L Þ can be expected to behave like the running coupling g 2 ðλ=λ 0 ; g 2 0 Þ of the continuum theory, is approximately given by λ ∈ ½3a; 0.25L, where the lower bound of 3a is somewhat optimistic. This renders the smallest lattices L ≤ 18 of very limited use.
We can optimize the gradient flow coupling by adding a shift τ to the flow time [83]: The effect of τ vanishes in the continuum limit (if it exists), and it can be tuned to get rid of most Oða 2 Þ lattice artifacts. Due to the very large finite size artifacts, we cannot use the method used in Refs. [49,83] to optimize τ. However, as it turns out that the gradient flow coupling obtained by our simulations are quite small, we can expect the scheme independent perturbative two-loop beta function to be a very accurate approximation to the true beta function at large enough flow time (cf. Fig. 2). For each β L , we tune τ by matching the largest volume L ¼ 30a gradient flow coupling to the two-loop perturbative coupling over the interval λ ∈ ½3a; 5a. We then use this value for τ for smaller volumes L=a ¼ 12, 18, 24. Examples for (11) at optimal τðβ L Þ are shown in Fig. 7 (N f ¼ 24) and Fig. 8  (N f ¼ 48) for two different values of β L , together with the corresponding two-loop running coupling. In Fig. 9, we show the step scaling function Σðu; s; L; cÞ ¼ g 2 GF;sL ðcsL; β L Þj β L ∶g 2 GF;L ðcL;β L Þ¼u ; ð12Þ which tells us how the coupling constant evolves when the length scale where it is evaluated changes from cL to scL.
In the following, we use c ¼ 0.22; i.e., the gradient flow time is fixed so that ffiffiffiffi 8t p ¼ 0.22L. We are forced to use a relatively small value for c (in comparison with conventional c ¼ 0.4…0.5, used in Ref. [49], for example) in order to avoid excessive finite volume effects. The step scaling is particularly well suited for doing a continuum extrapolation (if it exists) of the running coupling. For a conventional extrapolation, one would, however, need a set of lattice sizes A ¼ fL 1 ; …; L n g, from which one can form several pairs ðL 0 ; LÞ that correspond to the same ratio L 0 =L ¼ s. The lattice volumes available to us do not allow this. Furthermore, very large lattice artifacts at smaller volumes would make the extrapolation unreliable.
What we can do instead is compare the step scaling function for different s directly with the corresponding twoloop results in Fig. 9. The agreement is reasonable, except for s ¼ 1=2, which involves the smallest system size L ¼ 12 for which the finite size and finite volume effects are particularly strong.
A more transparent comparison of data from pairs of system size ðL 0 ; LÞ with different s ¼ L 0 =L can be obtained with the discrete beta function β ðsÞ L ðu;L;cÞ ¼ − which approaches the conventional beta function in the limits s → 0 or g 2 → 0. Because g 2 GF is small, this is a good approximation of the true beta function. The results are shown in Fig. 10 for c ¼ 0.22 and s ¼ 3=2 ¼ 18=12, s ¼ 2 ¼ 24=12, and s ¼ 30=18 ¼ 5=3. For comparison, the corresponding results for the discrete beta function, evaluated with the two-loop running coupling, are shown as well. As we can observe, the lattice result approaches the two-loop result as the volume increases.
From the results above, we see that the measured gradient flow couplings are very small in the region where we can trust the measurements, at most g 2 GF < 1.75. The couplings grow as the distance is reduced (λ decreased), until lattice artifacts take over at λ ≲ 3a. In order to reach strong couplings, we use small values for the inverse bare coupling β L (even negative, as discussed in Sec. III A). However, it turns out that if β L is small enough, the gradient flow coupling becomes almost independent of its value. This behavior is compatible with the Landau-pole-like behavior, as we will discuss in Sec. III C 1 below. Fig. 7 but for N f ¼ 48. 1. Effective g 2 0 I As we saw above, the gradient flow coupling makes sense only when the flow scale is large enough, i.e. larger than about 3a. Nevertheless, it is of interest to relate g 2 GF to some lattice ultraviolet scale effective coupling. The inverse of β L ¼ 4=g 2 0 does not make much sense, because it can be negative.

FIG. 8. As in
Here, we use the value of the plaquette (trace of the 1 × 1 Wilson loop) as a proxy; it is an UV quantity on the lattice, and its expectation value is readily measurable. To convert the plaquette values to effective couplings, we use an "inverse Monte Carlo" procedure: we perform new simulations using a pure gauge SU(2) lattice theory with Wilson plaquette action [i.e., using only S G ðUÞ in Eq. (7)] and tune the pure gauge inverse coupling β PG L so that the plaquette expectation value matches that of the original theory. The effective bare coupling is obtained by inversion This quantity is positive for all of our lattices. Because the plaquette expectation value is, in practice, independent of the volume, g 2 0;eff is only a function of β L and N f . It increases monotonically as β L is decreased. At our smallest β L values, −0.3 (N f ¼ 24) and −1.0 (N f ¼ 48), it reaches maximum values of about 6 and 5.3, respectively.
In Fig. 11, we plot the gradient flow couplings g 2 GF , measured at fixed flow scale λ ¼ 4a, against g 2 0;eff . The gradient flow coupling is almost independent of the volume, excluding L=a ¼ 12, where flow scale 4a is large enough so that finite volume effects become important. We shall ignore L=a ¼ 12 in the following.
As mentioned earlier, the pure gauge system undergoes a smooth crossover to a lattice bulk phase around g 2 0;eff ≈ 1.7-1.8. This artifact, not present in our simulations with fermions, does not invalidate the use of g 2 0;eff as proxy for the coupling at the lattice cutoff scale, but the functional dependency of the gradient flow running coupling on g 2 0;eff can start to deviate from the expected behavior when g 2 0;eff enters the bulk region, i.e., when g 2 0;eff > 1.7-1.8. This is shown by vertical lines in Fig. 11. Nevertheless, the excellent match of the lattice measurements with the perturbative curve indicates that the effect is small.
We observe the following behavior: as g 2 0;eff → 0, the ratio g 2 GF =g 2 0;eff → 1, as dictated by universality at weak coupling. On the other hand, as g 2 0;eff grows, g 2 GF seems to approach a constant. This is characteristic behavior for an UV Landau pole; as the UV scale approaches the Landau pole, the UV scale coupling diverges. Because g 2 GF is evaluated at length scales which are by a constant factor larger than the UV length scale, the value of g 2 GF will instead approach a constant value.
The couplings g 2 0;eff and g 2 GF have been obtained using different schemas, and thus their values cannot be compared without matching. Nevertheless, we obtain a reasonable fit of the data in Fig. 11 with the two-loop perturbative running coupling g 2 ðλ=λ 0 ; g 2 0 Þ by identifying g 2 0 ≡ g 2 0;eff and g 2 ðλ=λ 0 ; g 2 0 Þ ≡ g 2 GF ðλ=a; g 2 0;eff Þ. A good qualitative agreement is obtained with the matching coefficients λ 0 ¼ a=9 for N f ¼ 24 and λ 0 ¼ a=18 for N f ¼ 48. This implies that the effective lattice coupling g 2 0;eff corresponds to a reference scale, which is about 9, respectively 18, times smaller than the lattice scale a. The fit is compatible with the existence of the Landau pole, because the two-loop beta function also features one. However, we naturally cannot exclude the existence of an UVFP at stronger UV coupling than reached here.
The possible existence of the Landau pole gives a natural explanation for the small value of the gradient flow coupling, the lattice UV length scale is approximately a, whereas the gradient flow coupling g 2 GF is evaluated at length scale 3a or larger. Thus, in terms of energy, g 2 GF is evaluated at scale μ GF < μ LP =3, where μ LP is the scale where Landau pole appears. This gives ample room for the coupling to decrease. The actual value of the coupling depends on the details of the scheme.
2. Effective g 2 0 II Another possibility to define an effective coupling for the lattice theory comes as a side product of the method we used to determine the optimal value of the τ parameter in the definition (11) of the improved gradient flow coupling, i.e., FIG. 11. The relation between the gradient flow coupling g 2 GF ðL; β L ; cÞ at flow scale λ ¼ cL ¼ 4a and the effective bare coupling g 2 0;eff ðβ L Þ at N f ¼ 24 (left) and N f ¼ 48 (right). The dashed line is the perturbative (pert.) two-loop result, where the continuum couplings g 2 0 and g 2 ðλ=λ 0 ; g 2 0 Þ are, respectively, identified with the g 2 0;eff and g 2 GF ðλ=a; g 2 0;eff Þ from the lattice. The g 2 0;eff ¼ g 2 GF line is shown with dots, and the vertical dashed line indicates the value of g 2 0;eff above which the pure gauge system is in the lattice bulk phase.
After having determined in this way the optimal τ, one also has fixed the values ðλ 0 ; g 2 0 Þ ¼ ðλ 0 ðτÞ; g 2 GF;0 ðτÞÞ, for which (18) has within the given fit interval the best overlap with the τ-shifted lattice data. Solving now λ pt ðg 2 ; g 2 0 ; λ 0 Þ ¼ 1 ð20Þ for g 2 , one obtains an effective coupling for the lattice theory at the cutoff scale a, which we call g 2 0;eff2 ðβ L Þ. In Figs. 7 and 8, the value of g 2 0;eff2 ðβ L Þ can be read off from the value of the perturbative two-loop curve (dashed line) at λ ¼ 1.
From the derivation of the effective coupling g 2 0;eff2 ðβ L Þ, it is not a big surprise that one finds in Fig. 12 excellent agreement between lattice data and two-loop result, when plotting the gradient flow coupling at fixed flow scale λ (in lattice units) as function of g 2 0;eff2 ðβ L Þ, identifying the perturbative g 2 ðλ=λ 0 ; g 2 0 Þ and g 2 GF ðλ; β L Þ, with λ 0 ≡ a and g 2 0 ≡ g 2 0;eff2 ðβ L Þ.

IV. CONCLUSIONS
The UV behavior of gauge-fermion theories as a function of the number of matter fields poses a fundamental problem on our understanding of quantum field theory, both perturbatively and nonperturbatively. Recent theoretical developments suggest that at inifinite number of flavors an interacting UV fixed point could exist making these theories asymptotically safe.
In this paper, we have discussed in detail the categorization of gauge-fermion theories based on their possible UV behaviors. Then, we described our computational setting to investigate the renormalization group evolution of SU(2) gauge theory with 24 or 48 Dirac fermions. Finally, we presented the results of our first pioneering analysis.
We found that with the methodologies developed in our past work, the nonperturbative computations can be successfully carried out. We have demonstrated that our results match well with perturbation theory. On the other hand, we were unable to reach lattices where strong renormalized couplings could be controllably obtained. We used the gradient flow procedure to determine the renormalized coupling on the lattice. By construction, gradient flow coupling can be evaluated at length scales which are at least a few lattice spacings. This leaves room for a Landau pole or otherwise large couplings to exist at shorter distances (in the lattice ultraviolet limit). We observed indications of this behavior by defining an effective lattice ultraviolet coupling. The observed behavior is compatible with the existence of a Landau pole, but we cannot exclude an FIG. 12. The relation between g 2 GF and effective lattice coupling g 2 0;eff2 ðβ L Þ at N f ¼ 24 (left) and N f ¼ 48 (right), evaluated at flow scale λ ¼ cL ¼ 4a (upper panels) and 6a (lower panels). The dashed line is the perturbative (pert.) two-loop result, where the continuum couplings g 2 0 and g 2 ðλ=λ 0 ; g 2 0 Þ are, respectively, identified with the g 2 0;eff and g 2 GF ðλ=a; g 2 0;eff2 Þ from the lattice. The g 2 0;eff2 ¼ g 2 GF line is shown with dots.
ultraviolet fixed point at even stronger couplings than reached here. Consequently, the true continuum behavior of these theories in the deep ultraviolet remain undetermined. Nevertheless, our study provides an important milestone toward establishing the ultaviolet fate of gauge-fermion theories when asymptotic freedom is lost. The results further motivate investigations of gauge theories at a very large number of fermion fields.