Determination of α s from static QCD potential with renormalon subtraction

We determine the strong coupling constant αs(MZ) from the static QCD potential by matching a lattice result and a theoretical calculation. We use a new theoretical framework based on operator product expansion (OPE), where renormalons are subtracted from the leading Wilson coefficient. We find that our OPE prediction can explain the lattice data at ΛQCDr . 0.8. This allows us to use a larger window in matching, which leads to a more reliable determination. We obtain αs(MZ) = 0.1179 +0.0015 −0.0014. Today, facing frontier experiments of particle physics such as those at LHC and super B Factory, 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 sets a benchmark for such predictions. For instance, a precise value of αs will play crucial roles in measurements of Higgs boson properties, in searches for new physics, in high-precision flavor physics, etc. The current value of αs, given as the world-combined result by the Particle Data Group (PDG), reads αs(MZ) = 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. Nevertheless, most lattice QCD determinations have the “window problem” in an explicit or implicit way, as pointed out in the Flavor Lattice Averaging Group (FLAG) report [2]: It is difficult to find a region where both lattice QCD and perturbative QCD predictions are accurate. At short distances (Q ≫ ΛQCD), where perturbation theory is accurate, lattice data are distorted by ultraviolet (UV) cutoff effects due to the finite lattice spacing a, whereas at larger distances (Q ∼ ΛQCD), where finite a effects are suppressed, perturbation theory is not reliable. The method of the finite volume scheme combined with step-scaling was proposed to solve this problem [3, 4]. This method enlarges reliable energy region of lattice simulation. As a result, matching with perturbative prediction can be taken in a wide range at high energy, 10–100 GeV. In this Letter, we propose an alternative approach to the window problem: We enlarge the validity range of theoretical prediction to the region where lattice calculations are accurate, Q ≪ a. To this end we use operator product expansion (OPE) with Preprint submitted to Physics Letter B May 25, 2019 subtraction of renormalons. Accuracy of a perturbative prediction has a limitation due to renormalons (which specify certain divergent behaviors of perturbative series), and an O(ΛQCD/Q )-error is inevitable, for a dimensionless observable with typical scale Q. In OPE, the O(ΛnQCD/Q )-term is described by a nonperturbative matrix element (ME). Hence, we can enlarge validity range of theoretical prediction to lower energy by subtracting renormalons appropriately from perturbative prediction in the framework of OPE. Although OPE is a good and well-known framework, there is a difficulty in practical calculations. There has not been an established way to factorize the two components of OPE, Wilson coefficients and nonperturbative MEs, which are conceptually UV and infrared (IR) quantities, respectively. Although one may find in the literature that Wilson coefficients are calculated in usual perturbation theory, this procedure is not desirable since loop integrals in dimensional regularization contain both UV and IR modes. In particular, IR contributions cause renormalon uncertainties in a Wilson coefficient, which makes it practically impossible to distinguish a nonperturbative ME from the renormalon uncertainties. In Refs. [5, 6], a formulation to separate UV and IR contributions in OPE has been proposed. A Wilson coefficient is constructed as a UV quantity, free from renormalon uncertainties (i.e., renormalons are subtracted from the Wilson coefficient). The formulation concurrently defines a nonperturbative ME as an IR quantity. This prevents mixing of the nonperturbative ME with renormalon uncertainties in the Wilson coefficient, thus enabling us to perform OPE in an ideal way. In particular, for the static QCD potential, one can calculate a Wilson coefficient systematically from the fixed-order perturbative result. It is identified with the leading term of OPE in the solid framework of potential nonrelativistic QCD (pNRQCD) effective field theory [7]. We determine αs from the QCD potential by matching a lattice result and OPE. We can take the matching range down to relatively low energy scale ΛQCDr . 0.6–0.8 by subtraction of renormalons. This is in contrast to previous αs determinations using the QCD potential [8, 9], in which lattice results are matched with perturbative results in the region ΛQCDr . 0.2–0.3. We use the lattice result at cutoffs up to 4.5 GeV obtained by the JLQCD collaboration. Our theoretical calculation is based on OPE with renormalon subtraction and the next-to-next-to-next-to-leading order (NLO) result of perturbation theory. The unique feature of our method is to perform OPE avoiding the mixing of a Wilson coefficient and a nonperturbative ME. This clarifies their respective roles, and an estimate of theoretical error can be given clearly. In contrast, in many studies considering OPE, an estimate of perturbative error cannot be distinguished from an estimate of nonperturbative effects, since they are mixed in a naive calculation method. We will show that our OPE prediction can explain lattice data at r & 0.5 GeV (or r . 0.4 fm), where usual perturbation theory cannot work sufficiently. Our theoretical prediction for the QCD potential is based on multipole expansion within pNRQCD, which is an OPE in ~r. The QCD potential is expanded as [7] VQCD(r) = VS(r) + δEUS(r) + . . . , (1) where the explicit r dependence of each term is VS(r) ∼ 1 r and δEUS(r) ∼ r , and the dots denote the higher-order terms in r. (We suppress the r-independent part.) In the

subtraction of renormalons.Accuracy of a perturbative prediction has a limitation due to renormalons (which specify certain divergent behaviors of perturbative series), and an O(Λ n QCD /Q n )-error is inevitable, for a dimensionless observable with typical scale Q.In OPE, the O(Λ n QCD /Q n )-term is described by a nonperturbative matrix element (ME).Hence, we can enlarge validity range of theoretical prediction to lower energy by subtracting renormalons appropriately from perturbative prediction in the framework of OPE.
Although OPE is a good and well-known framework, there is a difficulty in practical calculations.There has not been an established way to factorize the two components of OPE, Wilson coefficients and nonperturbative MEs, which are conceptually UV and infrared (IR) quantities, respectively.Although one may find in the literature that Wilson coefficients are calculated in usual perturbation theory, this procedure is not desirable since loop integrals in dimensional regularization contain both UV and IR modes.In particular, IR contributions cause renormalon uncertainties in a Wilson coefficient, which makes it practically impossible to distinguish a nonperturbative ME from the renormalon uncertainties.
In Refs.[5,6], a formulation to separate UV and IR contributions in OPE has been proposed.A Wilson coefficient is constructed as a UV quantity, free from renormalon uncertainties (i.e., renormalons are subtracted from the Wilson coefficient).The formulation concurrently defines a nonperturbative ME as an IR quantity.This prevents mixing of the nonperturbative ME with renormalon uncertainties in the Wilson coefficient, thus enabling us to perform OPE in an ideal way.In particular, for the static QCD potential, one can calculate a Wilson coefficient systematically from the fixed-order perturbative result.It is identified with the leading term of OPE in the solid framework of potential nonrelativistic QCD (pNRQCD) effective field theory [7].
We determine α s from the QCD potential by matching a lattice result and OPE.We can take the matching range down to relatively low energy scale Λ QCD r 0.6-0.8 by subtraction of renormalons.This is in contrast to previous α s determinations using the QCD potential [8,9], in which lattice results are matched with perturbative results in the region Λ QCD r 0.2-0.3.
We use the lattice result at cutoffs up to 4.5 GeV obtained by the JLQCD collaboration.Our theoretical calculation is based on OPE with renormalon subtraction and the next-to-next-to-next-to-leading order (N 3 LO) result of perturbation theory.The unique feature of our method is to perform OPE avoiding the mixing of a Wilson coefficient and a nonperturbative ME.This clarifies their respective roles, and an estimate of theoretical error can be given clearly.In contrast, in many studies considering OPE, an estimate of perturbative error cannot be distinguished from an estimate of nonperturbative effects, since they are mixed in a naive calculation method.We will show that our OPE prediction can explain lattice data at r −1 0.5 GeV (or r 0.4 fm), where usual perturbation theory cannot work sufficiently.
Our theoretical prediction for the QCD potential is based on multipole expansion within pNRQCD, which is an OPE in r.The QCD potential is expanded as [7] V QCD (r) = V S (r) + δE US (r) + . . ., where the explicit r dependence of each term is V S (r) ∼ 1 r and δE US (r) ∼ r 2 , and the dots denote the higher-order terms in r. (We suppress the r-independent part.)In the following we consider the first two terms of OPE, shown explicitly in Eq. (1), unless stated otherwise.While the singlet potential V S is a UV quantity, δE US and higher correction terms are dominantly IR quantities determined by nonperturbative dynamics.In usual perturbative evaluation of V S , renormalon uncertainties appear, whose leading r-dependent uncertainty is O(Λ 3 QCD r 2 ).Ref. [7] has pointed out that this renormalon uncertainty is canceled against that of δE US .This observation suggests that one should subtract renormalons from V S to define it as an unambiguous object and also to make δE US free from renormalons.
We subtract renormalon uncertainties of V S as follows [5,6].The QCD potential is formally given by In this expression q varies from 0 to ∞.Since the singlet potential corresponds to UV part of V QCD , we define with a factorization scale µ f .V S does not have renormalon uncertainties since IR contributions are removed. 1In V S , µ f -dependent part is sensitive to IR dynamics, and when combined with δE US , it becomes independent of µ f [up to O(r 2 )].In contrast, µ f -independent part of V S corresponds to a pure UV contribution, which is accurately predictable within perturbation theory.We construct a µ f -independent part of V S (r; µ f ), denoted as V RF S (r), in the following manner.We utilize the perturbative result for α V (q) known up to O(α 4 s ) (N 3 LO) [10,11,12].We improve the fixed-order result by renormalization group (RG) using the 4-loop β-function.Up to here, the integrand of Eq. ( 3) is determined.Then, by deforming the integral path in the complex-q plane, we can extract a µ f -independent singlet potential V RF S (r) with N 3 LL (leading log) accuracy; see Ref. [5] for details.V RF

S
does not have renormalon uncertainties or factorization scale dependence.The factorization scale dependent part of V S is absorbed into the nonperturbative ME.By this, µ f dependence of the nonperturbative ME vanishes as well [13].In this way, one can resolve the mixing of the Wilson coefficient V S with the nonperturbative term δE US , and obtains where each term is free of renormalons and µ f .In our analysis, we regard δE RF US as the non-local gluon condensate 2 of order Λ 3 QCD r 2 .
1 More accurately, dominant renormalons which arise from the q-integral are removed.Renormalons contained in α V (q) are subdominant and have not been well studied, which we neglect in this analysis. 2Proper treatment of δE US depends on distance region: it is a perturbative contribution when the ultrasoft scale ∆V (r) = C A α s /(2r) (with C A = 3) satisfies ∆V ≫ Λ QCD , while it is a nonperturbative condensate when ∆V Λ QCD .Although ∆V ≫ Λ QCD is satisfied at very short distances, ∆V and Λ QCD have similar sizes at Λ QCD r 0.2.Since our fitting range extends to relatively long distances Λ QCD r < 0.6-0.8,we regard δE US as a nonperturbative contribution.(Validity of this treatment is shown in Fig. 2.) The OPE prediction is compared to the potential V latt calculated nonperturbatively in 3-flavor lattice QCD in the isospin limit [14].We employ the Symanzik gauge [15] and Möbius domain-wall quark actions [14,16], with which the leading discretization effect is O(a 2 ).The lattice simulations are carried out at three lattice cutoffs, determined as a −1 = 2.453(4), 3.610(9) and 4.496 (9) GeV from the Wilson-flow scale [17].The lattice sizes at these cutoffs are 323 ×64, 48 3 ×96 and 64 3 ×128, respectively, with the physical size roughly kept fixed.At each a −1 , we take a single combination of the light and strange quark masses (m latt ud , m latt s ), which roughly correspond to (M π , M K ) ∼ (300 MeV, 520 MeV).We extract V latt from Wilson loops with the spatial Wilson lines parallel to the spatial directions (1, 0, 0) and (1, 1, 0), denoted as directions 1 and 2, respectively.
We determine α s with the following strategy, by performing two analyses with different methods.The first analysis [Analysis (I)] consists of two steps: extracting a continuum limit of the lattice result and determination of α s by comparing the OPE prediction with the continuum limit.We proceed while checking (i) if the lattice data can be smoothly extrapolated to the continuum limit, and (ii) if V RF S (r) can explain the lattice result up to nonperturbative effects of O(r 2 ).After confirming these features, we perform a global fit to determine α s in the second analysis [Analysis (II)], without separating continuum extrapolation of the lattice data and extraction of α s .Analysis (II) is a first-principle analysis, which avoids introducing a model interpolating function, required in the first analysis for continuum extrapolation.Our final result will be adopted from Analysis (II), whose errors are well controlled and are smaller than that of Analysis (I).Analysis (I) makes up for a shortcoming of Analysis (II) that the output follows from the inputs without revealing detailed profiles at intermediate steps.Throughout our analyses, correlations among the lattice data are taken into account by the covariant matrices and the jackknife method.(See Ref. [18] for the details of the analyses.)Analysis (I) We first extract the continuum limit of the lattice data, specifically X latt (r) ≡ r 1 [V latt (r) − V latt (r 1 )], where r 1 is the scale defined by r 2 1 dV dr (r 1 ) = 1.To construct a sequence of X latt (r; a) at the same r but different a's, we first interpolate the lattice data at each a using the fitting form where d = 1, 2 and i = 1, 2, 3 specify the direction and lattice spacing, respectively.The first three terms are the Cornell potential.The other terms are included to take into account lattice artifacts.The 1/r 3 term accounts for the O(a 2 ) discretization effect whose mass dimension is one.The last term similarly accounts for the finite volume effect.Due to the lack of rotational symmetry on the lattice, the coefficients in Eq. ( 5) can depend on d.Hence, we interpolate the data separately for each (d, i).From the fit (5), we calculate X latt (r; a i ) at each a i and at reference values of r (in physical units) where the coarsest lattice has the original data.
We then extrapolate X latt (r; a) to a → 0 by linear fits in a 2 .The continuum result is shown in Fig. 1, where only the points which satisfy χ 2 /d.o.f.< 2 in extrapolation are adopted.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.
Before determining α s from the obtained lattice result, as a consistency check, we confirm that V RF S (r) can explain the lattice result up to nonperturbative corrections of O(r 2 ).Since the lattice result and V RF S are obtained in different units (r 1 and Λ MeV [1] and use the central value of r 1 = 0.311(2) fm [19] to convert the lattice result to that in Λ MS units.We see in Fig. 2 that the difference between the lattice data and V RF S can be fitted well by a constant plus an -term at Λ MS r 0.8.This is a first numerical observation which suggests correctness of OPE of pNRQCD in this distance range.
We also examine consistency of the lattice data at a = 0 with other predictions: V S in Ref. [8] with N 3 LL accuracy, and the fixed-order prediction at N 3 LO.These contain the all the points are smoothly extrapolated to a → 0 with χ 2 /d.o.f.< 2. In contrast, if we include the data at r = a, X latt does not obey a linear behavior in a 2 , since the interpolating function is seriously distorted by the data at r ∼ a. O(Λ 3 QCD r 2 ) renormalon, and the consistency is confirmed in a limited range Λ MS r 0.55 or strongly depends on the choice of the renormalization scale µ.In contrast, in our formulation without the O(Λ 3 QCD r 2 ) renormalon, the validity range is enlarged to Λ MS r 0.8 and stable against scheme choice for RG improvement, which shows an advantage of our framework.See [18] for details.
Our α s determination reduces to the problem to find an appropriate x = Λ MS r 1 where the lattice result agrees with the OPE prediction.We use OPE including up to the r 2 -term as the theoretical prediction, where A 0 is an r-independent constant and A 2 specifies the size of the leading nonperturbative effect (they are treated as fitting parameters).We obtain x = 0.496 ± 0.024(stat) from the data at Λ PDG MS r < 0.8, adopting the range in which OPE is reliable.The obtained 3-flavor Λ MS gives the 5-flavor coupling α s (M Z ) = 0.1166 +0.0010 −0.0011 (stat) through 4-loop RG evolution with the charm and bottom quark threshold corrections [20].The size of the nonperturbative effect is estimated as A 2 /Λ 3 MS = 0.04 ± 0.22(stat).To evaluate systematic errors, we perform the following re-analyses.(I-a) Finite a effect: An analysis including shorter distance points r > a is performed.(I-b) Interpolating function: The fitting function (5) does not contain log r corrections (dictated by RG) in the Coulomb part at small r.We use another interpolating function consistent with the 1-loop RG at small r. (I-c) Subtraction point: We extract r 1 [V latt (r) − V latt (0.8r 1 )], where we change the subtraction point of the potential.(I-d) Higher-order corrections to V RF S : We replace V RF S at N 3 LL by V RF S ± δV RF S , where δV RF S is the difference between the N 3 LL and N 2 LL results.(I-e) Matching range: To examine stability of OPE truncated at O(r 2 ), the continuum result satisfying Λ PDG MS r < 0.9 or 0.7 is used instead of 0.8.(I-f) Ultrasoft (US) contribution: α V (q) at 3-loop contains an IR divergence, which is canceled by an extra contribution from the US scale [21,7].In the main analysis we use the LO perturbative result for the US contribution, whereas in the error analysis we regard the US contribution as dominantly nonperturbative and introduce a cutoff at µ US = 3Λ MS or 4Λ MS .(I-g) Error of r 1 : The scale r 1 = 0.311(2) fm is varied within its error.We summarize the systematic errors in α s determination in table 1.
Analysis (II) Extrapolation to continuum limit and matching with the OPE prediction are performed by a global fit in one step.It is based on an idea that the OPE prediction should coincide with the lattice result at small r besides discretization effects.The lattice data, after correcting for discretization effects, are given by The second term is included to remove finite-a effects at tree level [κ = O(α s )], where 1 r is the LO perturbative result in lattice theory with finite a and L; the third term is included for removing the remaining O(α 2 s a 2 ) effect.We determine Λ MS in GeV units by comparing the above corrected lattice data to the OPE prediction V RF S (r) + A 2 r 2 , where each dataset is converted to GeV units using the estimated a −1 i [GeV].In this global fit, there are 16 parameters in total: Λ MS , six A 0 's, A 2 , six κ's and two f 's.Since we have more effective data than the first analysis, we shift the fitting range to shorter distances.It serves to reduce the higher order uncertainty, which is the dominant error in our analysis.We use the lattice data at Λ PDG MS r < 0.6.We obtain Λ MS = 0.334 ± 0.010(stat) GeV, giving α s (M Z ) = 0.1179 ± 0.0007(stat). 4For A 2 , we have A 2 = −0.0091± 0.0054(stat) GeV 3 .
We consider the following systematic errors.Since our final result is obtained from Analysis (II), we consider systematics errors more in detail than in Analysis (I).(IIa) Finite a effect: We drop the data at r < 2a, while the data at r a are used in the main analysis. 5(II-b) Mass corrections: The input (u, d, s) masses in each lattice simulation differ from the physical point.Since the nonperturbative correction due to these mass differences is unknown, we treat it as a systematic error.The lattice data at the physical point are estimated using perturbation theory as V latt (r; m latt ) → V latt (r; m latt ) + [V pt (r; m) − V pt (r; m latt )], where V pt is a finite mass effect evaluated in perturbative QCD at N 2 LO [22] with the MS masses m.We also substitute a constituent quark mass of 300 MeV for m to estimate the correction.Furthermore, since in the main analysis we use V RF S in the massless approximation, finite mass corrections are added.(II-c) Higher-order corrections to V RF S : An analysis parallel to the first one is performed.(II-d) Matching range: The upper limit of r is varied as Λ PDG MS r < 0.8 or 0.5.(II-e) US contributions: An analysis parallel to the first one is performed.(II-f) Scheme dependence: The µ findependent part of V S varies by a choice of scheme.A different scheme practically causes an O(r 3 ) difference in the OPE prediction (6).We add an r 3 -term in the fit to remove the scheme dependence and see how α s varies.(II-g) Lattice spacing: The lattice spacing is shifted by its uncertainty.We also take into account the error of the Wilson-flow scale.(See Ref. [18] for details.)We summarize the systematic errors in table 1.
We present the results of Analysis (I) and (II) in Fig. 3, where one can see that they are mutually consistent.Analysis (II) is superior to Analysis (I) in the sense that it is a first-principle analysis and that our dominant error, higher order uncertainty, is reduced thanks to the use of shorter distance range.Hence, we adopt the result of Analysis (II) as our final result.
In this Letter we determined α s from the QCD potential by comparing the lattice result and OPE prediction after subtracting renormalons from the leading Wilson coefficient.

Figure 1 :
Figure 1: Lattice result for the QCD potential after taking the continuum limit in Analysis (I).

Figure 2 :
Figure 2: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.

Table 1 :
Systematic errors in α s (M Z ) (in units of 10 −4 ) estimated from the variations of α s (M Z ).