Determination of $\alpha_s$ from static QCD potential: OPE with renormalon subtraction and Lattice QCD

We determine the strong coupling constant $\alpha_s$ from the static QCD potential by matching a theoretical calculation with a lattice QCD computation. We employ a new theoretical formulation based on operator product expansion, in which renormalons are subtracted from the leading Wilson coefficient. Its prediction for the potential turns out to be valid at the static color charge distance $\Lambda_{\overline{\rm MS}} r \lesssim 0.8$ ($r \lesssim 0.4$ fm), which is significantly larger than ordinary perturbation theory. With lattice data down to $\Lambda_{\overline{\rm MS}} r \sim 0.09$ ($r \sim 0.05$ fm), we perform the matching in a wide region of $r$, which has been difficult in previous determinations of $\alpha_s$ from the potential. Our final result is $\alpha_s(M_Z^2) = 0.1179^{+0.0015}_{-0.0014}$ with 1.3 % accuracy. The dominant uncertainty comes from higher order corrections to the perturbative prediction and can be straightforwardly reduced by simulating finer lattices.


Introduction
Today, facing frontier experiments of particle physics, such as the ones at LHC and Super B Factory (Belle II), there exist increasing demands for more accurate theoretical predictions based on QCD on various phenomena of the strong interaction.Precise determination of the strong coupling constant α s , which is a fundamental parameter of QCD, sets a benchmark for such predictions.In fact, many theoretical developments are required for improving accuracy of α s determination, and once α s is determined, it serves as an input parameter for various predictions.For instance, a precise value of α s will play crucial roles in measurements of Higgs boson properties, in searches for new physics, or in high-precision flavor physics.It is also demanded in the context of precise determination of the top quark mass, predicting running of the Higgs quartic coupling, etc.
Let us quote the current value of α s , given as the world-combined result by the Particle Data Group (PDG), α s (M 2 Z ) = 0.1181 ± 0.0011 [1].Dominant contributions to this value come from determinations by lattice QCD, which have smaller errors than other determinations using more direct experimental inputs.The Flavor Lattice Averaging Group (FLAG) reports an average of lattice determinations as α s (M 2 Z ) = 0.1182 ± 0.0012 [2] based on the studies in Refs.[3][4][5][6][7].The relative accuracies of these current values are 0.9-1.0 per cent.
In determinations of α s by lattice QCD, we need to pay attention to the so-called "window problem," as pointed out in the FLAG report [2].This is a problem that it is difficult to find a wide enough region where both lattice QCD and perturbative QCD predictions are accurate.A lattice simulation is carried out with a finite lattice spacing a, whose inverse plays the role of an ultraviolet (UV) cutoff scale.Hence, the lattice results are accurate in the energy region Q ≪ a −1 .On the other hand, perturbative calculations are accurate at Q 1 GeV(≫ Λ QCD ∼ 300 MeV).Determinations of α s are performed by matching of both results.It turns out that, for currently available lattice cutoff scales, the energy window 1 GeV Q ≪ a −1 cannot be taken widely.
The method of finite volume scheme combined with step-scaling [8][9][10][11] can resolve this problem even at currently available lattice cutoffs.In this method, discretization and finite volume effects are kept under control by a finite volume scheme, while lattice data after the step-scaling running can be matched with perturbation theory at sufficiently high scale.As a result, matching with perturbative prediction can be performed at 10-100 GeV.A recent determination based on this method gives α s (M 2 Z ) with 0.7 per cent relative accuracy [12] (not yet included in the above average values).
In this paper, we determine α s by taking an alternative approach to the window problem: We enlarge the validity range of a theoretical calculation to lower energy where lattice calculations are accurate due to Q ≪ a −1 .For this purpose, we use operator product expansion (OPE) as a theoretical framework.Its difference from perturbative calculations can be stated as follows.Perturbative predictions have an inevitable uncertainty known as renormalon uncertainty, which stems from a certain divergent behavior of perturbative series at large orders.(See [13] for a review of renormalon.)For a dimensionless observable R(Q) with typical scale Q, a renormalon uncertainty is estimated as O((Λ QCD /Q) n ) with a positive integer n (dependent on the observable).In the context of OPE of the same observable, given by the perturbative result is encoded in the leading Wilson coefficient C 1 .In fact, the renormalon uncertainty of C 1 generally has the same order of magnitude as the leading nonperturbative effect (the second term), which corresponds to dim[O 1 ] = n [14].It is expected that the renormalon uncertainty in the leading Wilson coefficient gets canceled when the nonperturbative matrix element is added.Hence, OPE may realize a wider validity range due to the absence of the renormalon uncertainty, in particular at lower energy side.
However, OPE cannot be made a maximal use as long as we naively calculate C 1 in the ordinary perturbation theory.This is because we do not know sufficiently about the nonperturbative matrix element.It is not obvious how to practically eliminate the renormalon uncertainty of C 1 using OPE.In the case where the renormalon uncertainty remains in C 1 , one encounters a difficulty that the nonperturbative effect cannot be estimated using OPE (1.1) since C 1 has an error comparable to this nonperturbative effect.In other words, a renormalon uncertainty causes a mixing between C 1 and the nonperturbative effect.Many studies considering OPE in the literature are not free from such a difficulty.
In Refs.[15,16], a method to cope with a renormalon uncertainty has been proposed.This method enables us to divide C 1 into a renormalon uncertainty and a renormalon free part.By this, we remove a renormalon uncertainty from C 1 before referring to the nonperturbetive matrix element.In this method, we first define C 1 as a UV quantity à la Wilson by introducing an IR cutoff scale µ f (corresponding to a factorization scale of an effective field theory).Then, we separate C 1 (Q2 ; µ f ) into its cutoff independent part and dependent part.While a cutoff dependent part exhibits a connection to the IR physics, a cutoff independent part is regarded as a genuine UV contribution.This cutoff independent part corresponds to a renormalon free part, determined within perturbation theory.Furthermore, by absorbing the cutoff dependent part into the leading nonperturbative matrix element, the nonperturbative matrix element can also be defined as a renormalon free quantity.Hence, we can define the leading Wilson coefficient and the leading nonperturbative effect such that they are clearly separated.This enables us to estimate the nonperturbative effect without suffering directly from the higher order uncertainty of C 1 .
We apply this calculation method for the static QCD potential following Ref.[15].The typical energy scale of the static QCD potential is r −1 , which is the inverse of the distance between the static color charges.OPE of the static QCD potential in r can be performed in the form of the multipole expansion within the effective field theory (EFT), potential non-relativistic QCD (pNRQCD) [17].Thanks to this solid basis, for instance, the so-called u = 3/2 renormalon cancellation 1 has been convincingly shown [17], which gives a solid basis to our OPE formula.We construct a renormalon subtracted Wilson coefficient (which will be denoted by V RF S below) based on the fixed order result, which is currently known up to the next-to-next-to-next-to-leading order (N 3 LO), i.e., O(α 4 s ) [18][19][20][21].A unique feature of our renormalon subtraction is that not only the leading renormalon (at u = 1/2) but also the next-to-leading renormalon (at u = 3/2) is removed from V RF S .In OPE, the leading term is given by V RF S ∼ O(1/r).The NLO term represents the leading nonperturbative effect and is O(r2 ).We include the NLO term with an unknown coefficient which is to be determined by a fit.We will show that the OPE structure can be clearly observed by comparing V RF S with a lattice result: the difference between them behaves as O(r 2 ).Our OPE prediction turns out to be valid up to Λ QCD r 0.8, corresponding to r −1 0.5 GeV.This shows that our theoretical prediction indeed has a wider validity range than the ordinary perturbation theory, which is valid at r −1 1 GeV.We determine α s from the static QCD potential by matching a lattice result with the above OPE where the renormalon uncertainty is subtracted.The lattice results that we use are obtained by the JLQCD collaboration at large cutoffs up to 4.5 GeV [22,23], which facilitate the matching between lattice and OPE calculations.Determinations of α s using the static QCD potential have been performed in Ref. [7] with 3-flavor lattice simulation and in Ref. [24] with 2-flavor lattice simulation.In these studies, perturbative calculations are matched with lattice results in the perturbative regime Λ QCD r 0.2-0.3.Our determination is carried out with the OPE calculation, and the matching range is taken as Λ QCD r 0.6-0.8.
The paper is organized as follows.In Sec. 2, we present our theoretical formula to subtract renormalons in OPE (partially supplemented in Appendix B).In Sec. 3, we determine α s by matching the theoretical calculation with a lattice result.Lattice results and the way to determine α s are also explained therein.Conclusions and discussion are given in Sec. 4. Some referential materials and supplementary arguments are given in Appendices.

Theoretical framework
Our renormalon subtraction formula [15] is constructed based on the EFT, potential nonrelativistic QCD (pNRQCD) [17].This EFT factorizes two typical scales of the static QCD potential.One of the scales is the soft scale ∼ 1/r, which is the inverse of the distance between the static charges.The other is the ultrasoft (US) scale ∼ ∆V (≪ 1/r), which is the energy difference between the octet and singlet bound states.pNRQCD enables us to investigate the US scale physics in a systematic expansion in r∆V ≪ 1.The Lagrangian of this EFT consists of the singlet and octet matter fields and the US gluon fields, while the soft scale contributions are integrated out and encoded in the Wilson coefficients.Our formula is based on the multipole expansion, which expands the static QCD potential in r using the hierarchy r∆V ≪ 1.
In Sec.2.1, we present our formula to subtract renormalons after a brief review of the multipole expansion, on which the formula is based.In Sec.2.2, we explain our treatment of the IR divergence in the three-loop result, which is related to the US scale dynamics.In Sec 2.3, we estimate the higher order perturbative uncertainty of the theoretical calculation, which is required in estimation of systematic errors in α s determination.

Formula to subtract renormalons
Our theoretical calculation is based on the multipole expansion, which gives an expansion of the static QCD potential in r [17]: The explicit r-dependence is given by V S ∼ O(1/r) and δE US ∼ O(r 2 ).(The dots denote higher order terms in r.)The singlet potential V S originates from the soft scale 2 ∼ 1/r (≫ Λ QCD ) and can be evaluated in perturbation theory.In terms of the pNRQCD Lagrangian, V S is a Wilson coefficient.Perturbative result in coordinate space is usually obtained through Fourier transform (FT) of the perturbative evaluation of α V (q 2 ), where the perturbative result of α V (q 2 ) is currently known up to three-loop order [19][20][21]: Here, P n is an n-th order polynomial and we denote its constant part by a n : (2.4) The logarithmic terms in P n can be calculated from the renormalization group (RG) equation and are expressed by a k with k < n and the coefficients of the beta function.δP n represents the IR divergence and associated logarithmic dependence.It is zero for n ≤ 2, and non-zero for n = 3; see Eq. (2.15) for δP 3 .We explain our prescription for regularizing this divergence in the next section 2.2.We collect the explicit expressions for a n in Appendix A.
The NLO term of Eq. (2.1), δE US , is dominantly determined by dynamics of the US scale ∼ ∆V .It is given by a correlator of the US fields in pNRQCD: where E a is the US chromoelectric field; see Ref. [17] for details.Despite the fact that conceptually the singlet potential is a soft quantity, the integration region is usually taken as 0 ≤ q < ∞ as shown in Eq. (2.2).In particular, IR region of the integral is known to cause renormalon uncertainties in V S , and it brings about a mixing between V S and δE US .To avoid this feature, we construct V S as a renormalon free quantity below, following Ref.[15].
We first introduce a factorization (cutoff) scale µ f to divide the energy region as Λ QCD ≪ µ f ≪ 1/r, and define V S as a soft quantity in terms of this cutoff scale: Since all the known renormalons stem from the low energy region of the q-integral, the above definition renders V S free from renormalons. 3In V S (r; µ f ), there is a cutoff dependent part by construction, which is regarded as an IR sensitive part of V S .Such a dependence vanishes only when it is combined with the IR quantities such as δE US .Hence, the mixing is induced through the factorization scale, and in this respect the cutoff dependent part can be regarded as a renormalon related part.In contrast, the cutoff independent part is determined within perturbation theory independently of IR contributions, and hence it can be regarded as a genuine renormalon free part.The renormalon free quantity, which we denote by V RF S (r) [= cutoff independent part of V S (r; µ f )], can be constructed systematically as follows.For α V (q 2 ) in Eq. (2.6), we adopt the next-to-next-to-next-to-leading log (N 3 LL) result, which is obtained by RG improvement of the N 3 LO fixed order result: where α s (q 2 ) is the running coupling constant, namely, the solution to the RG equation at four-loop (for consistency): We solve this equation numerically. 4The three-loop coefficient is originally IR divergent as mentioned and we regularize it with the prescription explained below [a Reg.I or II 3 (q) is given by Eq. (2.20) or (2.21) in Sec.2.2].Related to this divergence, a regularized a 3 has a q-dependence unlike the coefficients up to a 2 .We set n f = 3 and the corresponding light quarks (u, d, s) are treated in the massless approximation in our main analysis.(We consider the finite mass effects in Sec.3.3.)Up to here, the integrand of Eq. (2.6) is determined.Then, by deforming the integration contour of Eq. (2.6) in the complex qplane, we can separate a cutoff independent part from a cutoff dependent part.We explain this formulation explicitly in Appendix B, which is a brief review of Ref. [15].After this procedure, we obtain the following expression: where V RF S (r) is µ f independent and is a renormalon free quantity.The cutoff dependence of V S (r; µ f ) is encoded in the r 0 and r 2 terms (and in further higher order terms), which correspond to the u = 1/2 and u = 3/2 renormalons, respectively.V RF S consists of a Coulomb-like part and a linear part: is expressed with a one-dimensional integral, whose explicit form is given in Eq. (B.9).We evaluate this integral numerically.V C has a Coulomb-like form with logarithmic corrections at short distances.The coefficient of the linear part, C 1 , is proportional to Λ 2 MS , and it can unambiguously be calculated as5 (2.11) One obtains V RF S /Λ MS as a function of Λ MS r without free parameters, where Λ MS is the only dimensionful parameter in massless QCD.Here and hereafter, Λ MS is the Λ-parameter at four-loop in the MS scheme with n f = 3, unless otherwise stated: (See Appendix C for the definition of Λ MS .)We show V RF S (r) in Fig. 1.So far, we have concentrated on the perturbative part V S .Now let us see how the result is combined with the multipole expansion (2.1).Since we define the soft quantity V S with the IR cutoff scale µ f , it is natural to define the US quantity δE US with the UV cutoff scale µ f .Accordingly, the multipole expansion is written as It is confirmed in Ref. [25] that the cutoff dependence of V S (r; µ f ) of the r 2 -term gets canceled against the leading cutoff dependence of δE US (r; µ f ) at the LL level.This corresponds to the u = 3/2 renormalon cancellation, which was first reported in Ref. [17].Although the explicit confirmation at the N 3 LL level (which we consider) is still missing, we assume a parallel scenario.Hence, by using Eq.(2.9), we can perform the multipole expansion as Here, δE RF US is the sum of δE US and C 2 (µ f )r 2 ; hence it is µ f independent and free from renormalons.We omit the constant C 0 (µ f ), which does not have a significant meaning in comparing a theoretical calculation with lattice results; see Sec. 3. Eq. (2.13) is the central formula of our theoretical calculation.The first term is given by Eq. (2.10) and shown in Fig. 1.In our analyses, we regard δE RF US , which is an US quantity, as a nonperturbative object (non-local gluon condensate), and assume δE RF US ∼ Λ 3 QCD r 2 .This is because we focus on relatively long distances where ∆V ≫ Λ QCD is in general not assured.The perturbative result for δE US (obtained within pNRQCD) is used for a limited purpose.Hence, we treat the second term as δE US = A 2 r 2 where A 2 is a fitting parameter, showing the size of the (renormalon free) nonperturbative effect.
Let us state the unique features of V RF S , which is a central object in Eq. (2.13).First, let us clarify the difference from the usual RG improved predictions.Usual N k LL predictions for the static QCD potential are reliable at short distances, but they have an unphysical singularity around r ∼ Λ −1 MS , 6 which distorts the behavior around this region drastically.In contrast, V RF S (r) does not have an unphysical singularity, while N 3 LL accuracy is held at short distances.Therefore, reliable range of V RF S (r) on the low energy side is not limited a priori. 7Secondly, V RF S does not have any renormalons.In particular, it is free not only from the u = 1/2 renormalon but also from the u = 3/2 renormalon, and thus, it is free from the leading r-dependent renormalon uncertainty of O(Λ 3 MS r 2 ).Thanks to these features, V RF S is reliable at short to relatively long distances Λ MS r ∼ O(1).This allows our OPE prediction to have a wide validity range, as will be shown in Sec.3.2.2.

Treatment of US scale
We explain our prescriptions for regularizing the IR divergence in the three-loop coefficient.The IR divergence was first discovered in Ref. [26] and calculated in Ref. [27].In dimensional regularization with D = 4 − 2ǫ, it reads This IR divergence signals breakdown of naive perturbative expansion and is attributed to the dynamics at the US scale.The counterpart of the above divergence is provided from δE US (2.5).Namely, the FT of δE US (r) at O(α 4 s ), defined by is evaluated as [19] δE with

.19)
6 An N k LL prediction is given by where αs(1/r 2 ) is the solution to the (k + 1)-loop RG equation for αs.Due to the singularity of αs(1/r 2 ) around r −1 ∼ Λ MS , the prediction has an unphysical singularity. 7Nevertheless, the naive expectation for the validity range of V RF S is the region up to where the Λ 3 MS r 2 becomes significant, due to the structure of OPE.
= −δP 3 .Hence, the sum of the soft contribution (V S ) and US contribution (δE US ) at order α 4 s is finite.We consider the following two prescriptions.In the first one, we adopt the regularization where we add the US contribution: When RG improvement is applied, we also replace α s (µ 2 ) inside logarithm by α s (q 2 ), in the same way as in Eq. (2.7).In the second prescription, we set a cutoff scale µ US above the US scale to avoid the divergence.In this case, we use This can be compared to Eq. (2.20) by regarding the regulator C A α s as ∆V (r)r ∼ µ US /q.We choose µ US = 3Λ MS or 4Λ MS .The motivation to consider the above two prescriptions in α s determination is to check sensitivity to the treatment of the US contribution.We will see that it does not induce a significant effect.We use the regularization I, given in Eq. (2.20), in our main analysis.

Higher order perturbative uncertainty
V RF S , which we calculate at N 3 LL accuracy, receives higher order corrections.We assume that V RF S can vary to due to higher order corrections.We take δV RF S conservatively as In calculating V RF S | N 2 LL , we use the 4-loop β function8 in evaluating α V (q 2 ), while the fixed order results are used up to a 2 rather than up to a 3 .Hence, δV RF S reduces to the Coulomb+linear part originating from the a 3 -term alone (at N 3 LL).We show the theoretical uncertainty of V RF S in Fig. 2. Higher order uncertainty δV RF S takes a Coulomb+linear form (with log corrections) similarly to V S .Hence, it is qualitatively different from the nonperturbative effect, whose form is quadratic in r.This enables us to estimate the nonperturbative effect while distinguishing it from the perturbative uncertainty. 9It serves for a stable estimation of the nonperturbative effect.We will revisit this point in Sec.3.2.2.
We also remark that the higher order uncertainty is expected to become smaller as the order grows.Since such a property does not hold in the presence of renormalons, this is another non-trivial merit of renormalon subtraction.at N 3 LL (solid) and the higher order uncertainty corresponding to the width of the dotted lines.The r-independent constants are adjusted such that three lines take the same value at Λ MS r = 0.1.

α s determination
We determine the strong coupling constant at the Z boson mass scale, α s (M 2 Z ).We achieve this by matching the theoretical calculation [presented in the previous section in particular in Eq. (2.13)] with a lattice result.We perform two analyses.
The first analysis, which we call Analysis (I), consists of two steps.In the first step, we take the continuum limit of the lattice result (without referring to the theoretical prediction above).In the second step, we compare it with the theoretical prediction to extract α s (M 2 Z ).We proceed while checking (a) if the lattice result can smoothly be extrapolated to the continuum limit, and (b) if V RF S can explain the lattice result up to O(r 2 ) difference, consistently with the OPE structure of Eq. (2.13).After confirming these features, we determine α s (M 2 Z ).In the second analysis, Analysis (II), we perform a global fit to determine α s (M 2 Z ).The two tasks, i.e., an extrapolation to the continuum limit of the lattice data and a determination of α s by comparison with the theoretical prediction, are carried out at once.
Our final result will be adopted from Analysis (II), where we achieve a smaller error than Analysis (I).Analysis (II) is a first-principle analysis, which avoids model interpolating function, required in Analysis (I) for continuum extrapolation.Nevertheless, Analysis (II) is performed without revealing detailed profiles at intermediate steps.To fill the gap, Analysis (I) is performed, where the intermediate steps of the analysis are examined and exhibited explicitly.
We start with an explanation of lattice simulations in Sec.3.1.Subsequently we present Analysis (I) in Sec.3.2 and Analysis (II) in Sec.3.3.Necessary formulas for the analyses are given explicitly in Appendix D.  1. Lattice simulation parameters.For the quark masses m ud and m s , we list bare values in lattice units.The renormalization factor to the MS scheme is available in Ref. [31].

Lattice simulations
Our analysis is performed by using lattice QCD data V latt (r) obtained by the JLQCD collaboration [22,23].Their numerical simulations are carried out in three-flavor QCD in the isospin limit by employing the the Symanzik gauge [28] and Möbius domain-wall quark actions [29].A careful choice of the detailed structure of the quark action reduces the computational cost to simulate fine lattices remarkably while preserving chiral symmetry to good precision [22].Lattice data of V latt are available at three lattice cutoffs, which are determined as a −1 = 2.453(4), 3.610(9) and 4.496 (9) GeV from the Wilson-flow scale [30].
In the following, we denote the three lattice spacings by a 1 , a 2 and a 3 (a 1 > a 2 > a 3 ).Discretization errors of V latt start at O(a 2 ), since chiral symmetry forbids O(a) errors.
The lattice sizes are N 3 s × T = 32 3 × 64, 48 3 × 96 and 64 3 × 128 at a 1 , a 2 and a 3 , respectively.In order to control finite volume effects, their physical sizes are roughly kept constant L = N s a ≈ 2.6 fm, and sufficiently larger than the short distance region r 0.5 fm, where we perform the matching with the OPE.At each cutoff, we take lattice data V latt (r) at a single combination of the degenerate up and down quark mass m ud and the strange quark mass m s .While m s is close to its physical value, m ud corresponds to unphysically heavy pion mass M π ≈ 300 MeV.The correction to V latt (r) due to the unphysical quark masses is taken into account in Analysis (II), but turns out to be small (see Sec. 3.3).Gauge ensembles are generated by using the Hybrid Monte Carlo algorithm.The statistics are 5,000 Molecular Dynamics (MD) time at each simulation point.Simulation parameters are summarized in table 1.
The potential V latt (r) is extracted from the asymptotic behavior of the rectangular Wilson loop where r and t represent its spatial and temporal sizes, respectively.A gauge link smearing [32] is applied to the spatial Wilson lines to suppress excited state contaminations at reasonably small t.The spatial Wilson lines and, hence, the quark pair separation r are chosen to be parallel to the spatial directions (1, 0, 0) and (1, 1, 0), which we call direction d = 1 and 2 in the following.Throughout this study, we estimate the statistical error by the jackknife method.The bin size is chosen as 25 (a 1 ) or 50 MD time (a 2 and a 3 ) by inspecting the bin size dependence of the jackknife error of V latt .The number of bins is N 1 = 200 (a 1 ) or N 2 = N 3 = 100 (a 2 and a 3 ).The statistical correlation is taken into account in the fit (3.1) and subsequent analyses.

Analysis (I): Two-step analysis
In Analysis (I), we first take the continuum limit of the lattice data in Sec 3.2.1.Using the extracted result, we confirm the validity of OPE and compare it with conventional methods in Sec.3.2.2.Matching of OPE with the lattice result to determine α s is performed in Sec.3.2.3.

Continuum extrapolation
We take the continuum limit of the dimensionless combination Here, r 1 is the scale defined by r 2 1 dV dr (r 1 ) = 1.We fix the r-independent constant by subtracting the value at r = r 1 . 10o take the continuum limit X cont latt (r) = X latt (r; a = 0), we first construct sequences {X latt (r; a)} a=a 1 ,a 2 ,a 3 with physical distances r fixed.We choose these reference distances r as the physical points where the coarsest lattice has original data. 11To obtain X latt (r; a) for each a at a desired distance r, we interpolate the lattice data, which are originally discrete, and calculate X latt (r; a) via the interpolating function.An extrapolation to the continuum limit of the sequence {X latt (r; a)} a=a 1 ,a 2 ,a 3 can straightforwardly be performed once the sequence is constructed.
To interpolate the lattice data, we use the following function form: where d = 1, 2 and i = 1, 2, 3 specify the direction and lattice spacing, respectively.The first three terms represent the Cornell potential, which is consistent with the LO perturbation theory at short distances and consistent with the string model at long distances.If we assume this function form to be correct at the continuum limit, correction terms can arise due to finite a and L effects.The fourth term, 1/r 3 -term, is included to take into account the O(a 2 ) discretization error.(Note that the potential has mass dimension one.)The fifth term is similarly introduced for finite L effect to absorb a 1/L 3 -term.Furthermore, the lattice potential data are function of r rather than r, since the rotational symmetry is broken.Therefore, the coefficients can differ depending on the direction.We interpolate the lattice data separately for each direction.This is the reason why the subscript d appears in Eq. (3.2).
In interpolation, we use the lattice data in the range 2a < r < L/2.Namely, we do not use, for instance, the data at r = a.This aims at being free from serious finite a and L effects.(We show in Appendix.E that when we include the data at r = a, continuum extrapolation cannot be performed reasonably.)Indeed, the function form (3.2) is chosen assuming the hierarchy r/a ≫ 1 and r/L ≪ 1.
From the fit, we obtain an interpolating function in a units.We show an example in Fig. 3, where one sees that the interpolating function can indeed fit the lattice data.It is an easy task to calculate X latt (r; a) using the a-unit interpolating function: First, we calculate Inter | d=1,i=1 /d.o.f.= 6.9/(13 − 5).
the ratio r 1 /a from the (slope of) interpolating function, with which one can convert the function into r 1 units.Secondly, we read off a value of X latt (r) at a reference point.By repeating the above procedure for all the jackknife samples, we obtain the average of X latt (r; a) and its statistical error δX latt (r; a i ) for each a at the reference distances r. 12Note that in defining χ 2 of the interpolation [given in Eq. (D.1)], we take the covariance matrix [∆ latt d,i (r k , r j )] common to all the jackknife samples, although the potential V latt,d,i is certainly different for a different jackknife sample.The common matrix is the one calculated from all the samples.Our analysis using the jackknife method proceeds in the same way hereafter.
In table 2, we summarize the fitting parameters in interpolation of Eq. (3.2).Generally, we have smaller statistical errors for finer lattice since the number of the data points used in the interpolation gets larger.The lattice spacings obtained via r 1 /a (using the value r 1 = 0.311(2) fm [33][34][35]) are consistent with the ones determined from the Wilson-flow scale, although the former ones have much larger statistical errors.Now we are in a position to extrapolate the sequences {X latt (r; a)} a=a 1 ,a 2 ,a 3 to the continuum limit a → 0. In Fig. 4, we plot the data point of X latt (r; a) as a function of a 2 /r 2  1 , where we choose r = 3a 1 and r = 8a 1 from the analysis for d = 1.We extrapolate the data by a linear fit in a 2 , in accord with the O(a 2 ) discretization error.Namely, we extrapolate the lattice data at each reference distance to the continuum limit with where γ, δ are the fitting parameters.γ corresponds to X cont latt .In Fig. 4, one can see that the data follow this function and are extrapolated to the continuum limit.To see how smoothly  , which extrapolate the data to the continuum limit.χ 2 ex /d.o.f., which is the reduced χ 2 in this extrapolation, are 1.43 (left) and 0.15 (right).The black data at a = 0 are the extracted continuum limit with the shown statistical errors.
the data are extrapolated to the continuum limit, we show the reduced χ 2 [i.e.χ 2 ex /d.o.f. of Eq. (D.3)] at the reference distances in Fig. 5. Almost all the points are extrapolated to the continuum limit smoothly with the reduced χ 2 less than 2.Only the farthest data of d = 1, which has χ 2 /d.o.f.> 2, is not adopted as our continuum limit result.
In this way, we obtain the continuum limit X cont latt (r).It is shown in Fig. 6.We also list the numerical values in table 3. The covariance matrix for X cont latt (as well as its definition) is presented in Appendix D for the first 6 points in the short distance region.

Consistency checks and comparison with conventional methods
Using the extracted continuum result, before determination of α s , we check consistency of OPE as given in Eq. (2.13).First, we examine the perturbative part, V RF S .We check whether V RF S has a reasonable behavior as we go to higher orders.For this purpose, we construct V RF S at LL to N 2 LL in a parallel way to Sec. 2.1 14 and compare them with the current order prediction at N 3 LL.Note that, at N k LL, the prediction V RF S /Λ /Λ 4-loop MS .We determine these ratios by taking α s (Q 2 ) = 0.2 as an input regardless of the order of the running coupling (following Ref. [15]), which yields  for k = 0, 1, 2. The above condition α s (Q 2 ) = 0.2 assures that different order predictions have no large deviations at Λ 4-loop MS r ∼ 0.0642.This is legitimate since these perturbative predictions should be accurate at such a high energy scale.In Fig. 7, we plot each order prediction in Λ 4-loop MS units, where the lattice result is shown as well.The lattice result X cont latt is converted to Λ MS units from r 1 units by assuming Λ MS = Λ PDG MS ≡ 336 MeV [which corresponds to the current PDG central value of α s (M 2 Z )], and using the central value of r 1 = 0.311(2) fm.In plotting these theoretical predictions, the r-independent  constants are adjusted such that the different order predictions have a common value at Λ 4-loop MS r = 0.0642, and the N 3 LL prediction matches the shortest distance lattice data.From the figure, one can see that the perturbative part, V RF S , gradually approaches the lattice result at higher order. 15ow let us investigate a more detailed issue: we check if the difference between V RF S (at N 3 LL) and the lattice result is O(r 2 ) as OPE dictates.In Fig. 8, we show these two potentials in Λ MS units.The lattice potential is the same as the previous one.For the singlet potential, we add an r-independent constant so that the difference between them is zero at the origin.This constant is determined by a fit assuming that the difference follows const.+const.×r 2 . 16We show their difference by the red boxes.In the difference, a linearlike behavior with an O(Λ 2 MS ) coefficient, which is observed in the lattice and perturbative potentials, vanishes.In fact, they are consistent with an r 2 -term at short distances, as shown by the red line.From this figure, we confirm that the OPE turns out to be valid up to Λ MS r 0.8, corresponding to r 0.5 fm or r −1 0.5 GeV.We remark that this curve is almost unchanged even if we adopt the first 3 points in the fit, although we use the first 6 points in drawing the figure.
To clarify the impact of the above result, in Fig. 9 we compare the validity range of theoretical prediction with conventional ones.We first consider the case adopting the N 3 LL prediction used in the main analysis of Ref. [7], instead of V RF S .The prediction in Ref. [7] has the u = 3/2 renormalon and the unphyiscal singularity at Λ MS r ≃ 0.56 unlike V RF S , although it is free from the u = 1/2 renormalon. 17Due to this singularity, the prediction cannot be obtained at Λ MS r 0.56, and it starts to be distorted around this region as seen from the left panel of Fig. 9 (orange line).In the right panel, the difference from our continuum lattice result is shown (orange points).One cannot observe a const.+const.×r 2 behavior even at Λ MS r 0.6. 18econdly, we consider the fixed order perturbative prediction of V S at N 3 LO.It is free from the u = 1/2 renormalon (once a value at some distance is subtracted) and the unphysical singularity, while it has the u = 3/2 renormalon.Since it is a fixed-order potential, the prediction is reliable only around the region r ∼ µ −1 .We choose µ as Λ MS /µ = 0.4, where µ −1 is close to the smallest r among the lattice data points in the continuum limit.This also fixes the value of α s (µ 2 ) as α s (µ 2 ) = 0.59.In the right panel of Fig. 9, the difference from the lattice data is shown (dark blue points).We can observe the OPE structure up to a certain distance region: the first three points (Λ MS r 0.55) can be fitted reasonably by a const.+const.×r 2 function, while the first six points (Λ MS r 0.8) cannot be.However, we note that the result of the analysis is sensitive to a choice of the renormalization scale.If we make the renormalization scale twice (Λ MS /µ = 0.2), the range that OPE is applicable and the coefficient of an r 2 -term vary considerably.(We cannot take the scale 1/2, since the running coupling constant diverges above this scale.)This indicates that the OPE structure is not held stable against the higher order correction.In contrast, if we perform a parallel analysis using V RF S , 19 we always confirm that OPE is valid up to Λ MS r 0.8 and the variation of the coefficient of an r 2 -term is milder.Namely, the OPE structure is stably observed.This allows us to estimate the nonperturbative effect (coefficient of an r 2 -term) in a reasonable and reliable way.
The above arguments show that our theoretical calculation indeed allows us to use a wider range for matching than previous studies.We confirmed that the OPE structure is observed up to Λ MS r 0.8.This is achieved thanks to a stable and reliable prediction of V RF S at short to relatively long distances.This feature originates from the RG improvement, the absence of the unphysical singularity, and the u = 3/2 renormalon subtraction.The latter two features result from subtraction of IR contributions [see Eq. (2.6)] in constructing V RF S .This subtraction removes instability caused by IR dynamics.

α s determination: Matching between OPE and lattice result
We now determine α s (M 2 Z ) by matching the lattice result with OPE.Our determination of α s reduces to the problem to find an appropriate x = Λ MS r 1 such that OPE agrees with the lattice result.Once x is determined, we obtain Λ MS through the value r 1 = 0.311(2) fm.Then, we obtain α s (M 2 Z ) by solving the RG equation for α s (µ 2 ).We compare the lattice and theoretical potentials in Λ MS units by converting the lattice result to Λ MS units with x: The OPE prediction is given by where A 0 and A 2 are the fitting parameters. 20In the matching, we choose the lattice data points satisfying Λ PDG MS r < 0.8 in order for OPE to be valid.Hence, the first six points of Fig. 6 are used.The correlation matrix required in this analysis is presented in Appendix D. By a fit, we obtain where the statistical error is estimated by the jackknife method.The reduced χ 2 is χ 2 match /d.o.f.≈ 0.4/(6 − 3) [see Eq. (D.5) for definition of χ 2 match ].From Eq. (3.6), we obtain using r 1 = 0.311 fm.The size of the nonperturbative effect is estimated as The obtained 3-flavor Λ MS of Eq. (3.7) gives the 5-flavor coupling as This value is obtained as follows.First, we calculate α s (µ 2 )| n f =3 below the charm MS mass µ < m c = 1.3 GeV, from the obtained Λ MS using Eq.(C.1) in Appendix C. Secondly, we obtain the 4-flavor coupling at the charm MS mass (which we take as a matching scale) µ = m c , using the 3-loop matching equation [37].Then, we obtain α s (µ 2 ) for m c < µ < m b = 4.2 GeV by solving the RG equation for n f = 4. Similarly, we obtain the 5-flavor coupling at the bottom MS mass by the matching equation.Then, we obtain the coupling at the Z boson mass M Z = 91.187GeV, α s (M 2 Z ), by solving the RG equation for n f = 5.We solve the RG equation for α s (µ 2 ) numerically; see footnote 4.
For convenience, we summarize the conditions used in our main analysis, with which we determine the central value of α s (M 2 Z ).
• Controlling finite a and L effects: Lattice data in the range 2a < r < L/2 are used in interpolation • Interpolating function: Cornell type potential [Eq. ( • Extracted quantity: • Matching range (Used lattice data in the continuum limit): Λ PDG MS r < 0.8 • Conversion of x to Λ MS : Central value of r 1 = 0.311(2) fm We now estimate systematic errors of our determination.For this purpose, we perform re-analyses by changing the conditions as follows and examine variations of determined α s (M 2 Z ).
• Finite a effects: We use the lattice data of a < r < L/2 in interpolation such that the shorter distance points r a are included, although we still omit the data at r = a. 21 Interpolating function: The Cornell potential has a defect that it does not have a logarithmic correction in the Coulomb part at short distances where it should be 1/(r log(rΛ MS )) rather than 1/r.Such a logarithmic correction follows from the oneloop β function.We replace the Coulomb part by the one consistent with the one-loop β function: where V large-β 0 C (r) is the Coulomb potential calculated in the large-β 0 approximation according to the method of Ref. [15] or Appendix B. Its asymptotic form is given by 22 • Subtraction point: We take the continuum limit of where the subtraction point is changed.
• Higher order uncertainty: We replace V RF S in matching as with t = −1 or 1 in order to estimate higher order uncertainty; see Eq. (2.23) for δV RF S .
• US regularization: We adopt the regularization prescription II, given by Eq. (2.21).We choose µ US as 3Λ MS and 4Λ MS .
• Matching range: We vary the range of the lattice result used in the matching as to examine the stability of OPE truncated at O(r 2 ).
The estimated systematic errors are summarized in table 4. By taking the root-sum-square of the errors, we obtain from Analysis (I). 22In interpolating the lattice unit potential with the above fitting form, we introduce y = Λ 1-loop MS a as the fitting parameter in order to convert V large-β 0 C /Λ 1-loop MS to a units as

Analysis (II): Global fit
In Analysis (I), an interpolating function is assumed in order to take the continuum limit of the potential, although the exact form is unknown.This is a short-coming from the viewpoint of first principles.In Analysis (II), we perform a first-principle determination, without using such a model-like interpolating function.This is achieved by a global fit in which the continuum extrapolation and the matching with a theoretical calculation are performed at once.This analysis is based on the idea that the OPE prediction should be correct at short distances and coincide with the lattice data once the discretization errors are removed.Then, OPE is matched with the modified lattice data which can be regarded as the result in the continuum limit: Discretization errors contained in the original lattice data V latt,d,i are removed by the second and third terms (depending on i and d), and the last term adjusts the r-independent constant; 1 r is the LO result in the lattice perturbation theory, which deviates from a smooth 1/r-function due to finite a and L effects.Hence, the second term removes the discretization error at the tree-level.Note that the tree-level process is given by a one-gluon exchanging diagram and is order α s .Here, κ is regarded as an effective coupling of lattice perturbation theory, and is treated as a fitting parameter.The third term extrapolates the data to the continuum limit by removing the remaining error of order α 2 s a 2 .We perform matching by converting lattice and theoretical potentials to GeV units.Lattice data are converted to this unit using a's estimated by the Wilson-flow scale.The theoretical potential is converted with z = Λ MS [GeV], which is unknown in advance and thus is treated as a fitting parameter.Therefore, an OPE prediction used here is given by (3.18) Since an r-independent constant is already included in Eq. (3.17), it is not included here.
In matching, we adopt the lattice data in the range rΛ PDG MS < 0.6.Here, we choose a shorter distance region than in Analysis (I) since we have more available data points.In this analysis, we do not omit short distance data at r ∼ a, and in particular we include the data at r = a.(Note that the continuum extrapolation cannot be taken reasonably in the case including the data at r = a in Analysis (I), as discussed in Appendix E.) Thus, we take into account the tree-level correction, which is powerful to remove the discretization error at short distances (where perturbation theory works) and does not need the hierarchy r ≫ a.The number of the i-th lattice data used in the matching is 7, 10, 13 points for i = 1, 2, 3, respectively.
In this analysis, we determine 16 parameters in total: Λ MS , A 2 , six tree-level correction parameters κ's, two f 's, and six r-independent constants c 0 's.Due to the nature of this global fit, the lattice result in the continuum limit is determined such that it matches with the OPE prediction.In this respect, the continuum extrapolation is not taken within lattice i (size) simulation, but it is constrained by the OPE prediction. 23Thus, the lattice data in the continuum limit in Analyses (I) and (II) have qualitatively different meanings.
In this global fit, we obtain where the statistical error is estimated by the jackknife method.We summarize the other parameters in this global fit in table 5.The reduced χ 2 of this fit is χ 2 GF /d.o.f.= 8.7/(30 − 16) [see Eq. (D.7) for definition of χ 2 GF ], showing the validity of the analysis.f 's are consistent with zero, which suggests that the discretization error is quite small after the tree-level correction is taken into account.
To check if the tree-level correction works in a reasonable way, we show the determined values of κ's in Fig. 10.In this figure, we compare κ d,i with its naively expected value, C F α s (µ 2 ) taking the renormalization scale as µ = a −1 i .Note that C F = 4/3 is multiplied since the LO result in the continuum theory is V QCD (r)| tree = −C F α s /r.In plotting the running coupling, we assume Λ MS = Λ PDG MS and n f = 3.The determined κ's are consistent with the naively expected values within the statistical errors, which supports validity of our analysis.Large statistical errors for κ d=2,i stem from the small number of data for d = 2.
We show the lattice result in the continuum limit [Eq.(3.17)] and the OPE prediction [Eq.(3.18)] which are determined by the fit in Fig. 11.From the figure, one can see that the analysis is performed reasonably, and that the OPE calculation and the lattice result are mutually consistent in the examined region.
The obtained Λ MS in Eq. (3.7) gives The procedure to obtain α s (M 2 Z ) is the same as for Analysis (I).For convenience, we summarize the conditions used in our main analysis, with which we determine the central value of α s (M 2 Z ).
• Controlling finite a effects: The data at r ≥ a are used combined with the tree-level correction.For this V pt , we take µ = 3 GeV and α s = 0.25.In this way, we estimate both theoretical prediction and lattice result at the physical point.
• Matching range: We vary the range of the lattice result used in the matching as to examine the stability of OPE truncated at O(r 2 ).
• Factorization scheme: In extracting the renormalon free part V RF S , we rewrite the integrand of V S by a complex function; see (B.2) in Appendix B. In general, there can be other choices for this function, and in this regard, we have chosen a certain scheme.A different scheme practically causes an O(r 3 ) difference in the OPE prediction truncated at O(r 2 ); see Ref. [16] for details. 25 To see an effect of this scheme dependence, we add an A 3 r 3 -term in the fit so that this scheme dependence is absorbed.(Note that, in order to determine coefficients up to higher orders in r, a wider fitting range is required.We choose the range in this analysis as Λ PDG MS r < 0.8, where A 2 and A 3 are stable against variation of the range. 26) The estimated systematic errors are summarized in table 6.Some error sources included in Analysis (I) are absent thanks to the first-principle nature of this analysis.In addition, most of the systematic errors are reduced compared to Analysis (I).In particular, the higher order uncertainty is smaller since a shorter distance region is used; see Fig. 2. The mass effects turn out to be negligibly small even if we consider the constituent quark mass.This is because we are probing a sufficiently short-distance region.

Summary of results
We have performed two analyses.In Analysis (I), we first took the continuum limit of the lattice data, and then we matched the result with the OPE prediction.Although this analysis partially relies on a model-like assumption, we explicitly showed that (a) the our another analysis [Analysis (I)], where we examined intermediate processes step by step.
Although the energy region extends to lower energy side than conventional determinations using perturbative calculation, it is confirmed that varying the energy range does not induce a significant systematic error.Dominant error of our determination comes from systematic errors, in particular from the higher order perturbative uncertainty of the leading Wilson coefficient.We emphasize that a finer lattice simulation will straightforwardly reduce this error, since we can adopt a shorter distance range in the fit, where the uncertainty becomes smaller. 27,28 We believe that our analyses are useful in not only determining α s but also in promoting understanding on the OPE structure and lattice discretization errors.For instance, we observed that the difference between the renormalon subtracted Wilson coefficient and the lattice result behaves as O(r 2 ), as indicated from OPE.This is a first numerical evidence of the OPE structure, as predicted by pNRQCD, for the static QCD potential in this distance range.Concerning the lattice discretization error, we clarified that (i) the data at r = a indeed has a serious finite a effect (Appendix E), and (ii) once the tree-level correction is considered combined with the OPE calculation, the finite a effect can be removed with reasonable values of lattice effective couplings.The contours C a and C b are displayed in Fig. 13.The integral along C a is clearly independent of µ f .Although the integral along C b looks µ f dependent, it contains a µ f -independent part.We evaluate this integral as In expansion of the exponential factor, the real and pure imaginary coefficients appear in turn.
The terms with real coefficients satisfy the relation {f (z)} * = f (z * ).Owing to this, these parts can be calculated as where C Λ QCD is shown in Fig. 14.The coefficients C −1 and C 1 are µ f independent and real.Numerical evaluation of these coefficients is sufficient for our purpose.[C 1 is given by Eq. (2.11).]We remark that the analytical results up to N 2 LL can be found in Ref. [15].On the other hand, the terms with imaginary coefficients do not satisfy the relation {f (z)} * = f (z * ), and the above deformation cannot be applied.Therefore, we have where µ f dependence remains.Based on the above argument, we can construct a µ f -independent quantity V RF S .Note that a µ f -independent part is also given by the integral along C a .Then by collecting all the µ f -independent part, we obtain dq q e −qr Im α V (−q 2 + i0) − C −1 .(B.9) In the last line, we rotate the contour C a to the line along e iπ/2 q with real positive q.
Once the µ f -dependent part of V S (r; µ f ) is considered as well, one obtains the decomposition shown in Eq. (2.9).

C Definition of Λ MS
The definition of the scale Λ in the MS scheme, Λ MS , is given by where α s represents the coupling at the renormalization scale µ.We approximate the β function at four-loop as in Eq. (2.8), which gives the definition of Λ 4-loop MS , used extensively in this paper.

D χ 2 and covariance matrix
We present definitions of χ 2 and covariance matrices used in our analyses, which may be useful especially for non-expert readers.Interpolation [Analysis (I)] We define χ 2 in the interpolation with a covariance matrix as  latt,d,i (r) is defined in Eq. (3.2) and k, l run over the lattice points under consideration.The covariance matrix ∆ latt (r k , r l ) is calculated as in the jackknife method, where N i is the number of bins for the i-th lattice simulation; see table 1.If the subscript d is shown, it expresses a covariance matrix among the potentials V latt,i,d . 30ontinuum extrapolation [Analysis (I)] χ 2 in the extrapolation to the continuum limit is defined as where Q(a) is defined by Eq. (3.3).The covariance matrix for X cont latt is calculated as ∆ cont (r i , r j ) = (N tot − 1) (X cont latt (r i ) − X cont latt (r i ) ) • (X cont latt (r j ) − X cont latt (r j ) ) .(D.4) Note that in the continuum extrapolation, the jackknife samples with the size N tot = 3 i=1 N i = 400 are generated since we have three independent lattice measurements.We present the numerical result of ∆ cont in table 7.
Global fit [Analysis (II)] We define χ 2 in the global fit in Analysis (II) as Table 7. Covariance matrix for X cont latt , ∆ cont (r i , r j ).The first row is r i /r 1 and the fist column is r j /r 1 .The (i, j) component is the numerical value of ∆ cont (r i , r j ).Note that ∆ cont (r i , r j ) is a symmetric matrix, and hence, we only show the elements of the upper triangular part.
Here, the covariance matrix consists of three matrices of dimension 7, 10, and 13 in a block diagonal form: where the definition of each matrix is given by Eq. (D.2).
E Case including data at r = a in Analysis (I) In Analysis (I), we do not use the data at r = a in interpolating lattice data in our analyses, in order to suppress serious finite a effects.Here, let us see what happens if we include this shortest point.We interpolate lattice data including the ones at r = a, and obtain X latt (r; a) in the same way.In Fig. 15, we plot the data points of X latt (r; a) taking the horizontal axis as (a/r 1 ) 2 .One can see that they do not obey linear behaviors in a 2 .We remark that even the data for r = 7a 1 , where the finite a effect is considered to be well suppressed, cannot smoothly be extrapolated to the continuum limit.It shows that the data at r = a, which has a small statistical error, dominantly contributes to determining the interpolating function, and thus, the interpolating function is seriously distorted.In Fig. 16, we show χ 2 ex /d.o.f. in this case, corresponding to Fig. 4. Extrapolations to the continuum limit do not work for d = 1.For d = 2, where the shortest point is located at r = √ 2a, extrapolations work better than for d = 1.We conclude that the data at r = a has serious discretization error, and we should be cautious about treating it.

Figure 1 .
Figure 1.Renormalon free singlet potential V RF S /Λ MS as a function of Λ MS r (black solid line).V C (r)/Λ MS is shown by the blue dotted line, and the linear contribution C 1 r/Λ MS [Eq.(2.11)] is shown by the green dashed line.The results are obtained with regularization method I [Eq.(2.20)].

Figure 4 .
Figure 4. X latt (r; a) as functions of (a/r 1 ) 2 .We show them for r = 3a 1 (left) and r = 8a 1 (right), which appear in the analysis for d = 1, as examples.Black lines are linear functions in a 2 , which extrapolate the data to the continuum limit.χ 2 ex /d.o.f., which is the reduced χ 2 in this extrapolation, are 1.43 (left) and 0.15 (right).The black data at a = 0 are the extracted continuum limit with the shown statistical errors.

Figure 5 .
Figure 5.The reduced χ 2 in extrapolation [see Eq. (D.3) for definition] for the distances where the continuum limit are taken.As a benchmark, χ 2 /d.o.f.= 2 is shown by the red line.

Figure 6 .
Figure 6.Continuum limit of the lattice result, X cont latt .Blue points originate from d = 1 and red ones from d = 2.
Therefore, in order to plot the k-th order prediction in Λ 4-loop MS units, we need a conversion parameter Λ (k+1)-loop MS

Figure 7 .
Figure 7.Comparison of the lattice result X cont latt (black dots) with V RF Sat LL (red), NLL (green), N 2 LL (blue), and N 3 LL(purple), for the inputs n f = 3 and α s (Q 2 ) = 0.2.Λ PDG MS = 336 GeV is used to convert X cont latt to Λ MS units.r-independent constant of each potential is adjusted.

Figure 8 .
Figure 8.Comparison of the lattice result (cont.limit: blue circles) and leading OPE prediction (V RF S /Λ MS : blue line) using Λ PDG MS and adjusting r-independent part.The difference (red boxes) is fitted by const.×r 2 (red line) at small r.

Figure 9 .
Figure 9. (Left) Static potentials obtained from lattice (black), V RF S /Λ MS (purple), the N 3 LL prediction of Ref. [7] (orange), and the fixed order N 3 LO prediction with Λ MS /µ = 0.4 (dark blue).(Right) Differences between the lattice result and the theoretical predictions.Purple data represent the difference from V RF S , the orange ones from the N 3 LL prediction of Ref. [7], and the dark blue ones from the N 3 LO prediction.The curves in this figure are const.+const.×r 2 functions determined by fits.The purple solid line is determined with the first six points [χ 2 /d.o.f.= 2.4/(6 − 2)], and the purple dotted one with the first three points [χ 2 /d.o.f.= 0.19/(3 − 2)].The orange solid line is determined with the first three points [χ 2 /d.o.f.= 85/(3 − 2)], and the orange dotted one with the first two points [d.o.f.= 0].The dark blue solid line is determined with the first six points [χ 2 /d.o.f.= 141/(6 − 2)], and the dark blue dotted one with the first three points [χ 2 /d.o.f.= 0.002/(3 − 2)].

8 ±1Table 4 .
Estimates of systematic errors in Analysis (I) from variations of the central value of α s (M 2 Z ) in units of 10 −4 when varying the analysis conditions.In the upper row, variations are shown.(Detailed conditions are shown inside brackets).Assigned systematic errors are shown in the lower row.

Figure 11 .Table 6 .
Figure 11.Lattice result in the continuum limit (black points) and the OPE calculation (green) determined simultaneously by the fit in Analysis (II).The distance region used in this fit rΛ PDG MS < 0.6 is shown by the dotted line.For reference, r = r 1 is also shown.

Figure 13 .
Figure 13.Contour C a and C b in the complex q-plane.

Figure 15 .
Figure 15.X latt (r; a) as functions of (a/r 1 ) 2 when we include r = a in interpolation.We show them for r = a 1 (left) and r = 7a 1 (right), which are reference distances for d = 1.Black lines are linear functions in a 2 determined by fittings.χ 2 ex /d.o.f., which is the reduced χ 2 in this extrapolation, are 20 (left) and 4.3 (right).

Table 2 .
Fitting parameters in Eq. (3.2), r 1 /a, and the reduced χ 2 .The fitting parameters are shown as dimensionful quantities (except for α), which are originally obtained as dimensionless parameters normalized by proper powers of a.To make them dimensionful parameters, we use the Wilson-flow scale a. N i,d is the number of the lattice data used in interpolation.

Table 3 .
Numerical results of X cont latt (r).
0.0648, and Λ 4-loop MS /Q = 0.0642.By regarding Q as a common scale, we obtain the ratios Λ

Table 5 .
Fitting parameters in Analysis (II).N i,d expresses the number of data points for direction d of the i-th lattice.