Static force from generalized Wilson loops on the lattice using the gradient flow

The static QCD force from the lattice can be used to extract $\Lambda_{\overline{\textrm{MS}}}$, which determines the running of the strong coupling. Usually, this is done with a numerical derivative of the static potential. However, this introduces additional systematic uncertainties; thus, we use another observable to measure the static force directly. This observable consists of a Wilson loop with a chromoelectric field insertion. We work in the pure SU(3) gauge theory. We use gradient flow to improve the signal-to-noise ratio and to address the field insertion. We extract $\Lambda_{\overline{\textrm{MS}}}^{n_f=0}$ from the data by exploring different methods to perform the zero-flow-time limit. We obtain the value $\sqrt{8t_0} \Lambda_{\overline{\textrm{MS}}}^{n_f=0} =0.629^{+22}_{-26}$, where $t_0$ is a flow-time reference scale. We also obtain precise determinations of several scales: $r_0/r_1$, $\sqrt{8 t_0}/r_0$, $\sqrt{8 t_0}/r_1$ and we compare these to the literature. The gradient flow appears to be a promising method for calculations of Wilson loops with chromoelectric and chromomagnetic insertions in quenched and unquenched configurations


I. INTRODUCTION
The Standard Model of particle physics is one of the most precisely tested theories.A precise knowledge of the Standard Model parameters is a necessary condition to work out accurate perturbative predictions, and to compare them with high-precision experimental measurements.Quantum chromodynamics (QCD) is the sector of the Standard Model that describes the strong interaction.It is a field theory based on the gauge group SU(3) that depends on just one coupling, g, or equivalently α s = g 2 /(4π).The coupling may be traded at any time with an intrinsic scale; in the MS scheme, this is Λ MS .Once renormalized, α s is small at energy scales much larger than Λ MS , a property known as asymptotic freedom, but it becomes of order 1 at energy scales close to Λ MS .At high energies, we can rely on weak coupling perturbation theory to compute QCD observables.
The value of Λ MS , or equivalently α s , at a large energy scale can be determined by comparing some high-energy observable computed in weak coupling perturbation theory with data.A viable alternative is to replace data with lattice QCD computations i.e., the exact evaluation of the observable in QCD via Monte Carlo computations.For the current status of Λ MS extractions from lattice QCD, see for example the recent reviews [1,2].While in the last several years lattice extractions of Λ MS have been mostly done in QCD with dynamical quarks, the interest towards the running of the coupling in the pure gauge version of QCD i.e., without dynamical quarks, also called quenched QCD has been recently reignited [3].With modern lattice methods, the extraction of the coupling from the pure gauge theory can be done much more precisely nowadays than in the past, when quenched calculations were the only viable option.
The static energy is a well-understood quantity in lattice QCD [4,5] that can be used to set the lattice scale [6].Furthermore, the static energy is an observable that can be used to extract Λ MS by comparing its perturbative expression with lattice data at short distances.The QCD scale Λ MS has been extracted from the static energy both in pure gauge [7][8][9][10][11] and with dynamical fermions [12][13][14][15][16][17][18].
In lattice QCD, the static energy suffers from a linear divergence and needs to be renormalized.In dimensionally regularized perturbation theory, this divergence becomes a renormalon of mass dimension 1.For these reasons, it presents some advantages to look at the derivative of the static energy, which is the static force.The static force is not affected by linear divergences in lattice QCD or by renormalons of mass dimension 1 in dimensional regularization, yet it contains all the relevant information on the running of the strong coupling, as this is entirely encoded in the slope of the static energy.The force is known up to next-to-next-to-next-to-leading logarithmic accuracy in the coupling [19][20][21][22][23][24].
The derivative of the static energy performed numerically on lattice data introduces additional systematic uncertainties.Therefore, it may be advantageous to compute the force directly from a suitable observable [25][26][27].This observable consists of a Wilson loop with a chromoelectric field insertion.A difficulty related to this observable is that field insertions on Wilson loops when evaluated on the lattice have a bad signal-to-noise ratio and a slow convergence to the continuum limit originating from the discretization of the field components.This was studied in a previous work [28].In this work, to overcome the difficulty, we rely on the gradient flow method [29][30][31] to improve the signal-to-noise ratio and to remove the discretization effects of the field components.Wilson loops smeared with gradient flow have been studied before in terms of Creutz ratios [32][33][34].In Ref. [35], the Wilson line correlator in the Coulomb gauge was measured at finite T with gradient flow.To our knowledge, this was the first time the gradient flow was applied to the force directly.Furthermore, this study also serves as a preparation for the study of similar Wilson loops with field insertions appearing in the computation of several observables in the context of nonrelativistic effective field theories.New methods for integrating the gradient flow equations, based on Runge-Kutta methods, have been developed and implemented during the last few years [36].The force in gradient flow, besides that with lattice QCD, can also be computed analytically in perturbation theory.In perturbation theory, the static force at finite flow time is known at one-loop order [37].
The paper is structured as follows: In Sec.II, we discuss the theoretical background: we introduce the gradient flow and the static force at zero and finite flow time, in continuum and on the lattice.The lattice setup is described in Sec.III, and in Sec.IV we show our numerical results and perform the continuum limits.Finally, in Sec.V, we discuss our results and extract Λ n f =0 MS .Preliminary results based on these data have appeared before in conference proceedings [38,39].

II. THEORETICAL BACKGROUND A. The static force
The static energy V (r) in Euclidean QCD is related to a rectangular Wilson loop W r×T with temporal extent from 0 to T and spatial extent r [40] by where a is the lattice spacing, g the strong coupling, Tr the trace over the color matrices, and P is the path-ordering operator for the color matrices.In dimensional regularization, the static energy has a renormalon ambiguity of order Λ QCD , while on the lattice there is a linear divergence of order 1/a coming from the self-energy of the Wilson line.
Both the perturbative and lattice problems manifest as a constant shift to the potential.This may be renormalized by fixing the potential to a given value at a given point r * : Alternatively, taking the derivative removes the divergent constant, and we obtain a renormalized quantity, the static force.
The static force F (r) is defined as the derivative of the static energy: In perturbation theory, the static force is known up to next-to-next-to-next-to-leading logarithmic order (N 3 LL) [19][20][21][22][23][24].On the lattice, this derivative is evaluated from the static energy data either with interpolations or with finite differences, which leads to increased systematic errors.It is possible, however, to carry out the derivative of the Wilson lines at the level of Eq. ( 2), and rewrite the force as [25][26][27]] where the expression in the numerator consists of a static Wilson loop with a chromoelectric field insertion on the temporal Wilson line at position t * , and r is the spatial direction of the quark-antiquark pair separation; t * can be chosen arbitrarily.Since both expressions for the force represent the same renormalized quantity, it holds in the continuum that F E = F ∂V = F .

B. Gradient flow
We rely on the gradient flow method [29][30][31] for measuring Wilson loops with and without chromoelectric field insertions.The gradient flow is a continuous transformation of the gauge fields toward the minimum of the Yang-Mills gauge action along a fictitious flow time τ F : where B µ (τ F ) are the flowed gauge fields at flow time τ F with the SU(3) QCD gauge fields A µ as the initial condition at zero flow time, S YM is the Yang-Mills action evaluated with the flowed gauge fields; G µν is the field strength tensor evaluated with the flowed gauge fields; and D µ is the gauge covariant derivative.The flow depends on the local neighboring gauge field values through the derivative of the action with respect to the gauge field at position x, and its characteristic range is given by the flow radius √ 8τ F .This results in a smearing that cools off systematically the ultraviolet physics and automatically renormalizes gauge-invariant observables [41,42].Furthermore, we introduce the reference scale t 0 [31], defined implicitly through the expectation value of the action density as The gradient flow equation is adapted for flowed link variables V τ F (µ, x) on the lattice as where S Gauge (V τ F ) is some lattice gauge action evaluated with the flowed link variables, ∂ x,µ is the derivative with respect to V τ F (x, µ), and U µ (x) is the original SU(3) link variable.The flowed link variables of a gauge field configuration depend uniquely on the initial gauge field configuration U -i.e., V τ F = V τ F (U ) -and flowed observables are obtained by replacing the original link variables with the flowed link variables: O(τ F ) = O| U =Vτ F .The flowed expectation value of O is evaluated on the flowed gauge ensemble and can be written as which is still a path integral with the Euclidean action S E of the original zero-flow-time theory.Therefore, on the lattice, the expectation value is given by where N is the number of gauge fields.We solve the gradient flow on the lattice by an iterative Runge-Kutta implementation for the SU(3) matrices.We use either a fixed step-size algorithm [31] or an adaptive step-size algorithm [36,43].

C. The perturbative static force at finite flow time
The one-loop formula for the static force in gradient flow is [37] with , n f the number of flavors, and a 1 = 31C A /9 − 10n f /9.The functions F 0 (r, τ F ) and F L NLO (r, τ F , µ) are given analytically with where γ E is the Euler-Mascheroni constant and M (a, b, z) is the confluent hypergeometric function defined by with (x) k = Γ(x + k)/Γ(x), and For r > √ τ F , we approximate F F NLO with the polynomial where C a = 109.358,C b = 43.8438,C c = 404.790,and c n values are listed in Table I.In practice, the computation of M (1,0,0) takes a considerable amount of time, which is a problem for fitting this function.Those terms depend only on the flow-time ratio τ F /r 2 ; hence, we precompute them on a fine flow-time ratio grid, and we use spline interpolations for further calls of the perturbative formula.The one-loop formula has an explicit dependence on the renormalization scale µ in the form of log(µ 2 r 2 ) and an implicit dependence through the perturbative strong coupling α s (µ).At zero flow time, the only scale is the distance r; therefore, setting µ = 1/r is a natural choice.At large flow time, r is negligible with respect to √ 8τ F , and therefore, the natural choice is µ = 1/ √ 8τ F .A parametrization of µ that interpolates between these two limiting cases is µ = 1/ √ r 2 + 8τ F [37].Because in our lattice calculations we are not collecting data at large flow time, we adopt in this work the more general parametrization   (20).from three loops, an ultrasoft (us) scale of order α s /r also enters the static force equations at zero flow time [19].For the scope of this paper, we set the ultrasoft scale to be µ us = C A α s (1/r)/(2r).
We renormalize the coupling in the MS scheme; hence, both α s and the scale Λ MS are defined in that scheme.The scale Λ MS can be obtained by comparing the perturbative expression of the force with lattice data.Since we work in the pure SU(3) theory, the comparison provides Λ 0 ≡ Λ n f =0 MS .The small-flow-time expansion of r 2 F reads with c L = −22/3.We remark that [−12β 0 − 6C A c L ] = 8n f , which is 0 in this study (n f = 0).This means that at small flow time, the static force approaches a constant behavior, and that corrections to it are smaller than α 2 s τ F /r 2 ; r 2 F (r, τ F = 0) is r 2 times the one-loop perturbative force at zero flow time.
The static force at zero flow time is known up to N 3 LL accuracy [19][20][21][22][23][24], and the higher loop contributions are crucial for the extraction of the Λ 0 parameter.To benefit from this knowledge, we model the flowed force with the 1-loop expression at finite flow time, and we demand it to converge to the expression at arbitrary order at zero flow time.In terms of equations, our model function is given by where r 2 F (r, τ F = 0) is the static force at a given order at zero flow time, and r 2 F 1-loop is the full one-loop expression in Eq. (15).In this way, we correct for the change of the force due to the flow time up to one-loop order.The accuracy of the flow-time correction is consistent with the three-loop accuracy of the zero-flow-time part, as long as the flow-time correction subleading to α 2 s τ F /r 2 is small compared to α 4 s .This appears to be the case in our study, where we invoke the restriction τ F /r 2 ≲ 0.05.
For the rest of the paper, we refer to Eq. ( 23) when dealing with the flowed static force at higher orders.We label the one-loop-order force [next-to leading order (NLO)] as F1l, the two-loop force [next-to-next-to leading order (N 2 LO)] as F2l, the two-loop force with leading ultrasoft logarithms resummed [next-to-next-to leading logarithmic order (N 2 LL)] as F2lLus, the three-loop force [next-to-next-to-next-to leading order (N 3 LO)] with F3l, and the threeloop force with leading ultrasoft logarithms resummed as F3lLus.For the reasons discussed in [15], we also restrict the present study to the F3lLus force, although the force at three-loop order with next-to-leading ultrasoft logarithms resummed [next-to-next-to-next-to leading logarithmic order (N 3 LL)] would be available.

D. The force on the lattice
The Wilson loops W r×T are constructed as the closed, path-ordered product of link variables, consisting of two straight spatial Wilson lines in the spatial plane separated by T in the temporal direction.The spatial Wilson lines have the length r in the direction r.The ends of both spatial Wilson lines are connected by two straight temporal Wilson lines.The static force can be obtained as the numerical derivative of Eq. (1) from the symmetric finite difference NS NT β a [fm] t0/a 2 N conf Label 20 3   Other methods of defining the derivative of Eq. (1) consist, for example, in using the derivative of interpolating functions; these methods, however, add additional systematic uncertainties.
The main purpose of this work is to obtain the static force directly by computing P W r×T gE j (r = r ĵ, t * ), which consists of inserting a discretized jth chromoelectric field component into the path-ordered product at the temporal position t * in one of the temporal Wilson lines.In general, t * is arbitrary.Nevertheless, we choose t * = T /2 for even-spaced separations, and an average over t * = T /2±a/2 for odd-spaced separations.This reduces the interactions between the chromoelectric field and the corners of the Wilson loop.We use the clover discretization for the fieldstrength tensor, where U µν is a plaquette in the µ-ν plane.This symmetric definition of the chromoelectric field corresponds to the symmetric center difference according to Eq. ( 25) at tree level.Finally, we replace a 2 F µν with a 2 F µν − 1 Tr(a 2 F µν )/3, which makes the components of the field-strength tensor traceless and corresponds to an a 2 improvement [44].The chromoelectric field components are accessible through the components a The direct determination of the force F E on the lattice follows from Eq. ( 6) and the discretized version of P W r×T gE j (r = r ĵ, t * ).The finite extent of the chromoelectric field through its discretization introduces additional self energy contributions with a nontrivial lattice spacing dependence.These self-energy contributions slow down the convergence to the continuum limit; they have been studied in lattice perturbation theory [45].They are absent in the force obtained through the derivative, F ∂V .Since both calculations provide the same physical quantity, we may set where the constant Z E reabsorbs the additional self-energy contributions at finite lattice spacing.If Z E = 1, no self-energy contributions are present, and we can assume that the quantity behaves in a trivial way in the continuum limit.Z E from the static force was investigated nonperturbatively in a former study [28], and it was found that Z E has only a weak r dependence.For a different discretization of the chromoelectric field insertion needed for determining transport coefficients, the renormalization constant Z E was computed up to NLO in lattice perturbation theory [46].More recent studies [47,48] rely on the gradient flow method to renormalize the field insertions.In this study, we also use gradient flow to show the renormalization property explicitly, and to improve the signal-to-noise ratio.In the rest of this work, the default force measurement is given in terms of the chromoelectric field -i.e., F ≡ F E -while we still call the force obtained from the derivative of the static energy ∂ r V ≡ F ∂V .

III. LATTICE SETUP AND TECHNICAL DETAILS
On the lattice, we measure a 2 F (t, r, T ), which, multiplied by r 2 /a 2 , yields the dimensionless quantity r 2 F (t, r, T ).We have data on a periodic N 3 s × N t (N s : spatial lattice extent, N t : temporal lattice extent) grid; in our case, the available grids are 20 3 × 40, 26 3 × 52, 30 3 × 60, and 40 3 × 80. Table II shows our lattice parameters.We use the scaling from [4], which is based on the scale r 0 , to fine-tune the simulation parameters.With this scaling, our lattices share a physical size of approximately (1.2 fm) 3 × 2.4 fm.We produce the lattice configurations using over-relaxation and heat-bath algorithms with Wilson action and periodic boundary conditions.
We solve the gradient flow equation ( 12) with a fixed step-size integrator for the coarsest lattice (20 3 × 40) and an adaptive step-size integrator for the finer lattices.In all gradient flow integrations, we use the Symanzik action.We measure pure Wilson loops and loops with chromoelectric field insertions.That way, we can determine both the static potential and its numerical derivative for obtaining the static force, and the force directly.In Appendix B, we show briefly the impact of gradient flow on the bare Wilson loops with and without chromoelectric field insertions.To find the reference scale t 0 , we use the clover discretization in Eq. ( 26) to measure ⟨E⟩ and solve Eq. ( 11) for t 0 .We use this reference scale to express our quantities r and τ F in units of √ t 0 and t 0 , respectively, and to perform the continuum limit as a 2 /t 0 → 0.
The adaptive step-size integrator changes the gradient flow step sizes after every integration step, dependent on the lattice configuration.This means that the lattice measurement is done at different flow-time grids for each individual lattice configuration.Therefore, we need to interpolate the data to a common flow-time grid among the different lattice configurations.We use simple spline interpolations, since the gradient flow is a continuous transformation of the fields that produces a continuous function of the flow time.The interpolation can be done to a fixed flow-time grid in physical units, or to a fixed flow-time-ratio τ F /r 2 grid.Data along a flow-time grid in physical units at a given fixed r can be easily presented on a flow-time-ratio axis by setting the x axis to τ F /r 2 .
We check the fluctuation of the topological charge and observe no full freezing of the topology.At our largest lattice (L40), the fluctuation of topology slows down, which increases the autocorrelation times of the topological charge.We inspect the autocorrelation of the L40 lattice more closely in Appendix C. The static energy (and consequently, its derivative) is known to be less affected by topological slowing down [49] than a bare gradient-flowcoupling measurement would be [50].However, to be safe, we block our data such that the block size is larger than the topological autocorrelation times we see on any of the ensembles and settle for 30 jackknife blocks per ensemble.While the static force measurement is only weakly correlated with the topological charge, the scale t 0 can have a stronger dependence on it.To counter this, we measure t 0 only on a subset of configurations, having longer Monte Carlo time in between the measurements.
We have performed the simulations at a constant physical box size and have not tested for finite-volume effects from varying the physical volume.In [51], the finite-volume effects were studied and found to be minimal for hybrid static energies for almost the same set of lattice parameters.We expect the finite-volume effect to be equally small for the static force.

IV. PREPARATORY ANALYSES
In this section, we analyze and prepare the raw lattice data, and we perform the continuum limits needed for the Λ 0 extraction.This preparation contains the plateau extraction for the static energy and the T → ∞ limit for the direct force measurement, which is covered in Sec.IV A. Renormalization properties of gradient flow are discussed in Sec.IV B, followed by the continuum extraction, worked out in Sec.IV C. To finalize this section, we investigate the behavior of the various scales r 0 , r 1 , and t 0 .

A. Plateau extraction
We extract the plateaus in the T → ∞ limit based on a procedure from [52] that relies on an Akaike information criterion.We summarize it here briefly.This procedure is applied to data at every fixed r and τ F combination for each lattice ensemble.
In the first step, we perform constant fits for all possible continuous ranges within T /a = 1 and a certain T max with a minimum of at least three support points, each fit minimizing χ 2 as where i 1 and i 2 are the indices defining the specific lower and upper limits of the range; f represents the model function, which, in our case, is the constant function f (x, a) = a 0 ; and a is the parameter vector, which consists of only one component for the constant fit.The dataset D is specified by y i , the force measurement at T = x i , and C ij is the covariance matrix along the T axis.
In the next step, we define, for a given fit and for a specific range i 1 to i 2 , the Akaike information criterion (AIC) as where k is the number of parameters of the fit function, k = 1 for the constant fit, and i 1 − i 2 is the relative number of discarded data points within the total range.In theory, we are required to have the absolute number of discarded data points, which is given by N tot − (i 1 − i 2 ); however, this corresponds to a global shift for all AIC values, which can be eliminated in the next step.
In the third step, we find the model probability of a specific fit range as where Z is a normalization constant such that the sum of the probabilities over all considered fit ranges reduces to 1: At this stage, we observe that a global shift of all AIC values has no effect on the model probability, since it is absorbed into a redefinition of Z. Thus, we perform the global shift AIC i1,i2 , because it makes the computations of the model probability on the most important ranges numerically more stable.
In the last step, we compute model expectation values and deviations with the given model probabilities.The final results and their deviations from the plateau fit are thus given by This procedure can be generalized to more complex model functions, which we use in this work for the continuum limit and the Λ 0 extraction.Furthermore, to achieve a better understanding of the process, we compute the average fit ranges and their deviations, defined equivalently to Eq. (34).
A crucial part of this procedure is the selection of a certain T max .In principle, T max can be chosen to cover the whole temporal range of the Wilson loops, since the information criterion should eventually select the important ranges.However, especially at larger T , the statistical errors for the covariance matrix are underestimated.These give inaccurate model probabilities for insignificant fit ranges.Therefore, we have to select a proper T max value to prevent this behavior.To find a suitable T max , we perform the whole procedure for several T max values and select the one where, for a small variation of T max , the final result stays invariant.This is justified by the assumption that the important fit ranges, selected by the information criterion, should be fully captured by T max .Hence, a small variation of T max should not modify the important fit ranges.The T max selection can be automatized by taking data up to a T max where the relative statistical error is less than 10 %, and then decreasing T max iteratively until an invariant range of T max is found.However, the automatized procedure does not work in every case, and sometimes a selection has to be made manually.To do so, it is enough to identify a representative selection of a few flow times and distances r, and use piecewise linear interpolations to cover the whole dataset.
The Akaike information procedure also provides us an error estimate of the plateau extraction.Nevertheless, to propagate the statistical error, we use jackknife resampling with 30 blocks and perform the plateau extraction for every jackknife block.If the resulting jackknife error is comparable with the fit error, we use only the jackknife samples to propagate the error.In several cases, however, the choice of the fit range is a significant systematic error source, and we need to take it into account.In these cases, we label the systematic error by the Akaike information criterion of the corresponding observable as σ2 AIC , and the statistical error as σ 2 stat .

B. Implications of the gradient flow
Gradient flow has an impact not only on the signal-to-noise ratio improvement, but also on reducing discretization effects that occur through self-interaction contributions.Therefore, we are also interested in the renormalization factor Z E of the chromoelectric field insertion from Eq. ( 28).We find Z E nonperturbatively by solving Eq. (28) for Z E , FIG. 1: Z E for all lattice sizes with increasing flow time.We see that for flow radii larger than one lattice spacing, the factor Z E becomes practically 1.
which gives the ratio of the numerical derivative of the static potential and the direct force measurement according to Eq. ( 6): In [28], it was shown that Z E (r) has a weak r dependence.Hence, we extract Z E as a plateau fit, similar to the T → ∞ limit discussed in Sec.IV A, over r at fixed τ F ; we keep r approximately between 0.3 r 0 and 0.65 r 0 .Figure 1 shows the result for Z E for all lattice sizes against the flow time.At minimal flow times, we obtain that Z E > 1, meaning that at minimal flow times, the direct force measurement is affected by self-energy contributions originating from the chromoelectric field discretization.Furthermore, we obtain that Z E = 1 within 1 % deviation for flow radii larger than one lattice spacing, which is required for reliable continuum limits.We recognize a small bump within the 1 % range for flow radii within 1.5 < √ 8τ F /a < 5 for all lattice sizes.This is expected, as the ∂ r V (r) part is a simple finite difference and hence only approximates the static force.We observe that this systematic difference between the definitions of the force seems to vanish at larger flow radii ( √ 8τ F /a > 6).In conclusion, we find that a minimum amount of flow time ( √ 8τ F > a) has to be applied to be in the regime where the gradient flow has practically eliminated the nontrivial discretization effects.

Continuum interpolation
Preparing for the continuum limit, we interpolate r 2 F on all ensembles at a fixed flow time to a common r range in √ t 0 units.We use polynomial interpolations at different orders: with a k,n being the coefficients.We have different coefficients for different fixed-order polynomials.An interpolation at fixed order corresponds to a single fit, where we minimize the χ 2 .
In addition, interpolations have the advantage that small fluctuations within the data get smoothed.Polynomial interpolations can have fluctuations, especially higher-order polynomial fits, and hence, we average over different polynomial orders to reduce those fluctuations: nmax n=nmin where w n are normalized, adjustable weights.
For the L20, L26, and L30 lattices, we use equal weights for all polynomial orders.We choose polynomials of orders n = 4, 5, 6, 7 for L20; orders n = 7, 8, 9, 10 for L26; and orders n = 8, 9, 10, 11 for L30.For L40, we use an Akaike average [52] over the orders from 4 to 12, because at some flow times the plateau extraction at larger r gives fluctuating results.This results from underestimated systematic effects and causes a change in the effective polynomial orders, which is considered through the Akaike average.The weights are found analogously to Eqs. ( 30) and (31) through where the final values of w n are fixed by the normalization condition.To propagate the statistical error, the continuum interpolation is done for every jackknife block.This procedure works for most of the data.An exception is the data at small r (up to r/a = 3) for the L26 lattice.We obtain a miscarried interpolation in this case due to large-r effects in the interpolation.For the L26 lattice, it turns out that spline interpolations up to r/a = 3 and changing to the polynomial interpolation for larger r works properly.

Tree-level improvement
To reduce the effects of finite lattice spacings, we apply a tree-level improvement procedure to the data at finite flow time by dividing out the leading lattice perturbation theory result.Following [53], the static energy in lattice perturbation theory at finite flow time can be written as where we have assumed a = 1 for simplicity and D is the lattice propagator: with c w = 0 for the Wilson action, which we use for the simulation part, D MC , and c w = 1/3 for the Symanzik action, which we use in the gradient flow equations, D GF .Similarly to the zero-flow-time case that was derived in [28], the static force coming from the chromoelectric field insertion reduces to a symmetric finite difference: Now, we can tree-level-improve the measured force F meas by dividing out the leading lattice expression and multiplying by the continuum tree-level gradient flow expression where F tree-level cont is defined as the O(α S ) contribution to Eq. ( 15).3: The r 0 scale along the flow-time axis.

Continuum extrapolation
We use the interpolated and tree-level improved data for the continuum limit, which is obtained from extrapolations linear and quadratic in a 2 /t 0 at fixed physical distances r and fixed physical flow times τ F : where m, c l , A, B, and C l are the fit parameters, and c l and C l are the continuum limits of the linear and quadratic extrapolations, respectively.We take an Akaike average, defined for polynomials as in Eqs. ( 41) and (42).Figure 2 shows a working example for the continuum limit at two fixed distances r and for each distance at two different flow times.
Although we use tree-level improved data and Akaike average over linear and quadratic continuum limits, often both χ 2 /d.o.f.'s are too large.Hence, we restrict to data which accomplish that at least one of the extrapolations gives χ 2 /d.o.f.< 3.0 and that the flow radius fulfills √ 8τ F > a.The remaining, filtered data represent reliable continuum limit results, with which we continue the further analyses.

D. r0 and r1 scales
In terms of the force F (r), we define a reference scale as  III: Results for r 0 and r 1 in the zero-flow-time limit at finite lattice spacing.The table includes the statistical error ("stat" subscript), the systematic error by choosing different fit ranges ("AIC" subscript), and the errors added in quadrature.
with a dimensionless number c. Common choices are r 0 (τ F ) ≡ r(c = 1.65, τ F ), r 0 ≡ r 0 (τ F = 0) [6], and r 1 (τ F ) ≡ r(c = 1, τ F ), r 1 ≡ r 1 (τ F = 0) [54].In the original proposals, the scales r(c) are defined with the force at zero flow time, however, in our case, the force is computed at different flow times.Thus, we obtain flow-time-dependent r 0 and r 1 as shown in Figs. 3 and 4. To find r 0 and r 1 , we perform multiple polynomial interpolations of r 2 F (r, τ F ) for larger r values along the r axis and at fixed flow-time ratio τ F /r 2 , and we find the roots of the individual interpolations as (r 2 F ) inter − c = 0.The final scales are given in terms of an Akaike average over the roots of the polynomial interpolators.
Both scales, r 0 and r 1 , approach a constant plateau within a recognizable flow-time regime and start deviating from this plateau at larger flow times.We perform plateau fits within this range for a zero-flow-time extraction of both scales, and we find them comparable to the scales fixed by [4].The results are shown in Table III.We see that the error of the scale setting is dominated by the statistical fluctuation rather than by the Akaike error of the polynomial interpolators.
0.000 0.005 0.010 0.015 0.020 0.025 0.030 FIG.6: Left: An example for the continuum limit of the ratio √ 8t 0 /r 0 at a fixed flow-time ratio.The straight line corresponds to a continuum limit linear in a 2 /t 0 , and the curved line corresponds to a polynomial up to quadratic order in a 2 /t 0 .The red cross and bar shows the final continuum limit of a weighted average of both extrapolations.
Right: The flowed continuum ratio with the constant zero-flow-time limit.The ratio of the scales, r 0 (τ F )/r 1 (τ F ), is shown in Fig. 5.We take the a 2 /t 0 continuum limit followed by a constant zero-flow-time limit.Finally, we obtain r 0 r 1 = 1.380 (14).(50) To our knowledge, there are no prior direct determinations of the scale ratio r 0 /r 1 in pure gauge.However, in Ref. [55] this ratio was roughly estimated to be ∼ 1.37 based on the curves shown in [4], which agrees well with our extraction within errors.The extracted value is about 9% larger than the ratio in unquenched theories with 2+1 or 2+1+1 fermion flavors [1,47].Such a shift is to be expected, due to the effect of unquenching the quark flavors.A similar effect between quenched and unquenched scales has been seen for the gradient-flow scale ratios in Ref. [56].We repeat the same procedure for the ratios √ 8t 0 /r 0 and √ 8t 0 /r 1 .Figure 6 shows an example for the ratio √ 8t 0 /r 0 .The final results for the ratios of the scales are √ 8t 0 r 0 = 0.9569( 66) , (51) Our extracted ratio √ 8t 0 /r 0 agrees within errors with the previous determinations [31,[57][58][59][60][61][62][63],1 albeit the mean value is slightly above most of the existing results.We show the comparison to the existing literature in Fig. 7.The previous FIG.8: Examples of the flowed force in a regime where a constant behavior can be obtained.The left side shows the data for the smallest r where a constant regime can be obtained.The right side shows the flowed force at larger r, we can apply a constant fit in a proper range of τ F /r 2 .
results are somewhat correlated, since most of them focus on t 0 calculation and use the data Refs.[4,5] for at least part of their dataset for r 0 .Again, the quenched ratios are larger than the unquenched ones [1] as has been previously seen in Ref. [56].To our knowledge, this is the first direct measurement of the scale ratio √ 8t 0 /r 1 in pure gauge.

V. ANALYSES OF THE CONTINUUM RESULTS
After having worked out the continuum limit in Sec.IV C, we compare the lattice results with the perturbative expressions to extract Λ 0 .Since gradient flow introduces another scale next to 1/r, we have several possibilities to compare the lattice results with the perturbative expressions.In the first approach, we use the simple expression of the flowed force in Eq. ( 22), which turns out to be also applicable to large-r results.In the second approach, we compare the lattice results either by keeping the scale 1/r fixed and inspecting the behavior along the flow time, or vice versa, by keeping τ F fixed and inspecting the behavior along the distance r.We show plots only for the perturbative one-loop (F1l) and three-loop with leading ultrasoft resummation (F3lLus) expressions; however, in the result tables, we also provide extractions based on the two-loop with and without leading ultrasoft resummation (F2l and F2lLus), and based on three-loop without ultrasoft resummation (F3l) expressions.
A. Constant zero-flow-time limit at fixed r We know from Eq. ( 22) that the static force has a constant behavior at small flow times.Physical quantities are defined at zero flow time; hence, we need to perform the zero-flow-time limit, τ F → 0, while we keep r fixed.In the constant regime, we perform this by a constant fit at fixed distance r. Figure 8 shows data where we obtain a constant behavior of the flowed force.The left side shows the data for the smallest r before the smaller flow time comes into conflict with the √ 8τ F > a condition.The right side shows data at larger r values, where the condition √ 8τ F > a is fulfilled at even minimal flow-time ratios.The small-flow-time expansion is performed in the ratio τ F /r 2 , thus, small flow time is defined in the sense of small flow time ratio, which is a dimensionless quantity.The condition √ 8τ F > a is given in terms of flow times in physical units; hence, considering this condition in terms of ratio moves it to smaller flow-time ratios for larger r, since the r in the denominator decreases the ratio.
In Fig. 9, we see the final result of the constant zero-flow-time limits.We identify an increasing behavior with small errors up to r/ √ t 0 ≈ 3.0.Around the distance of r/ √ t 0 ≈ 2.25, we are faced with difficulty extrapolation to the continuum, and we obtain χ 2 /d.o.f. of order 4 and larger.
We perform Cornell fit to the data from r/ √ t 0 ≈ 1.1 to r/ √ t 0 ≈ 3.0, and we obtain The extracted force from the constant zero-flow-time limits with two Cornell fits.In a window from r/ √ t 0 ≈ 2.1 to r/ √ t 0 ≈ 2.3, the continuum limit turns out to be difficult to reach.We lose the signal at larger r around r/ √ t 0 ≈ 3.0.
The string tension σ is a quantity that is dominated by the large-r regime; hence, in addition, we perform the fit only for data beyond the region where the continuum limits are problematic up to r/ √ t 0 ≈ 3.0.In this case, we obtain σt 0 = 0.154 (6). ( The uncertainty for A is 5 times larger now, which is to be expected, since it is a small r quantity.The results for σ from both ranges agree within their uncertainties.With the result in Eq. ( 51), we can express the string tension in r 0 units: In the past, the string tension was found to be σr 2 0 = 1.353(3) [64] at finite lattice spacing with β = 6.0, which is in good agreement with our result.Continuum results were obtained in [34,65] in another reference scale r.With the ratio √ 8t 0 /r from [34], these results become σt 0 = 0.133(2) and σt 0 = 0.143(2), respectively.Nevertheless, the ratio √ 8t 0 /r in [34] is only an approximation over several lattice sizes, and it is not extrapolated to the continuum limit.Reliable continuum results for the string tension can be found in Refs.[66,67] in units of the critical temperature T c .With the conversion factor T c √ t 0 from [58], these results become σt 0 = 0.1484 (22) and σt 0 = 0.156(3), which agree well with our result.
At small r, we extract Λ 0 from the data by fitting the perturbative force at the available orders.We solve the renormalization group equation for the running of the coupling α s numerically by using the RunDec package [68][69][70], which takes Λ 0 as an input parameter.Setting the scale µ to a specific choice, we remain with a fit function depending solely on the parameter Λ 0 , which serves as the fit parameter.The perturbative coefficients can be found in the literature [19][20][21][22][23][24], and an explicit equation for the force can be found in Eqs.(10) and (11) of [15].Figure 10 shows examples of the fit for two different orders.Table IV shows the fit results.We observe that the error is dominated by the statistical fluctuations rather than by the systematic uncertainty due to the choice of the fit ranges.

B. Flow-time behavior of the force at fixed r or fixed τF
In the very small-r regime, the requirement √ 8τ F > a moves the data points along the τ F /r 2 axis outside the regime where they are flow-time independent even for the smallest possible flow time.Therefore, we compare the lattice data with the full expression of the force, Eq. ( 23) -i.e., without expanding for small τ F .
We inspect the small-r flow-time behavior first at fixed distances r, which corresponds to the classical zero-flow-time limit.In a second approach, we fit along the r axis at fixed flow times to extract √ 8t 0 Λ 0 .Since the dependence on τ F of the numerical extraction turns out to be negligible within the distance and flow-time regions used for the fits to the lattice data, getting √ 8t 0 Λ 0 provides in practice its zero-flow-time limit.
Figures 11 and 12 show the flow-time behavior of the force at two different fixed r values along the flow-time axis.We compare our lattice data with the perturbative expressions of the force at different orders and fit them to the data.Λ 0 serves as the fit parameter.In the figures, we show the fit results for the case of different, but fixed choices    of b in Eq. ( 21).The fit range starts at the smallest possible flow time.For the upper bound, we take an Akaike average [52] over different fit ranges to reduce the systematics due to the fit range choice.
From the figures, we get that the choice of the scale parameter b has a definitive effect on how well the perturbative curve describes the data.A value of b = 1 guarantees the correct asymptotic behavior at large flow time [37].However, we see that b = 1 is the worst of our b choices at describing the actual lattice data in the range of flow times we are interested in.Hence, we use negative values of b, which still lead to valid scaling in our range of flow times, as FIG.14: √ 8t 0 Λ 0 in the zero-flow-time limit from a constant fit including the statistical and fit errors.We compare different orders and different choices of b.The blue points with the label "const zftl" are the results from Table IV.discussed in Appendix A. At small r (left-side plots in Figs.11 and 12), the slope along the flow time is strong, while at the largest r to which we can reasonably apply fixed-r fits (right-side plots in Figs.11 and 12), the data points seem already sensitive to the constant behavior of the force expected at small flow times.At small flow times within our considered flow-time range, the fits with b = 0 and b = −0.5 agree reasonably well with the data.Table V shows an incomplete part of the fit results.
We conclude that for smaller r, which requires going to larger flow-time ratios τ F /r 2 , a negative value of b describes the data better in the large-flow-time regime than b = 0 or 1, whereas for larger r all choices might describe the data within the given uncertainties.That means that fixing b introduces another source of uncertainty which has to be considered.In the zero-flow-time limit, all choices of b give µ = 1/r, which is the natural scale choice for the static force and energy at zero flow time.

Fixed τF
In the next step, we use data at fixed flow times τ F and fit the perturbative force along the r axis.We use an Akaike average [52] over different fit windows for reducing systematics by choosing the right fit window.We perform one-parametric fits at fixed b for √ 8t 0 Λ 0 .Figure 13 shows an example fit for b = 0, 1, −0.5 for F1l and F3lLus at the same flow time.The left vertical line with the dimmer band corresponds to the average lower fit limit for the b = 0  FIG.15: The results for √ 8t 0 Λ 0 .Left: Comparison of all fit results for all methods and for all orders.The inner error bar corresponds to the combined statistical and AIC errors, when applicable.The outer error bar represents the total error, including the systematic uncertainties from the scale variation.The row "large r" shows the results from the first method, where we perform the constant zero-flow-time limit of the force first.Right: Our final result compared to the literature.The points above the dashed line were determined with gradient-flow-based scale settings, t 0 (squares) or w 0 (triangle), while the lower points (circles) are converted from r 0 units with either our ratio in Eq. ( 51) (filled points), or with the ratio from [63] (hollow points).
fit, and the right vertical line to its average upper limit.We observe that the b = 0 fit describes the data over a wide r range from small to larger r.The b = 1 fit describes the data around r/ √ t 0 = 1 and up to larger r in the same way as the b = 0 fit, but it deviates from the data at smaller r in contrast to the b = 0 fit.This matches with its lower fit range limit being at larger-r, but the upper fit range limit being the same as in the b = 0 fit, indicating that the effective fit range for b = 1 is more on the larger r side.The b = −0.5 fit describes the data around r/ √ t 0 = 1 like the b = 0 and b = 1 fits, but in contrast to the b = 1 case, it fits better to the data at small r and deviates from the data at larger r.This matches with its upper fit range limit being at smaller r, but the lower fit range limit being the same as in the b = 0 fit, indicating that the effective fit range for b = −0.5 is more on the smaller-r side.We conclude that the range for b from −0.5 to 1 fits well to the data, but with different effective fit ranges.The b = 0 fit has the widest effective fit range and will serve as our default choice.
Figure 14 shows the fit results for √ 8t 0 Λ 0 at the valid flow-time positions for b = 1, 0, −0.5 at F1l and F3lLus.We observe that for a fixed choice of b, the values of Λ 0 are constant along the flow-time axis.This indicates that the flow-time dependence of the static force at finite flow-time in the distance and flow time ranges explored in this fixed τ F analysis is well captured by a constant one-loop gradient flow correction to the static force.We then extrapolate Λ 0 to the zero-flow-time limit with a constant function.The final results of the constant zero-flow-time limits are shown in Table VI.

C. Estimate of the perturbative systematic uncertainties and final results
Up to this point, we have presented results for Λ 0 with error estimates that include only the statistical errors and the systematic errors from choosing different fit ranges.We still need to include the perturbative uncertainty from the unknown higher-order terms in the perturbative expansion.We can do this by varying the scale µ (21).In previous studies of the static energy [14,15,17], the zero-flow-time scale µ = 1/r was varied by a factor of √ 2. We make here the same choice and vary the parameter s in Eq. ( 21) from s = 0.5 to s = 2.We vary the s parameter only in the zero-flow-time part of Eq. ( 23) and keep it fixed at s = 1 in the finite-flow-time part.In principle, we could vary the scale by a factor of 2 instead of √ 2, but it was noted in Ref. [17] that this requires access to quite small distances r.Our current data do not contain small enough distances to allow for this wider variation.
As already stated in the previous sections, the finite-flow-time part of the static force has a considerable dependence The vertical lines represent the characteristic flow-time window of this study.The left limit on the x axis corresponds to the infinite-flow-time limit, while the right limit corresponds to the zero-flow-time limit.Right: The same ratio, but as a function of τ F /r 2 zoomed around the region of interest.
flow time, since for it, the log(µr) terms vanish.However, this choice does not capture log(µ √ τ F ) terms that become important at large flow time.This is shown by the ratio F L NLO (r, τ F , µ)/F 0 (r, τ F ) becoming large at large flow times for this choice of parameters.The choice s = 0, b = 1 is the natural choice at large flow time, since for it log(µ √ τ F ) terms vanish.However, this choice does not capture log(µr) terms that become important at small flow times.This is shown by the ratio F L NLO (r, τ F , µ)/F 0 (r, τ F ) becoming large at small flow times for this choice of parameters.The choice of s = 1 and b = 1 interpolates between these two extreme and provides a small F L NLO (r, τ F , µ) correction with respect to the leading gradient flow term F 0 (r, τ F ) over the whole range of flow times.Also, the overall scale dependence turns out to be weak with this scale choice.
Nevertheless, our lattice data explore a very specific and limited region of flow-time values, the region in between the vertical lines in Fig. 16.A zoomed-in view of this region is shown in the right plot that makes manifest that different choices of b, keeping s = 1, provide, indeed, even smaller and more stable corrections F L NLO (r, τ F , µ) in the region of interest than the choice b = 1.In particular, this is the case for b = 0 and b = −0.5, which indeed best fit our lattice data, as we have discussed in the main body of the paper.More negative values of b further reduce the relative size of F L NLO (r, τ F , µ), but they make it more scale dependent.Hence, the one-loop expression of the gradient flow expression of the force suggests that for 0.03 ≲ τ F /r 2 ≲ 0.05, the ideal choice of the parameter b is in between 0 and a negative number larger than −1.This is confirmed by the lattice data.Clearly, also, a parametrization with negative b must smoothly go over µ = 1/ √ τ F at large flow time.However, the specific form of the parametrization at large flow times, τ F /r 2 > 0.05, cannot be explored with the present data.autocorrelation times τ int ≲ O(10) for all the relevant observables including the topological charge Q.However, we have less statistics for the L40 lattice, which leads to only 110 configurations per block.Moreover, we observe slower fluctuation of topological charge for the L40 lattice with τ int = 101.6 for the unblocked data.To see whether the block size is large enough also for this ensemble, we plot τ int (Q 2 )3 as a function of the block size in Fig. 19.We observe that even the topological charge has a τ int of order 1, at our chosen block size.More importantly, we also note that the static force has autocorrelation times below 1 for all possible blocksizes.This is expected, as the static energy is known to be moderately unaffected by the topological slowing down [49].Similar topology independence has also been observed for other methods of determining the strong coupling constant [78].
FIG.9:The extracted force from the constant zero-flow-time limits with two Cornell fits.In a window from r/ √ t 0 ≈ 2.1 to r/ √ t 0 ≈ 2.3, the continuum limit turns out to be difficult to reach.We lose the signal at larger r around r/ √ t 0 ≈ 3.0.
FIG.10:The √ 8t 0 Λ 0 extraction with constant zero-flow-times extrapolation at the smallest possible distances r.The vertical lines with the dimmer bar represent the lower and upper averaged fit range according to Eqs.(35) and(36).

2 FF1l 2 FF1lFIG. 11 :
FIG.11: Comparison of different zero-flow-time limits at one-loop order at fixed r.We compare fits with different scale choices obtained by varying b according to Eq. (21).

2 FF3lLus 2 FF3lLusFIG. 12 :
FIG. 12: As in Fig.11, but with the fit function now being the force at three loops with leading ultrasoft resummation.
Fit results at fixed r along the flow-time axis for three different choices of the perturbative order of the static force.We start at the smallest possible flow time and use an Akaike range fit for different upper limits.No value for σ AIC indicates that there are not enough points along the flow time to perform an Akaike range fit.The given χ 2 /d.o.f.corresponds to the most likely fit range.τ F /r 2 Max gives the Akaike averaged upper flow-time limit.
FIG. 13:√ 8t 0 Λ 0 fit at fixed flow time for different orders.The left vertical line corresponds to the average lower fit limit for b = 0, and the right line to its average upper fit limit.The bands around the lines represent the errors of the fit window.

8 F = 1 / r 2 + 8 F = 1 / r 2 8 0. 5 F = 1 / r 2 8 F 8 F = 1 / r 2 + 8 F = 1 / r 2 8 0. 5 F = 1 / r 2 8 FFIG. 16 :
FIG.16: Left: The ratio F L NLO (r, τ F , µ)/F 0 (r, τ F ) [see Eq. (15)] as a function of r/ √ τ F plotted for various values of s and b.The vertical lines represent the characteristic flow-time window of this study.The left limit on the x axis corresponds to the infinite-flow-time limit, while the right limit corresponds to the zero-flow-time limit.Right: The same ratio, but as a function of τ F /r 2 zoomed around the region of interest.

FIG. 19 :
FIG.19: Integrated autocorrelation time τ int as a function of the block size for the largest lattice ensemble L40.The curve with a solid line shows the autocorrelation time for the static force, while the dashed curve shows τ int for the topological charge squared, and the vertical line indicates our chosen blocking.

TABLE I :
Numerical values of the coefficients c n appearing in Eq.

TABLE II :
The parameters for the lattice ensembles.The scale t 0 was set on a smaller subset of lattice configurations.
Examples of continuum extrapolation for two different separations r.The figure shows the linear and quadratic continuum limits: the line corresponds to the linear limit, and the curve to the quadratic one.

TABLE VI :
Results of the √ 8t 0 Λ 0 extraction from the constant zero-flow-time limit of Λ 0 at fixed flow time.The error includes the statistical and the AIC error from the fit.