Λ parameter of the SU(3) Yang-Mills theory from the continuous 𝛽 function

Nonperturbative determinations of the renormalization group 𝛽 function are essential to connect lattice results to perturbative predictions of strongly coupled gauge theories and to determine the Λ parameter or the strong coupling constant. The continuous 𝛽 function is very well suited for this task because it is applicable both in the weakly coupled deconﬁned regime as well as the strongly coupled conﬁned regime. Here we report on our results for the 𝛽 function of the pure gauge SU(3) Yang-Mills theory in the gradient ﬂow scheme. Our calculations cover the renormalized coupling range 𝑔 2 GF ∼ 1.2 −27 , allowing for a direct determination of √ 8 𝑡 0 Λ MS in this system. Our prediction, √ 8 𝑡 0 Λ MS = 0 . 622(10) , is in good agreement with recent direct determinations of this quantity.


I. INTRODUCTION
The renormalization group (RG) function is defined as the logarithmic derivative of the renormalized running coupling 2 with respect to an energy scale and describes the scale dependence of the renormalized coupling of a 4-dimensional gauge-fermion system.Precise determination of the function is essential to understand the nonperturbative running of the gauge coupling and to predict, e.g., the Λ parameter in quantum chromodynamics (QCD) where 2 ( ) is the renormalized coupling at the energy scale and the constants 0 and 1 denote the universal 1-and 2-loop coefficients of the RG function.Beyond 2-loop the function is scheme dependent, but it is straightforward to connect the Λ parameter obtained using a different scheme to the phenomenologically favored MS scheme using the 1-loop perturbative relation between the corresponding running couplings.
In the chiral limit of QCD, the Λ parameter sets the physical scale.Its accurate determination allows the prediction of the strong coupling constant ( ) at ≥ , the mass of the boson [1,2].Λ parameters of systems with different flavor numbers can be connected by nonperturbative decoupling relations.In particular, it is possible to predict ( ) in = 3 QCD from the pure Yang-Mills Λ parameter [3].A precise calculation of the Λ parameter in the pure Yang-Mills system is therefore desirable.
The RG transformation of an asymptotically free system has an ultraviolet (UV) fixed point (FP) on the 2 0 = 0 critical surface, where 2 0 denotes the (marginally) relevant bare coupling.The UVFP and the renormalized trajectory (RT) emerging from the UVFP describe continuum physics.The renormalized gauge coupling in Eq. ( 1) "measures" the flow along the RT.The gradient flow (GF) transformation [4][5][6] is a particularly promising choice to define a renormalization scheme, determine the running coupling 2 GF ( ), and calculate a nonperturbative function using lattice simulations.While on its own the GF transformation does not describe an RG transformation because it does not reduce the degrees of freedom of the system, it can be modified to define a complete RG scheme if a coarse-graining step is incorporated when defining expectation values [7,8].In this setup, the energy scale is set asymptotically by the GF flow time For a generic RG flow, there is a one-to-one correspondence between running renormalized coupling and local observables that do not renormalize [9].In the GF scheme the renormalized coupling can be defined in terms of the continuum flowed energy density ( ) as since ( ) has no anomalous dimension [6].The normalization factor  = 128 2 ∕(3 2 − 3) in Eq. ( 4) is chosen such that the renormalized coupling 2 GF matches the renormalized coupling of the MS scheme at tree-level.
The GF coupling offers a convenient low-energy hadronic scale 1∕ √ 8 0 , with the flow time 0 defined where the GF coupling takes a specified value, 2 GF ( 0 ) ≡  × 0.3 [6].The 0 scale has been used extensively in simulations to set the lattice scale and its continuum value is known in Yang-Mills systems [6,10] and in QCD [1] with 2+1 flavors [11][12][13] and 2+1+1 flavors [14][15][16].Further, we can use 2  GF to define the RG function in the GF scheme GF ( 2 GF ).If GF is determined up to 2 GF ( 0 ) ≈ 15.8, we can obtain √ 8 0 Λ GF by integrating Eq. ( 2) up to 2  GF ( 0 ).The finite volume step-scaling function [10,17,18] is a commonly applied method to determine the function.How- ever, it requires that the lattice size is the only dimensionful quantity of the system, preventing its application in the confining, chirally broken regime of QCD.Since the 0 scale corresponds to a low energy hadronic scale in the confining regime, the step-scaling method cannot determine the function at strong enough gauge couplings to directly predict the Λ parameter at the 0 scale.An alternative approach is to determine the infinite volume continuous function as described in Refs.[19][20][21].While this method requires separate infinite volume and continuum limit extrapolations, it is applicable in the confining regime.
In this paper, we report on our findings on the continuous function up to and even beyond the renormalized gauge coupling2 GF ( 0 ) in the pure gauge SU(3) Yang-Mills theory.Using the continuous function method we are able to determine the scale dependence of the running coupling in the confining regime.We do not need any other hadronic observable to determine the Λ parameter at the 0 scale.We compare our value of the Λ parameter to other determinations from Refs.[10,22] and values reported in the FLAG 2021 review [1].Preliminary results of our work were reported in Ref. [23].A similar calculation of the continuous function of SU(3) Yang-Mills theory reports preliminary results in Ref. [22].
This paper is organized as follows.In Sec.II we summarize the details of our numerical simulations and explain our determination of the continuous function in Sec.III.In the subsequent section we discuss the determination of the Λ parameter and close by discussing our results in Sec.V.

II. NUMERICAL DETAILS
Our study is based on simulations performed using the tree-level improved Symanzik (Lüscher-Weisz) gauge action [24,25].We consider nineteen bare gauge couplings and generate configurations with periodic boundary conditions (BC) in all four directions using the hybrid Monte Carlo (HMC) update algorithm [26] as implemented in GRID [27,28].We set the trajectory length to = 2 molecular dynamics units (MDTU) and save configurations every 10 trajectories (20 MDTU).Subsequently, we use QLUA [29,30] to perform gradient flow measurements.
Table I lists the number of thermalized configurations analyzed for each bare coupling and volume as well as the HMC acceptance rates which all range between 70% and 90%.In the strong coupling regime1 ( ≤ 4.9) we have four volumes, (( ∕ )4 = 20 4 , 24 4 , 28 4 , 32 4 ), while at weaker couplings ( ≥ 5.0) we add a fifth, 48 4 volume. 2 The Polyakov loop expectation value in Fig. 1 shows that most ensembles for ≥ 5.50 are in the deconfined phase, while most ≤ 5.0 ensembles are confining. 3Since = 5.30 sits in the transition region, we discard it from our main analysis but use it to check for systematic effects later on.We study the autocorrelation of the renormalized coupling 2  GF using the Γ-method [31] and typically find integrated autocorrelation times of less than four measurements (80 MDTU).However, near the transition from the deconfined to the confined regime, the integrated autocorrelation times increase up to fifteen measurements (300 MDTU).
While simulations in the deconfined regime have zero topological charge, nonzero topological charges are expected in the confined region.As we decrease for a fixed volume, we indeed observe that nonzero topological charges arise and that their fluctuations increase.We study the effect of nonzero topological charge on 2  GF by "filtering" configurations according to topological sectors and compare the corresponding values of 2 GF .Our data suggest that for the fast-running pure gauge system, the function is sufficiently large that the impact of nonzero topological charges is not statistically resolved. 4

A. Definition of the GF coupling and function
We define the finite volume gradient flow coupling on the lattice [6,18] as  4 .Configurations are separated by 20 MDTUs.All ensembles with > 5.5 are in the deconfined phase, while most with ≤ 5.50 either transition between confined and deconfined or are confined.The values for = 5.30 are set in italics because these ensembles are not part of the main analysis and are only used to estimate systematic effects.
where 2 0 = 6∕ is the bare gauge coupling.The term ( , ) = 1 + ( , ) corrects for the zero modes caused by the periodic BC of the gauge fields [18].The energy density ⟨ ( )⟩ can be approximated by local gauge observables.In general, 2 GF ( ; , 2 0 ) also depends on the gauge action, the gradient flow, and the operator chosen to approximate ⟨ ( )⟩.In this work we use Symanzik improved gauge action and consider two different gradient flow transformations, Wilson (W) and Zeuthen (Z) flow.We determine the Wilson plaquette (W), the clover (C), and the tree-level improved Symanzik operator (S) to estimate ⟨ ( )⟩.It is possible to calculate the treelevel lattice corrections to 2 GF ( ; , 2 0 ) for a given action-flowoperator combination and include those in ( , ) [33].In our analysis, we consider ( , ) with and without tree-level corrections, and we refer to the former as "tree-level normalization" (tln).We use a shorthand notation to distinguish the different flow-operator combinations, e.g.ZS refers to Zeuthen flow and Symanzik operator.When using tln in the definition of 2 GF ( ; , 2 0 ) we prepend "n", i.e. nZS refers to tln improved Zeuthen flow and Symanzik operator.In the continuum limit, the RG function is independent of the bare action, GF transformation or operator choice, and the comparison of the different combinations can serve to estimate systematic effects.Continuum extrapolations on the lattice are more stable when lattice artifacts are small, and we chose the nZS combination as our preferred analysis.
From the GF coupling in Eq. ( 5), we derive the GF func- in finite volume by discretizing the flow-time derivative d∕d with a 5-point stencil.The interval of the time derivative is set by the time step used to integrate the gradient flow equation for the gauge field.We choose = 0.04 and explicitly verified that this choice has no impact on our analysis by repeating the GF measurements using = 0.01 for selected ensembles.
Our numerical analysis starts by determining the renormalized couplings 2 GF ( ; , 2 0 ) at all flow times for our three operators.We determine 2 GF ( ; , 2 0 ) for both data sets obtained with Zeuthen and Wilson flow, respectively, and for our entire set of gauge field ensembles.Next, we numerically calculate the derivative as specified in Eq. ( 6).The uncertainties are propagated using standard correlated error propagation techniques implemented in the software packages gvar [34] and lsqfit [35].
The analysis of the RG function proceeds by first taking the infinite volume limit ∕ → 0 of both 2 GF ( ; , 2 0 ) and GF ( ; , 2 0 ) independently at fixed and ∕ 2 .This is followed by an interpolation of GF ( ; 2 0 ) in 2 GF ( ; 2 0 ) at fixed ∕ 2 .The last step of our analysis is to take the 2 ∕ → 0 continuum limit at fixed 2 GF .These steps are detailed in the rest of this section.

B. Infinite volume extrapolation
In a 4-dimensional gauge-fermion system the volume is a relevant parameter, and the RG equation in finite volume includes a term describing its effect on the running of the renormalized coupling.We prefer to avoid complications arising from this a priori unknown quantity and define the RG function in the infinite volume limit.This requires first extrapolating to infinite volume before considering the continuum limit.
In Fig. 2 we show typical infinite volume extrapolations for 2 GF ( ; , 2 0 ) and GF ( ; , 2 0 ) at relatively weak coupling ( = 6.0), at intermediate couplings ( = 4.9 and 5.50), and in the strong coupling regime ( = 4.35).Each panel shows the ( ∕ ) 4 extrapolation at five different lattice flow time values that cover the range we use in the continuum limit extrapolation (cf.Sec.III E).In the strong coupling regime, we observe very mild volume dependence.This is consistent with the expectation that below the confinement scale conf , i.e. at renormalized couplings stronger than 2 GF ( conf ), the confinement scale provides an infrared cutoff and the volume dependence is suppressed.Thus we find it sufficient to use volumes with ∕ = {20, 24, 28, 32} in this regime.At renormalized couplings that correspond to energy scales above the confinement scale, the volume dependence is more significant, as the only infrared cutoff is due to the finite volume.Here we drop the ∕ = 20 ensembles and add 48 4 lattices in the infinite volume extrapolation.At = 5.5 and 6.0 volumes ∕ ≤ 32 appear deconfined, while ∕ = 48 is transitioning to deconfined at = 6.0 and is confined at = 5.5.Similarly the ∕ = 32 at = 5.5 ensemble is transitioning and exhibits very long autocorrelation times.Likely these long autocorrelations result in underestimated statistical errors causing the 2 deviation.At = 5.0 the Polyakov loop expectation value shows that volumes ∕ ≥ 28 are confining, while ∕ = 20 and 24 are in the transition region.∕ = 24 shows especially large autocorrelation times.Nevertheless, at = 5.0 we use all five volumes in the infinite volume extrapolation, though it is dominated by the larger volumes.We always perform the infinite volume extrapolation using a linear fit ansatz in ( ∕ ) 4 .For flow times 2.0 ≲ ∕ 2 ≲ 4.0 contributing to the continuum extrapolations in Sec.III E, most fits have good -values (≳ 10%).Notable exceptions are = 5.0 and 5.5, where large autocorrelations make it difficult to estimate the 1.00 1. 25  errors correctly.When adding = 5.3 to our analysis, the infinite volume extrapolations have vanishing -values.The situation did not improve despite extending the Monte Carlo simulations and using larger thermalization cuts.We therefore conclude that = 5.3 suffers from sitting in the transition region and discard it from our main analysis.

C. Different flows, operators and tree-level improvement
In Figs. 3 and 4 we show the infinite volume extrapolated GF coupling 2 GF ( ; 2 0 ) as a function of the gradient flow time ∕ 2 with and without tln improvement at = 6.0 and 4.35, respectively.We consider both Zeuthen and Wilson flow and all three operators.In the continuum limit the different flows and operators must agree, so any difference between them at finite bare coupling points to cutoff effects.Figure 3 showcases the improvement due to tln at weak coupling.As the left panels show, there are significant differences between the three operators even at flow time ∕ 2 = 4.0.Comparing the top and bottom panels on the left, large differences for 2 GF ( ; 2 0 ) between Zeuthen (top) and Wilson (bottom) flow can be seen.Both the operator and flow dependence are largely removed by tln normalization, which is demonstrated in the panels on the right.
Figure 4 shows the same quantities at a strong bare coupling, = 4.35.The gauge coupling runs fast and the same flow time range covers nearly 20 times the range of 2 GF at = 4.35 compared to = 6.0.While the different operators and flows without tln correction may appear to be closer at = 4.35 than at = 6.0, the absolute difference is much larger.Treelevel correction works similarly in the strong coupling, removing most flow and operator dependence.The somewhat larger spread after tln suggests that the perturbatively motivated improvement is not as effective in the strong coupling regime as at weak coupling.The effect of tln improvement at strong coupling was recently studied in the context of gradient flow scale setting for the case of 2+1 flavor gauge field configurations generated with Iwasaki gauge action and Shamir domain-wall fermions [38].That work also showed that tln removed most flow and operator dependence, but the improved predictions still had significant cutoff effects.However, here we study the pure-gauge system where no effect due to fermion loops can spoil the improvement and, furthermore, we use the perturbatively improved Symanzik gauge action in comparison to Iwasaki gauge action.In Sec.III E we will show that in our case tln improvement not only reduces the difference between different flows and operators but also reduces cutoff effects.
Based on the improvement from tln, we designate the treelevel improved Zeuthen flow with Symanzik operator combination (nZS) as our preferred analysis.In the remainder of this paper we therefore concentrate on analyzing this combi- nation.Other tln improved flow-operator combinations are consistent within uncertainties and considered as part of our discussion on systematic uncertainties in III G.Moreover, we repeat our analysis without tln improvement, which increases discretization errors.Within these more sizeable uncertainties, the unimproved combinations are consistent with our preferred nZS prediction.

GF
The last task before we can take the 2 ∕ → 0 continuum limit is the identification of { GF ; 2 0 ), ∕ 2 } pairs at fixed 2 GF ; 2 0 ).We achieve this by interpolating GF ( ; 2 0 ) in terms of 2 GF ( ; 2 0 ) at fixed values of lattice flow times ∕ 2 .Our interpolating form has to be able to describe the weak coupling regime where we expect perturbative behavior ( 2 GF ) ∼ − 0 4 GF +  6 GF as well as the strong coupling regime where it appears that GF ( 2 GF ) ≈ 0 + 1 2 GF .A polynomial is not appropriate to cover both regimes.Instead, we choose a form that is similar to the ratio of polynomials used in Padé approximations such that we enforce the leading order GF ( 2 GF ) ∼ 4 GF behavior.Since cutoff effects can change even the asymptotic behavior of the function, this ansatz could constrain the lower limit on the flow time ∕ 2 used in the analysis.However with = 4 we find that  2 GF provides a good description with -values between 17% -32% for flow times 2.0 ≲ ∕ 2 ≲ 4.0.The two panels of Fig. 5 illustrate the interpolation at various flow times in the range ∕ 2 ∈ [2.0, 4.0].The colored data points on both panels correspond to infinite volume limit ( 2 GF , GF ) pairs at each bare coupling and five selected flow time values.For reference we show the universal 1-and 2loop and the GF 3-loop [39] perturbative lines in gray using dashed, dotted, and dash-dotted lines, respectively.The colored bands describe the interpolation according to Eq. ( 7) with = 4.The top panel of Fig. 5 shows GF vs. 2 GF in the entire range of our data.The data points at different bare couplings and flow times form a smooth curve, indicating that cutoff effects are mild.The corresponding interpolating curves basically overlap, creating the purple-hued band.This plot also illustrates that the function is approximately linear in the strong coupling regime.We will discuss this feature further in Sec.III H.In the insert of the top panel we zoom in on the weak coupling regime.To enhance the small 2 GF behavior, we plot GF ( 2 GF )∕ 4 GF vs. 2 GF on the bottom panel.At large flow times and small couplings, GF ( 2 GF ) has to be consistent with the perturbative predictions.As the panel shows, both the data and the interpolating forms are close to the perturbative values in our selected flow time range, though we do not constraint the intercept of GF ( 2 GF )∕ 4 GF .

E. Continuum limit extrapolation
To obtain the continuous function in the continuum limit, we need to perform one last step, the 2 ∕ → 0 continuum extrapolation.We use a linear extrapolation in 2 ∕ → 0 and the flow time ∕ 2 has to be chosen large enough for the RG flow to reach the renormalized trajectory where a linear extrapolation describes the remnant cutoff effects (within statistical uncertainties).We also have to limit the maximum flow time because the infinite volume extrapolation is reliable only when the finite volume corrections follow the leading order 2 ∕ 4 scaling form.In practice, we vary min ∕ 2 and max ∕ 2 and monitor the quality of the linear extrapolation.In Fig. 6 we show examples of the continuum extrapolation from our weakest available coupling of 2 GF up to the 0 scale of  are highly correlated.We obtain the final continuum limit prediction by selecting a discrete set of ∕ 2 values.For this set we perform an uncorrelated fit and estimate its statistical uncertainty by repeating this fit using the central values shifted by ±1 .This avoids the complication of inverting a poorly conditioned correlation matrix and ensures we are not underestimating the statistical uncertainty.

F. The nonperturbative GF function
Figure 7 shows our result for the nonperturbative GF function for pure gauge SU(3) Yang-Mills in the gradient flow scheme with statistical errors only.In addition to the continuum limit prediction shown with a salmon-colored band, the plot also shows the infinite volume extrapolated nZS lattice data at flow times ∕ 2 ∈ [2.0, 4.0] where, for better visibility, we "thin" the data and only show every fifth data point, i.e. flow time values are separated by Δ ∕ 2 = 0.2.The nZS combination shows very little cutoff dependence and the raw lattice data sits on top of the continuum extrapolated value.Overall, our results span the coupling range from 2 GF ≳ 1.2 up to 2 GF ≲ 27, well into the confining regime of the system.

G. Systematic uncertainties of the GF function
In addition to the statistical uncertainties shown for the final result of our GF function in Fig. 7, we check for systematic effects by considering variations of our preferred nZS analysis.We discuss the different variations below and show the outcome as relative changes w.r.t.our nZS analysis in Fig. 9.
• For the interpolation in 2  GF we use a ratio of polynomials similar to Padé approximation.We alter the functional form by changing the order of the polynomials in Eq. ( 7) from = 4 to = 2 or 6, but find this is only a subleading effect on our resulting GF function.In Fig. 9 we show the (larger) effect by taking the absolute difference between interpolations based on = 4 and = 2.
• Our simulations at = 5.30 sit right at the deconfined/confined transition and we observe poor -values when performing the infinite volume extrapolation at this bare gauge coupling.We therefore discarded = 5.30 from our main analysis, but we use it now to check for systematic effects by including it as an additional data point in a sensitive region of the GF function.As expected, the impact of adding = 5.30 is largest around 2 GF ∼ 5. • To validate our infinite volume extrapolation, we consider two variations: 1. We drop the smallest volume and perform a linear fit to our three largest volumes.
2. We note that in Fig. 2 our largest volume ∕ = 32 ( ∕ = 48) for < 5.00 ( ≥ 5.00) is very close to the extrapolated infinite volume value.Therefore we simply repeat our analysis using only the largest volumes.
Both variations result in comparable uncertainties, but just using the largest volumes has a slightly larger effect.Hence we use the latter to estimate finite volume effects.
As can be seen in Fig. 9 this could be our dominant systematic uncertainty for strong couplings 2 GF ≳ 11.
• We test the continuum limit extrapolation by varying the range of the flow times entering the linear fit in 2 ∕ at fixed values of 2 GF .Keeping max ∕ 2 = 4.0 fixed, we vary min ∕ 2 from 1.52 to 2.0.Similarly, we vary max ∕ 2 from 4.0 to 5.0 while min ∕ 2 = 2.0 is kept fixed.Comparing these to our preferred analysis we see at most a variation of (0.3%) in the central value of our GF function.We show the maximum of these variations as "fit range" uncertainty in Fig. 9. Further, we repeat the fits using our preferred fit range but reduce the number of data points fitted by increasing the separation in flow time.As expected this variation results only in minuscule changes and can be neglected.

• In the continuum limit different flows and operators
should predict the same renormalized GF function.
Qualitatively this is the case, as Fig. 8 demonstrates.There we show six different tln improved flow-operator combinations to determine GF which all sit on top of each other forming nearly a single band.Looking at the relative changes in Fig. 9 we do, however, see deviations of (0.8%).Since for couplings in the range 7 ≲ 2 GF ≲ 11 that effect is larger than other systematic effects, we conservatively include these variations when obtaining our systematic uncertainty.
• We further check for consistency by analyzing our data without using tln improvement.We find that removing the tln improvement increases the discretization errors noticeably, but within the larger uncertainties, the results are consistent with our preferred analysis.
Using the information compiled in Fig. 9, we obtain the total uncertainty of our GF function by adding our statistical error

H. GF in the strong coupling regime
As can be seen in particular in Fig. 8, the GF function in the strong coupling regime ( 2 GF ≳ 13) is approximately linear.Studying the derivative d GF ∕d 2  GF , we observe a plateau within errors in the range 20 ≲ 2 GF ≲ 27.In that range we can approximate GF ( 2 GF ) ≈ 0 + 1 2 GF with 0 = 5.8(3), 1 = −1.32(1).This is very different from the perturbative prediction that would suggest a polynomial form with terms of 4  GF and higher order.Clearly, nonperturbative effects are at play here.We may use the above observation to predict the flow time dependence of the flowed energy density ⟨ 2 ( )⟩ ∼ 0 + 1 ∕ ̃ ) − 1 , where 0 = − 0 ∕ 1 ≈ 4.4 (2) and ̃ is an integration constant.We do not have a physical explanation for this behavior, though it is tempting to think that the form of ⟨ ( )⟩ as → ∞ is related to instantons, the only objects in the vacuum after UV fluctuations are removed by the flow.We mention that Refs.[41][42][43][44] hypothesize that in confining systems the RG function is linear in the confining, strongly coupled regime.References [41,42] attempt to solve the strong CP problem by connecting ⟨ ( )⟩ to the instanton density at large flow time.However, this work considers only one bare gauge coupling, does not take the continuum limit and uses much larger flow times than we do on similar volumes.While our results are consistent with a linear function in the strong coupling range, our slope 1 = −1.32(1) is significantly smaller than -1 obtained in Refs.[41,42].More- 5 Our final result for GF is also provided as an ASCII file [40].
GF in the weak coupling region.The salmoncolored band shows our nonperturbatively determined GF with combined statistical and systematic uncertainties.We match to the 3-loop GF function using Eq. ( 8) in the range 2  GF ∈ [1.4,1.8] indicated by the grey hatched area.Shifting our nonperturbative values by ±1 , we obtain the magenta bands providing the upper and lower limits of the resulting matched function shown in blue.
over, we observe a nonzero intercept.It would be interesting to understand the nonperturbative origin of the linearity of the GF function in the strong coupling confining region both in pure Yang-Mills systems and, if it persists, with dynamical fermions.

A. Matching the perturbative regime
In order to determine the Λ parameter according to Eq. ( 2), we need to know GF ( 2 GF ) down to 2 GF = 0.In the weak coupling regime we expect to recover the perturbatively predicted function.In Fig. 11 we show GF ∕ 4  GF to emphasize the weak coupling behavior of our nonperturbative result.The numerical data, shown by the salmon-colored band, is close to the 3-loop GF perturbative curve at our weakest gauge coupling 2 GF ≈ 1.2, but it does not yet connect smoothly.To remedy this limitation we extend our numerically determined GF function to the region below our weakest coupling data point ( 2 GF ≈ 1.2) using the parameterization where 0 , 1 , 2 are the 1-, 2-and 3-loop GF perturbative coefficients from Ref. [39] and is free.We determine by integrating the inverse function using both our numerically determined GF function and also using 4 parametrized by for ( ).We determine by equating these two values.The integration limits are set to cover the regime where we attempt to match GF and 4 .In the example shown in Fig. 11, we choose to match in the range 2 = 1.4, 2 = 1.8.In Fig. 12 we demonstrate that varying 2 and 2 has negligible impact on our value of √ 8 0 Λ MS .To account for the combined statistical and systematic uncertainty of our nonperturbative GF function, we repeat the matching shifting the central values by ±1 .That way we obtain the purple bands resulting in the upper and lower end of the blue band connecting our nonperturbative result to 2 GF = 0.

B. Calculating the Λ parameter
The final step of this analysis is to calculate √ 8 0 Λ GF by integrating Eq. (2) using the combination of our nonperturbative GF ( 2 GF ) for 2 GF ≥ 1.4 and the matched 4 ( 2 GF ) function for 2 GF < 1.4.The upper integration limit is set according to 2 GF ( 0 ) ≈ 15.8.Our prediction is where the error accounts for the uncertainty of our GF function as well as the uncertainty encountered due to the matching procedure.As a last step, we convert our prediction to the MS scheme using the exact 1-loop relation and find √ 8 0 Λ MS = 0.622 (10).

C. Comparison of √ 8 0 Λ MS determinations
The Yang-Mills Λ parameter has been studied previously using different approaches.Only the recent gradient flow studies in Refs.[10,22] directly predict the combination √ 8 0 Λ MS .In addition, we compare our value to determinations of 0 Λ MS listed by the flavor lattice averaging group (FLAG) [1] to meet the quality criteria to enter the average.These determinations are obtained using Schrödinger functional step-scaling methods [45,49], Wilson loops [46,48], or the short distance potential [47].We use the values quoted by FLAG 2021 for 0 Λ MS and convert them to √ 8 0 Λ MS using √ 8 0 ∕ 0 = 0.948(7) [6] (open symbols) or √ 8 0 ∕ 0 = 0.9414(90) [10] (filled symbols).Following the FLAG convention, we refer to the different results in Fig. 13 using either the name of the first author or, if applicable, the name of the collaboration and the two-digit year.
Given the spread in the values of √ 8 0 Λ MS , further scrutiny and understanding are needed before obtaining an average.We note, however, that the three most recent predictions are all mutually consistent.The high-precision result of Ref. [10] was re-affirmed in Ref. [50] using an alternative approach with better control over the continuum extrapolation.The estimate given in Ref. [42] is also consistent with these predictions.A possible source of difference to the older determinations is the conversion of 0 to √ 8 0 .

D. Nonperturbative matching of different schemes
In our analysis, a considerable systematic uncertainty arises from the weak coupling limit 2 GF ≲ 1.4.Our gradient flow setup is not efficient at weak coupling.It would be more economical to use data from existing calculations, e.g. the high precision Schrödinger functional data of Ref. [10] in the 0 < 2 GF < 1.4 regime and match it nonperturbatively to our data.
Such a matching requires finding the relation between our turns Eq. ( 12) into a straightforward fitting problem with undetermined coefficients.The only constraint is to identify and use the renormalized coupling range in the fit where the two schemes overlap.Such a nonperturbative matching and combination of different schemes could lead to a significantly improved prediction.Although we do not explore this method in the present analysis, it is worth considering in the future.

V. DISCUSSION
In this paper, we present a nonperturbative determination of the renormalization group function for the pure gauge Yang-Mills action.Using the gradient flow based continuous RG function, we present results for a wide range of values of the renormalized running coupling.Our results span the range of the perturbative weak coupling region 2 GF ≈ 1.2 up to the strongly coupled regime at 2 GF ≈ 27.This showcases the advantage of the continuous RG function because the continuous infinite volume function can be extended without limitation to the confining region.We also demonstrate the effectiveness of tree-level improvement of the gradient flow even in the strong coupling regime.
We investigate various sources of systematical uncertainties.For most of the 2 GF range covered, our statistical uncertainties are around 0.6%.In the weak coupling region, statistical and systematic errors are of similar size.At intermediate to strong coupling, we observe an increase in the systematic errors to 1.5% due to enhanced finite volume effects.
While in the weak coupling our results are close to the perturbative values, we observe in the confining regime that the GF function depends approximately linearly on the running coupling, implying a scaling relation of the flowed energy density ⟨ 2 ( )⟩ ∼ 0 + 1 ∕ ̃ − 1 with exponent 1 ≈ −1.32 (1).
This observation could be related to the topological structure of the vacuum, a possibility that warrants further investigation.
In the weak coupling regime, we are able to match our numerical results to the 3-loop GF function by extending the perturbative expression with a single 10 GF term.This matching allows us to predict the Λ parameter in the GF scheme.Using the perturbatively determined relation of the GF coupling 2 GF and the MS coupling, we obtain √ 8 0 Λ MS = 0.622 (10), where the error combines statistical and systematic uncertainties.This value is in good agreement with recent direct determinations of √ 8 0 Λ MS [10,22].A significant source of the systematic uncertainties in determining √ 8 0 Λ MS arises from the weak coupling regime.We outline a nonperturbative matching procedure to combine existing high-precision data in the weak coupling and our nonperturbative GF function that extends into the confining regime even beyond 2 GF ( 0 ).Such a combined determination could lead to a sizable reduction of the uncertainties on √ 8 0 Λ MS .

= 20 L/a = 24 L/a = 28 L/a = 32 L/a = 48 FIG. 1 .
FIG.1.The expectation value of the Polyakov loop on our ensembles.We measure the Polyakov loop at flow time ∕ 2 = ( ∕ ) 2 ∕32, and average the absolute values over all four directions.

FIG. 3 .
FIG. 3. Infinite volume extrapolated gradient flow coupling without tln (left) and with tln (right) as a function of ∕ 2 at = 6.0.Different colors correspond to different operators.The energy density is obtained for Zeuthen flow in the top and Wilson flow in the bottom panels.

FIG. 5 .
FIG. 5. Interpolation of GF in terms of 2GF at fixed flow times.Different colored symbols correspond to (2GF , GF ) pairs at flow times ∕ 2 = 2.0 -4.0.The colored bands are the interpolating functions based on ratios of polynomials given in Eq.(7).In the lower panel, we plot GF ( 2 GF )∕ 4 GF to enhance the weak coupling regime.

FIG. 6 .
FIG.6.Continuum 2 ∕ → 0 extrapolation between 2 GF = 1.2 to 15.8.The results are shown in two separate panels to accommodate the increasingly faster running of the coupling.In all cases, we show all three operators with Zeuthen flow after tln improvement, though the different operators overlap and are barely distinguishable in the plot.The open symbols are not included in the extrapolation fit.They are shown to illustrate the linear behavior of the data even outside the region used in the fit.

2 FIG. 7 .
FIG.7.The predicted GF function (salmon colored band) overlayed with the infinite volume extrapolated data at different bare coupling (colored data points) for our main analysis based on nZS for flow times ∕ 2 ∈ [2.0, 4.0], separated by Δ ∕ 2 = 0.2.The insert magnifies the weak coupling region.The nZS combination shows very little cutoff dependence and the raw lattice data sit on top of the continuum extrapolated value.

FIG. 8 .
FIG. 8. Comparison of different continuum limit results obtained for tln improved flow operator combinations.All results overlap and form a single band.

β
FIG.10.Final result for GF as a function of2  GF for the coupling range relevant to determine the Λ parameter.The yellow inner band shows only the statistical uncertainty, whereas the red outer band shows the combined statistical and systematic uncertainties.

2 GF coupling and the coupling 2  2 
of another scheme  expressed as 2 GF = .The relation of the corresponding functions can be obtained using the chain rule applied to the derivative of 2 GF with respect to 2 , which leads to the simple relation GF Each color corresponds to a fixed lattice flow time between ∕ 2 = 2.0 (yellow) to ∕ 2 = 4.0 (red).Black symbols indicate data points that are included in the infinite volume extrapolation, whereas symbols shown in gray are not included and just shown for reference.
but show additional data points both at smaller and larger flow times to illustrate the linear behavior in ∕ 2 .It is worth pointing out that the flow time is technically a continuous variable and the corresponding GF ( 2 GF ) values 2 GF ( 0 ) = 0.3  ≈ 15.8.In all cases, we use ∕ 2 ∈ [2.0, 4.0] Systematic uncertainties with respect to our main analysis based on nZS.By varying different parts of our analysis one after the other, we calculate the relative changes of the central value and compare the size of the different systematic uncertainties (colored lines) to our statistical uncertainty (salmon-colored band).
FIG.12.Systematic uncertainty in our perturbative matching procedure due to choosing 2 or 2 .The value for our preferred choices for 2 and 2 is highlighted in red.