Continuous renormalization group β function from lattice simulations

Motivated by the connection between Wilsonian real-space renormalization group and the gradient flow transformation we present an approach to calculate the continuous renormalization group β function in non-perturbative lattice simulations. While our method requires an additional extrapolation compared to the traditional step scaling calculations, it has several advantages that compensates for this extra step. We illustrate our approach by calculating the β function of 2-flavor QCD and show that lattice predictions from individual lattice ensembles, even without the required continuum and finite volume extrapolations, predict the β function very close to the result of the full analysis. Thus our method provides a non-perturbative framework and intuitive understanding into the structure of strongly coupled system, in addition to being complementary to existing lattice calculation.


INTRODUCTION
The renormalization group (RG) β function encodes the energy dependence of the running coupling. While the β function is scheme dependent, the number of its zeros, corresponding to infrared and ultraviolet fixed points (IRFP, UVFP), as well as the slope around the zeros are universal. The characteristic structure of the β function distinguishes conformal vs. confining, asymptotically vs. infrared free, trivial vs. asymptotically safe systems [1][2][3][4][5][6][7]. The perturbative β function of 4-dimensional nonabelian gauge-fermion systems are known up to 5-loop level in the MS scheme, but the perturbative expansion is unreliable at strong couplings [5,[8][9][10]. The β function calculated non-perturbatively is essential to describe strongly coupled systems both for QCD and within the conformal window.
We demonstrate the method by calculating the β function of a QCD-like system with two flavors and SU(3) gauge group. The final prediction of this study with relatively low statistics is shown by the grey band in Fig. 1. It is consistent with existing finite volume step scaling function calculations of 3-flavor QCD [21] in that it is close to the 1-loop perturbative prediction. More interesting is that the colored data points in Fig. 1, corresponding to raw lattice data in finite volume, differ only slightly from the result of the full analysis. The continuous β function approach predicts the running of  the renormalized coupling in a transparent way where cut-off and finite volume effects are clearly identifiable. This property could be particularly helpful when analyzing near-conformal systems, the properties of which are less understood than QCD. GF measurements of existing configurations of step-scaling studies can be reanalyzed to predict the continuous β function without additional computational cost. Hence this method provides an alternative to test systematical errors.
We find it helpful to describe the lattice calculation of the continuous β function using the language of realspace Wilsonian RG transformation [32][33][34]. The connection between GF and RG is discussed in detail in Ref. [35]. Gradient flow is a continuous smearing transformation that is appropriate to define real-space RG blocked quantities, but it is not an RG transformation as it lacks the crucial step of coarse graining. However, coarse graining can be incorporated when calculating expectation values. In particular, expectation values of local operators, like the energy density that enters the definition of the GF coupling, are identical with or without coarse graining. When the dimensionless GF time t/a 2 is related to the RG scale change as b ∝ t/a 2 , the GF transformation describes a continuous real-space RG transformation.
The topology of RG flow on the chiral m = 0 critical surface in an asymptotically free gauge-fermion system is sketched in Fig. 2. g 1 represents the relevant gauge coupling at the Gaussian fixed point (GFP), while g 2 refers to all other irrelevant couplings. The GFP is on the g 1 = 0 (lattice spacing a = 0) surface and the renormalized trajectory (RT) emerging from the GFP describes the cut-off independent continuum limit at finite renormalized coupling. Numerical simulations are performed with an action characterized by a set of bare couplings. If this action is in the vicinity of the GFP or its RT, the typical RG flow approaches the RT and follows it as the energy scale is decreased from the cut-off towards the infrared as indicated by the blue lines. RG flows starting at different bare couplings approach the RT differently, but once the irrelevant couplings have died out, they all follow the same 1-dimensional renormalized trajectory and describe the same continuum physics. The RT of chirally broken systems continues to g 1 → ∞, while conformal systems have an IRFP on the RT that stops the flows from either direction. While the topology of the RG space is universal, the location of the fixed points and their corresponding RTs depend on the RG transformation.
The RT is a 1-dimensional line, therefore, a dimensionless (zero canonical and zero anomalous dimension) local operator with non-vanishing expectation value can be used to define a running coupling along the RT. The simplest such quantity in gauge-fermion systems is the energy density multiplied by b 4 (or t 2 ) to compensate for its canonical dimension. This is indeed the quantity defined in Ref. [13] as the gradient flow coupling g 2 GF (t; g 2 0 ) ∝ t 2 E(t) . E(t), the energy density at flow time t, can be estimated through various local lattice operators like the plaquette or clover operators. At large flow time irrelevant terms in the lattice definition of E(t) die out. In that limit g 2 GF approaches a continuum renormalized running coupling and its derivative is the RG β function The above definition is valid in infinite volume only. In a box of finite length L the RG equation contains the term L(dg 2 GF /dL), a difficult to estimate quantity. In our approach we extrapolate L/a → ∞ at fixed t/a 2 which also sets the renormalization scheme c = 0. The continuum limit of the β function is obtained at fixed g 2 GF while taking t/a 2 → ∞. In QCD-like systems this automatically forces the bare gauge coupling towards zero, the critical surface of the GFP.
The Wilsonian RG description suggests that lattice simulations at a single bare coupling can predict, up to controllable cut-off corrections, a finite part of the RG β function. In practice the finite lattice volume limits the range where the infinite volume β function is well approximated. Chaining together several bare coupling values, we can cover the entire RT while the overlap and deviation between different volume and bare coupling predictions characterizes the finite volume and finite cut-off effects as illustrated in Fig. 1.
Once the GF coupling is determined and its derivative is calculated as the function of the GF time, the continuous β function calculation requires two steps: A) Infinite volume extrapolation at every GF time.
B) Infinite flow time extrapolation at every g 2 GF .
Step B) removes cut-off effects and replaces the L/a → ∞ continuum limit extrapolation of the step-scaling function approach.
Step A) is new in the continuous β function approach but is compensated by several advantages. In all GF analysis the flow time has to be chosen large enough to remove all but the largest irrelevant operator even on the smallest volume considered. In traditional step-scaling calculations the flow time grows with L 2 which leads to large statistical errors on the largest volumes. In the continuous β function approach the flow time is independent of the volume. This significantly reduces statistical errors. Finally the continuum limit of the continuous β function is obtained by extrapolating a continuous function of the flow time. While the data are highly correlated, a continuous function nevertheless allows control to determine the functional form, e.g. the scaling exponent of the irrelevant operator. The correlations themselves can be handled by a fully correlated analysis similar to fitting subsequent time slices of a 2point correlation function. The RG β function is defined in the chiral limit. Finite fermion mass introduces a relevant operator that drives the RG flow away from the critical surface. Thus simulations with am = 0 are essential to avoid yet another am → 0 extrapolation. This is always possible in chirally symmetric conformal systems and can be enforced in chirally broken models by limiting the simulation volumes to be finite in physical units, below the inverse of the critical temperature. The same constraint is present in step-scaling studies.
A continuous β function based on GF around the GFP has been considered before. The only published result is a prediction at one renormalized gauge coupling [26] that assumes perturbative scaling. Motivated by Wilsonian RG transformation, our fully non-perturbative approach determines β(g 2 GF ) over a wide range in g 2 GF . It is equally valid near a conformal fixed point and can predict a nonperturbative scaling exponent. The consistency of direct lattice data at different bare couplings shown in Fig. 1 provides an intuitive understanding of the structure of the investigated system. In addition we set the 5th dimension of DW fermions to L s = 12 which leads to residual masses am res < 10 −6 for all gauge couplings considered. The same combination of actions has been used in recent works [25,29] and properties for QCD simulations are reported in [38][39][40][41]. Measurements for three different flows, Wilson (W), Symanzik (S), and Zeuthen (Z), have been implemented in Qlua [42,43] and we analyze three operators, Wilson plaquette (W), clover (C), and Symanzik (S) to estimate the energy density [44,45]. Our data analysis is performed using the Γ-method which is designed to estimate and account for autocorrelations [46].
The GF coupling is defined as The normalization ensures that g 2 GF matches the MS coupling at tree level, and the term 1/(1 + δ) corrects for the gauge zero modes due to periodic gauge boundary con-  ditions [14]. On L 3 × L t volumes where ϑ is the standard Jacobi elliptic function [14]. We calculate β(g 2 GF ) = −t dg 2 GF (t)/dt using a symmetric 4point numerical approximation for the derivative.

A. Infinite volume extrapolation
The finite volume effects depend on t/L 2 , and at leading order are proportional to t 2 /L 4 . Figure 3 shows g 2 GF as the function of a 2 /t for our set of bare coupling values on all three volumes for the Symanzik operator with Zeuthen flow (ZS). The L/a = 16 ensembles exhibit growing finite volume effects for a 2 /t 0.2, but the two larger volumes remain close up to a 2 /t 0.1. We monitor both g 2 GF and β(g 2 GF ) closely and restrict the GF time such that the finite volume corrections remain very small  5. Representative a 2 /t → 0 continuum limit extrapolation at g 2 GF = 4.8. We show results for two infinite volume extrapolations and three operators, fitting filled symbols in the range 2.00 ≤ t/a 2 ≤ 3.64 (0.500 ≥ a 2 /t ≥ 0.274). The (uncorrelated) fits are independent but predict consistent a 2 /t = 0 continuum values. and the leading order t 2 /L 4 term is sufficient to take the infinite volume limit.
The L/a → ∞ limit has to be taken at fixed t/a 2 flow time and coupling g 2 GF . We therefore determine g 2 GF (t) and its derivative on every ensemble and interpolate β(g 2 GF (t); L) with a 4th order polynomial for each lattice volume to predict the finite volume β function at arbitrary renormalized coupling. The top panels of Fig. 4 illustrate this at t/a 2 = 2.2 and 4.2 (a 2 /t = 0.455 and 0.238) for the ZS combination. We predict the infinite volume β(g 2 GF ) using a linear a 4 /L 4 extrapolation of the interpolated β(g 2 GF ) values. The lower panels of Fig. 4 show examples for g 2 GF values spanning the accessible range in our numerical simulation. Finite volume effects are negligible at small flow time and the extrapolations are mild, well described by a linear a 4 /L 4 dependence even at t/a 2 = 4.2. As a consistency check we compare extrapolations using all three volumes to extrapolations using the two largest volumes only. While the errors of the infinite volume predictions change, the values are consistent. Other flow and operator combinations show similar volume dependence.

B. Infinite flow time extrapolation
The final step is the a 2 /t → 0 continuum limit extrapolation at fixed g 2 GF . The GF time is a continuous variable but the range of t/a 2 values has to be chosen with care. First, the flow time must be large enough for the RG flow to reach the RT where irrelevant operators are suppressed. Assuming one irrelevant operator with scaling dimension α < 0 dominates the cut-off effects, an extrapolation in (t/a 2 ) α/2 predicts the continuum limit. Around the GFP α = −2 and we find that our data is well described by a linear a 2 /t dependence when t/a 2 2.0. Second, the upper end of the flow time range must also be restricted. When finite volume effects depending on t/L 2 are large, the L → ∞ extrapolations become unreliable. Figure 3 suggests that t/a 2 4.0 is sufficient to control this. Any change of the continuum limit prediction when varying the minimal and maximal flow time values can be incorporated as systematical uncertainty.
We show an example for the a 2 /t → 0 continuum extrapolation at g 2 GF = 4.8 in Fig. 5 where we fit the data (filled symbols) in the range 2.00 ≤ t/a 2 ≤ 3.64 (0.500 ≥ a 2 /t ≥ 0.274). In principle the flow time t is a continuous variable; in practice we evaluate g 2 GF (t) with finite step-size ε = 0.04 and choose to dilute the data by plotting and fitting only every third t value. Further we perform uncorrelated fits neglecting that values in t are correlated which could easily be accounted for in a bootstrap or jackknife analysis. Once sufficiently large flow times are reached, the lower flow times in the fit range impact the size of the uncertainties in the continuum limit and largest value of the coupling g 2 GF which can be reached on a given data set. Hence statistical precision and reach in g 2 GF has to be balanced. Higher flow times in the fit range are subject to finite volume effects and we truncate to avoid uncontrolled systematics.
We compare the β(g 2 GF ; t) functions obtained using Zeuthen flow and Wilson plaquette, Symanzik, and clover operators. (The clover operator does not reach g 2 GF = 5.2 and is therefore missing in the second panel.) We consider two different infinite volume extrapolations, using all three volumes or only the largest two. For illustration we show additional data at larger flow time using open symbols. The linear extrapolations in a 2 /t shown by the colored bands are obtained independently for each operator. Their excellent agreement at the a 2 /t = 0 limit is a strong consistency check of the GF time range and the infinite volume extrapolation. Further consistency checks are possible by considering different flows. The scaling exponent α of the leading irrelevant operator could also be extracted when performing simultaneous fits to several operators. DISCUSSION We summarize the final result of our calculation by showing in Fig. 6 the continuous GF β function predicted from nine different flow and operator combinations. For reference we add the universal 1-and 2-loop perturbative predictions. The different flow/operator combinations are barely distinguishable, only the Wilson flow with clover operator (WC) prediction is off by ∼1σ. All others are in nearly perfect agreement. The non-perturbative prediction follows the universal perturbative curves up to g 2 GF ≈ 2.0, but at stronger couplings is closer to the 1-loop prediction. A similar trend is observed in Ref. [21,47] for 3-flavor QCD, suggesting that the GF coupling runs slower than the MS coupling. Since the RG β function is scheme dependent, the GF and MS schemes do not have to agree. However the precise nonperturbative running is needed to determine α s , the Λ s parameter, or connecting lattice simulations to continuum results [47,48].
The continuous β function described here works in both conformal or infrared free systems. The relation between GF and Wilsonian RG is especially useful in strongly coupled systems. The method complements the step scaling function studies, providing a new handle on systematic errors. Our proposed analysis can be repeated with existing GF measurements generated for step-scaling function studies without additional costs [49].

ACKNOWLEDGMENTS
We are very grateful to Peter Boyle, Guido Cossu, Anontin Portelli, and Azusa Yamaguchi who develop the Grid software library providing the basis of this work and who assisted us in installing and running Grid on different architectures and computing centers. We are indebted to Daniel Nogradi for extending his original calculation on symmetric volumes to asymmetric lattices and sharing the result prior to publication. We thank Alberto Ramos for many enlightening discussions during the "37th International Symposium on Lattice Field Theory", Wuhan, China, and Rainer Sommer and Stefan Sint for helpful comments. We benefited from many discussions with Thomas DeGrand, Ethan Neil, David Schaich, and Benjamin Svetitsky. A.H. and O.W. acknowledge support by DOE grant DE-SC0010005. A.H. would like to acknowledge the Mainz Institute for Theoretical Physics (MITP) of the Cluster of Excellence PRISMA+ (Project ID 39083149) for enabling us to complete a portion of this work. O.W. thanks the CERN theory group for their hospitality during the final stages of preparing this manuscript and acknowledges partial support by the Munich Institute for Astro-and Particle Physics (MI-APP) of the DFG cluster of excellence "Origin and Structure of the Universe".
Computations for this work were carried out in part on facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy and the RMACC Summit supercomputer [50], which is supported by the National Science Foundation (awards