Determination of $\Lambda_{\overline{\textrm{MS}}}^{(n_f=2)}$ and analytic parameterization of the static quark-antiquark potential

While lattice QCD allows for reliable results at small momentum transfers (large quark separations), perturbative QCD is restricted to large momentum transfers (small quark separations). The latter is determined up to a reference momentum scale $\Lambda$, which is to be provided from outside, e.g. from experiment or lattice QCD simulations. In this article, we extract $\Lambda_{\overline{\textrm{MS}}}$ for QCD with $n_f=2$ dynamical quark flavors by matching the perturbative static quark-antiquark potential in momentum space to lattice results in the intermediate momentum regime, where both approaches are expected to be applicable. In a second step, we combine the lattice and the perturbative results to provide a complete analytic parameterization of the static quark-antiquark potential in position space up to the string breaking scale. As an exemplary phenomenological application of our all-distances potential we compute the bottomonium spectrum in the static limit.


Introduction
This article finalizes our attempts started in [1], and subsequently improved in [2], to determine Λ MS by matching the static quark-antiquark potential obtained form perturbation theory to lattice data. While this endeavor seems almost trivial on first sight, it actually is not and requires to deal with several problems and tricky issues. To keep our series of articles selfcontained and to allow for a fair comparison of the different approaches and their results, also in the present article we stick to quantum chromodynamics (QCD) with n f = 2 dynamical quark flavors. However, the described approach can be readily adopted to other setups.
Perturbation theory is based on series expansions in the strong coupling α s , which thus is required to be small. The physical coupling α s ≡ α s (µ) generically depends on a momentum scale µ, which can be considered as a measure of the typical momentum transfer in a given process. Due to asymptotic freedom of QCD, α s (µ) 1 for large values of µ, while α s (µ) 1 for small values of µ. In momentum space we have µ ∼ p, while in position space µ ∼ 1/r. Correspondingly, perturbative calculations of the static potential in QCD are limited to small quark-antiquark separations r or large relative momentum transfers p, respectively. They are conventionally carried out in momentum space, where the static potential is presently known up to O(α 4 s ). Lattice simulations are tailored to the manifestly non-perturbative regime of QCD. They are naturally performed in position space and allow for controlled insights into the static potential from a minimum distance of a few times the lattice spacing a. Hence, the basic idea to determine Λ MS from the quark-antiquark potential amounts to fitting the perturbative static potential to lattice data in the intermediate regime of separations (momentum transfers) where both approaches are expected to allow for trustworthy insights, using Λ MS as fitting parameter.
In our initial article [1], we pursued this strategy in position space. To this end we transformed the perturbative static potential from momentum space to position space via an ordinary Fourier transform. As this Fourier transform also receives contributions from momenta where the perturbative expression is no longer trustworthy, the resulting perturbative potential in position space suffers from uncontrolled contributions, significantly worsening its convergence behavior as compared to the original momentum space potential. It can, however, be shown that the introduction of an additional momentum scale can remove these pathologies [33][34][35][36][37][38][39][40], thereby facilitating a reliable extraction of Λ MS . Proceeding along these lines, in [1] we obtained Λ (n f =2) MS = 315 (30) MeV.
Subsequently, in [2] we argued that this analysis can be performed more reliably in momentum space. To obtain the lattice potential in momentum space, we employed a discrete Fourier transform in three dimensions. In order to increase the number of available data points for the discrete Fourier transform, governing the resolution in momentum space, we used a fitting function to extrapolate the long distance behavior inferred from the lattice simulations up to distances of several hundreds of the lattice spacing a. The values of the extrapolating function were then stored for the sites of our extended three dimensional lattice serving as supplementary input for the discrete Fourier transform. Pursuing this strategy, in [2] we found Λ (n f =2) MS = 331 (21) MeV. Note, that the results of both analyses [1] and [2] are compatible with each other, the latter exhibiting a smaller error.
In the present article, we further improve and streamline the procedure to extract Λ MS in momentum space. The main difference to [2] is that we do not perform a discrete Fourier transform of lattice data from position to momentum space. Instead, we immediately parameterize the discrete lattice data points for the static potential in position space by a continuous function, thereby providing us with a "continuous lattice potential". More specifically, we choose the parameterizing function such that its Fourier transform to momentum space can be performed analytically, implying that the construction of the lattice potential in momentum space becomes essentially trivial. Another improvement to our previous articles [1] and [2] is that in the meantime the knowledge of the QCD β-function has been improved by an additional order in α s [41]. This results in a more precise relation between the strong coupling α s (µ) and the dimensionless ratio µ/Λ MS , prospectively further diminishing the error of our extracted value for Λ MS .
However, the present article does not only aim at an accurate and efficient determination of Λ MS . The extracted value of Λ MS is immediately used to construct a complete analytic parameterization of the static quark-antiquark potential in position space up to the string breaking scale. This potential encodes both perturbative and manifestly non-perturbative information and has various phenomenological applications. For instance, it can be employed to study the spectrum of heavy quarkonia. Here, we exemplarily adopt it to bottomonium in the static limit.
More specifically, our article is organized as follows.
Section 2 is devoted to the lattice computation of the static quark-antiquark potential V lat for QCD with n f = 2 dynamical quark flavors in position space. Here, our main goal is to parameterize the discrete data points for the potential obtained from Wilson loop averages in Sec. 2.2 by a continuous function. In Sec. 2.3, we confirm that a simple three parameter fit of the Cornell form, V lat = V 0 − α/r + σr, with offset V 0 and parameters α and σ, already accurately describes the lattice data points. We do not only extract the values of α and σ and their errors, but also account for their correlations. Section 3 focuses on the perturbative static potential. After briefly reviewing the known contributions to the static potential in momentum spaceṼ pert in Sec. 3.1, we discuss its position space analogue V pert in Sec. 3.2. Section 3.3 details on the relation of the perturbative strong coupling α s (µ) at a given momentum scale µ to Λ MS . In general, the constituting equation relating the dimensionless ratio µ/Λ MS to α s (µ) cannot be solved analytically for α s (µ). However, an analytical expression for α s (µ) can be extracted in the limit of µ/Λ MS 1.
In Sec. 4 we determine Λ MS for QCD with n f = 2 dynamical quark flavors. To this end, we fitṼ pert to the Fourier transform of the continuous lattice potential V lat in the intermediate momentum regime where both lattice simulations and perturbation theory are expected to allow for trustworthy results, treating Λ MS as fitting parameter.
Using the extracted value of Λ MS , in Sec. 5 we construct a complete analytic parameterization of the static potential in position space V , interpolating between both the perturbative and the manifestly non-perturbative regime. To this end, V pert is smoothly connected to the continuous lattice potential V lat at intermediate quark-antiquark separations. In order to allow for an analytic representation of V pert , we employ the analytic expression for α s (µ) in the limit of µ/Λ MS 1.
In Sec. 6 we study the bottomonium spectrum in the static limit. This analysis serves as an exemplary application of the all-distances potential V constructed in the preceding section.
Finally, we end with conclusions and a brief outlook in Sec. 7.
2 Lattice computation of the static potential

Lattice setup
We use the same n f = 2 gauge link configurations as in our previous articles concerned with the determination of Λ (n f =2) MS [1,2]. These gauge link configurations were generated by the European Twisted Mass Collaboration (ETMC) [42][43][44]. The gluon action is the tree-level Symanzik improved gauge action [45], with b 0 = 1−8b 1 and b 1 = −1/12, and the quark action is the Wilson twisted mass action [46][47][48][49], with ∇ µ and ∇ * µ are the gauge covariant forward and backward derivatives, m 0 and µ q are the bare untwisted and twisted quark masses, τ 3 is the third Pauli matrix acting in flavor space, and χ = (χ (u) , χ (d) ) represents the quark fields in the so-called twisted basis.
The twist angle ω is given by ω = arctan(µ R /m R ), where µ R and m R denote the renormalized twisted and untwisted quark masses. For the ensembles of gauge link configurations considered in the present study ω has been tuned to π/2 by adjusting m 0 appropriately. This ensures automatic O(a) improvement for many observables including the static potential (cf. [43] for details).
The considered gauge link configurations cover several different values of the lattice spacing (cf. Table 1, which also provides the corresponding pion masses m PS , spacetime volumes (L/a) 3 × T /a and numbers of gauge link configurations used for the computations of the static potential). The lattice spacing in physical units has been set via the pion mass and the pion decay constant using chiral perturbation theory. The resulting value for the hadronic scale 1 r 0 is r 0 = 0.420 (14) fm (cf. also Sec. 5 of [43] and Table 8 of [44]). For further details on the generation of these gauge field configurations as well as on the computation and the analysis of standard quantities, such as lattice spacing and pion mass, we refer to [43,44].

Extracting the lattice static potential from Wilson loop averages
We extract the static potential in position space V lat ( r) from the exponential decay of Wilson loop averages W ( r, t) with respect to their temporal extent t, while keeping their spatial extent 1 The hadronic scale r0 is defined via r 2 0 F (r0) = 1.65, with F (r) = dV (r)/dr [50].  r fixed [51]. To this end we first compute In a second step the t-independent quantity V lat ( r) is obtained by performing an uncorrelated χ 2 minimizing fit to V (effective) lat ( r, t) in a suitable t range. This range is chosen such that the contributions of excited states are strongly suppressed, while statistical errors are still small.
We perform two independent computations on each of the ensembles listed in Table 1.
• no-HYP computation: Temporal links remain unchanged, i.e. are not smeared. The resulting static potential has small discretization errors, in particular at small quark-antiquark separations r = | r|, but large statistical errors at large separations. To obtain a fine resolution at small r, we consider both on-and off-axis Wilson loops (for a detailed explanation regarding the construction of off-axis Wilson loops, cf. [1]). Spatial links are APE smeared, to improve the ground state overlap and, hence, to be able to extract the static potential more precisely (N APE = 20, α APE = 0.5; cf. [52] for details). To further reduce discretization errors we employ a tree-level improvement technique, which is explained in [1].
• HYP computation: Temporal links are HYP2 smeared, which corresponds to using the HYP2 static quark action [53][54][55]. The resulting static potential has large discretization errors at small quarkantiquark separations, but the reduced self energy of the static quarks leads to significantly smaller statistical errors at large separations. We consider only on-axis Wilson loops. Spatial links are again APE smeared (N APE = 60, α APE = 0.5).
While the no-HYP results have already been used in our previous determinations of Λ (n f =2) MS [1,2], the HYP results have been generated for this study.

Parameterization of the discrete lattice data by a continuous function
In contrast to our previous determination of Λ (n f =2) MS in momentum space [2], we parameterize the discrete lattice data for the static potential by a continuous function before transforming to momentum space. This has several advantages. For example rotational symmetry is restored already at an early stage, thereby avoiding technical problems like performing a cylinder cut. Moreover, this approach is technically less complicated and the uncertainty of the final result for Λ MS is somewhat smaller; see below. Finally, the continuous function used to parameterize the lattice potential forms an important constituent of the analytic all-distances potential constructed in Sec. 5.
For quark-antiquark separations r min ≤ r ≤ r max , with r min > ∼ 0.13 fm and r max ≤ 0.79 fm (the maximum range of separations, where results are available) the lattice potential computed on all four ensembles can be parameterized consistently by the Cornell potential, Here, V 0 is a physically irrelevant constant shift of the potential and σ is the string tension. While α = π/12 ≈ 0.26 for large r in the bosonic string picture [56,57], lattice simulations with n f = 2 quark flavors have extracted α ≈ +0.3 . . . + 0.5 [58], which is in agreement with our results. We have explicitly checked, whether accounting for additional terms ∼ ln m (r)/r, with m ∈ {1, 2, 3}, which can be Fourier-transformed analytically, improves the parameterization of the static potential (5). We find that at our current level of statistical precision such terms are not needed. Their prefactors are zero within statistical errors.  Table 2: Parameterization of the lattice static potential via V lat (r) = V 0 − α/r + σr. To allow for a straightforward comparison, the results for σ at fixed β (upper eight lines) have been converted from lattice units to 1/fm 2 without accounting for the lattice spacing errors.
The parametrization (5) does not account for string breaking, which is happening at quarkantiquark separations r ≈ r sb , where V (r sb ) = 2m B , with m B denoting the mass of the lightest heavy-light meson (quantum numbers J P = 0 − , 1 − ; cf. e.g. [52,59]). The string breaking distance has been determined using lattice QCD in [60], yielding r sb = 1.13(10)(10) fm. In this section we do not use any lattice data for r > r sb , and hence do not consider string breaking; cf. Sec. 6 for a study of the effect of string breaking on the bottomonium spectrum.
In the present study, we determine the parameters V 0 , α and σ by performing uncorrelated χ 2 minimizing fits of Eq. (5) to the discrete lattice QCD data of the static potential for separations r min ≤ r ≤ r max ; cf. Table 2 for the explicit values of r min and r max . The minimum distance r min is needed to exclude data points with sizable lattice discretization errors, and r max is required to exclude unwanted artifacts of the spatial periodicity of the lattice. In detail we proceed as follows: (i) For each of the four ensembles of β and both X ∈ {no-HYP, HYP} computations, we determine the averages α β,X andσ β,X , whereσ = σa 2 is the string tension in units of the lattice spacing. The errors ∆α β,X and ∆σ β,X are computed via the jackknife method. Moreover, we determine the correlation corr(α β,X ,σ β,X ). For two generic quantities α and σ, the latter is defined as The results of these calculations are collected in Table 2. Note that results obtained for the same β but different X ∈ {no-HYP, HYP} may disagree within statistical errors, because of different discretization errors.
(ii) For both no-HYP and HYP results we perform continuum extrapolations of α and σ to linear order in a 2 , which is the leading order of discretization errors in Wilson twisted mass lattice QCD at maximal twist (cf. Sec. 2.1). Correlations of α and σ are properly taken into account by using jackknife samples (α, σ) from step (i). Moreover, we account for the lattice spacing errors listed in Table 1 when converting a dimensionlessσ to a dimensionful σ. Since the lattice spacing errors constitute the dominant source of uncertainty for σ, they also reduce the correlation corr(α, σ) significantly; cf. the upper eight lines to the lower three lines of Table 2. Furthermore, notice that pion masses and spacetime volumes for different lattice spacings are similar, but not identical (cf. Table 1). We do not consider this as problematic, since the dependence of the potential on the pion mass and the spacetime volume is negligible within statistical errors [1]. This is also supported by the small reduced χ 2 of the continuum extrapolations of both our no-HYP and HYP computations. These extrapolations are shown in Figure 1, and the corresponding results for α X , ∆α X , σ X , ∆σ X and corr(α X , σ X ) are collected in the 9th and 10th line of Table 2. As expected, the continuum extrapolated no-HYP and HYP results for both α and σ are in agreement within statistical errors.
(iii) We combine the continuum extrapolated no-HYP and HYP results from step (ii) by performing constant fits. Correlations between α and σ are properly taken into account by using jackknife samples (α, σ) from step (ii). The fits are shown in Figure 2, and the results for α, ∆α, σ, ∆σ and corr(α, σ) are collected in Table 2. Figure 2 shows that the no-HYP results constrain α more precisely than the HYP results. Still, the HYP results reduce the final uncertainty of α by around 15%. The uncertainties of both the no-HYP and HYP results for σ are dominated by the lattice spacing error. Consequently, the combination of both results reduces the final uncertainty of σ only marginally; cf. also Figure 2. In any case, we consider it as reassuring to have two independently computed datasets based on different values of β, yielding perfectly compatible results.
The continuum extrapolation of the analytic parameterization (5) of the lattice static potential combining both no-HYP and HYP computations (cf. last line of Table 2) is shown in Figure 3. To cross-check our results with previous ETMC analyses, we determine the continuum extrapolated r 0 from our results for α and σ via yielding r 0 = 0.420 (15). This value is in perfect agreement with r 0 = 0.420(14) (extrapolated to the continuum and the chiral limit) from [44].
We note that there is an anti-correlation between α and σ, cf. corr(α, σ) = −0.17 for our final continuum result. This can be explained by the fact that both increasing α and σ results in a larger slope of V lat (r), implying that these two parameters have a similar effect on the shape of V lat (r). Hence, for precise statistical analyses based on the static potential, e.g. the determination of Λ MS or the computation of the bottomonium spectrum, as done in Sec. 4 and Sec. 6, this anti-correlation should be taken into account. In Sec. 5 we discuss in detailed how to include this anti-correlation in the computation of any observable, which makes use of the analytic parameterization (5) of the static quark-antiquark potential.  Table 2. The dark green lines correspond to our final values for α and σ, and the green bands depict our final errors ∆α and ∆σ, respectively.   Table 2), while the gray error band has been generated from the jackknife samples of step (iii).
3 The perturbative static potential

The perturbative static potential in momentum space
In perturbation theory the static quark-antiquark potential V is conventionally determined in momentum space. For gauge group SU (3), it can be expressed as with p = | p| > Λ QCD . The latter condition implies an explicit restriction to the perturbative momentum regime of QCD. The dimensionless quantity α V is a function of both the strong coupling α s (µ) evaluated at the renormalization scale µ (cf. also Sec. 3.3 below) and ln(µ 2 /p 2 ). It is explicitly known up to O(α 4 s ). The static potential is a renormalization group invariant, implying invariance of V under a change of µ. In perturbation theory this means that when evaluating a result known up to O(αk) for two different choices of µ, the differences among these results are relegated to O(αk +1 ), such that for small enough α s and large enoughk the specific choice of µ eventually becomes irrelevant.
Note, that if µ is chosen as µ = cp, where c denotes a proportionality constant, the logarithm in the argument of α V becomes independent of p. For c = 1 we have ln(µ 2 /p 2 )| µ=cp = 0, such that Eq. (8) can be written in a particularly compact form and α V becomes a function of α s (p) only. Adopting the latter choice, the known terms of the static potential can be represented as In the remainder of this article we will refer to Eq. (9), utilizing the identification µ = p, as the static potential in momentum space. A truncation of Eq. (9) accounting for terms up to O(α 1+n s ) is subsequently referred to as (next-to-) n leading-order or N n LO, respectively. The coefficients a 1 [61, 62] a 2 [63][64][65] and a 3ln [66,67] are known analytically, while some contributions to a 3 are only known numerically [68][69][70][71][72]. For gauge group SU(3), n f = 2 dynamical massless quark flavors and in the MS scheme [73,74], they read , a 3ln = 144π 2 , a 3 = 8783.16 (38) .
The running of α s (µ) with µ is governed by the QCD β-function, Its perturbative expansion in α s is presently known with the following accuracy, Equations (11) and (12) imply that α s (p) can be expressed in terms of α s (µ) and β 0 ln(µ 2 /p 2 ). Formally expanding α s (p) in powers of α s (µ) and solving these equations order by order in α s (µ) one obtains With the help of this identity the couplings in Eq. (9) can be promoted to any other renormalization scale. Upon insertion into Eq. (9), we recover the structure of Eq. (8), with α V known explicitly up to O(α 4 s ).

The perturbative static potential in position space
When choosing µ as independent of p, Eq. (8) can be straightforwardly transformed to position space by means of a standard Fourier transform in three dimensions, However, note that the Fourier integral (15) naturally includes momenta p Λ QCD for which perturbation theory is no longer trustworthy, potentially inducing uncontrolled contributions.
Contributions of this kind are already present in the perturbative potentialṼ pert (p). In standard perturbation theory loop diagrams come along with integrations (2π) 4 of the loop fourmomentum q over the full momentum regime, i.e., also receive contributions from outside the perturbative momentum regime. The leading uncontrolled contribution toṼ pert (p) is quadratic in Λ QCD and scales as ∼ − 4παs , translating to a term ∼ − αs r (rΛ QCD ) 2 in V pert (r). Obviously, the condition for momentum transfers p to be described reliably in perturbation theory, Λ QCD p 1, corresponds to the restriction rΛ QCD 1 in position space.
However, the latter problem can be cured by manifestly restricting the Fourier integral to the perturbative momentum regime p ≥ µ f > Λ QCD , where µ f denotes a momentum cutoff still in the perturbative regime; cf. [36,76,77]. The position space potential as defined by a restricted Fourier transform, does not suffer from enhanced pathological terms in comparison to V (p). In this article, we will adopt this definition of the perturbative potential in position space. Equation (16) can be evaluated analytically [77].
To allow for a compact representation of the explicit expressions for V pert (r) and δV pert (r, µ f ), it is convenient to introduce the following polynomials of degree k [78], with dimensionless expansion coefficients ρ km . For 1 ≤ k ≤ 3 we have ρ k0 = a k , and ρ kk = β k 0 . Moreover, we use the shorthand notation P (n) k (L) = ∂ n ∂L n P k (L). Setting the renormalization scale to µ = 1/r, we obtain and In the latter expression we limited ourselves to the leading term in an expansion in powers of rµ f [36]. This is completely sufficient as the terms explicitly accounted for in Eq. (20) are exactly those canceling the pathological contribution ∼ − αs r (rΛ QCD ) induced by the Fourier integral (15) in Eq. (16). Uncontrolled higher-order terms cannot be fully eliminated along these lines anyway.

The perturbative coupling α s (µ) and its relation to Λ MS
Up to now, we did not specify how the strong coupling α s (µ) is promoted to an explicit numerical value. This identification involves the definition of a reference scale, which -in the MS scheme -is denoted by Λ MS . More specifically, the scale Λ MS is introduced in terms of specific initial conditions in the integration of Eq. (11) [79] (cf. also [2]), This scale cannot be determined within perturbation theory, but has to be provided as an external input parameter. In this study, we aim at determining Λ (I) either plugging the perturbative β-function (12) at the best known accuracy into Eq. (21) and performing the integration over α s numerically, (II) or adopting a Taylor expansion of the integrand in Eq. (21) and performing the integral analytically. To this end only those terms whose coefficients are known explicitly are kept, i.e., These two choices may also serve as a consistency criterion as in the manifestly perturbative regime both options should, of course, yield compatible results.
Terms beyond O(1/ 5 ) also include higher, to date unknown coefficients b i≥5 of the QCD βfunction.
Even though the options (I) and (II) allow for reliable results for α s (µ) in a wider range of µ, Eq. (23) is specifically useful when aiming at an analytic expression of the static potential in the limit of µ Λ MS . Hence, in our determination of the value of Λ MS by matching the perturbative static potential to lattice results we will exclusively resort to the options (I) and (II). Contrarily, for the analytic expression of the position space potential we will adopt Eq. (23).

Determination of Λ MS for n f = 2 massless quark flavors
In this section we determine Λ MS by matching the perturbative static potentialṼ pert in momentum space, Eq. (9), to the Fourier transform of our analytic parameterization (5) of the lattice static potential V lat . Recall that the lattice static potential in momentum space allows for reliable insights below a certain momentum, while the perturbative results is applicable above a certain momentum. In turn, the matching procedure has to be done in the momentum regime where the validity of both results overlaps.

Matching perturbative and lattice QCD results for the static potential
The analytic parameterization of the lattice potential V lat (r) in Eq. (5) can be straightforwardly transformed to momentum space, where the physically irrelevant constant V 0 has been set to zero. In this section, we exclusively adopt the "continuum/combined" values for the parameters α and σ listed in the last line of Table 2.
For the perturbative static potentialṼ pert (p), we adopt the expression given in Eq. (9), with the implicit equation (21) inverted by one of the options (I) or (II) for the strong coupling α s (p). In turn,Ṽ pert (p) is determined up to a single parameter, namely Λ To matchṼ pert (p) andṼ lat (p), we minimize their squared difference, with respect to Λ (n f =2) MS in a given momentum interval p min ≤ p ≤ p max . The lower and upper momenta, p min and p max , respectively, are chosen such that bothṼ pert andṼ lat exhibit small systematic errors, i.e. errors due to the lattice discretization and the truncation of the perturbative series. The minimization is done numerically using standard integration and root finding techniques. The small difference in the results obtained by options (I) and (II) is included in the systematic error for the final result (cf. the discussion below). For an exemplarily plot, cf. In our previous determination of Λ (n f =2) MS in momentum space [2] we additionally had to account for a constant offsetṼ 0 in the matching ofṼ pert (p) andṼ lat (p). Such a constant was needed in [2], because the matching was performed with discrete lattice data points for the static potential in momentum space obtained by means of a discrete Fourier transform from the corresponding lattice data in position space. Due to lattice discretization errors at small separations r, the discrete lattice data points in position space do not exhibit a singularity for r → 0, but form a negative peak, saturating at a finite value. The shape of the peak depends on the lattice spacing a and becomes more pronounced for smaller a. Such a peak in position space translates into an a-dependent constant in momentum space. Contrarily, in the present work we parametrize the lattice data points in position space with the function (5), intrinsically excluding lattice data points for small r exhibiting large errors, and thus do not need to account for such a constant shift. A further argument against the necessity of such a constant shift in the present study is the fact that bothṼ lat (p) andṼ pert (p) exhibit the same asymptotic behavior for p → ∞. For consistency, we checked this numerically and foundṼ 0 = 0 within statistical errors.

Variation of input parameters, final result and uncertainty of Λ MS
To allow for a fair comparison of the results, we determine the systematic error of Λ (n f =2) MS in the same way as in our previous articles [1,2]. More specifically, we perform the matching procedure outlined in Sec. 4.1 for 20 000 times, while varying the input parameters as follows: • 50% of the matching is done withṼ pert as defined in Eq. (9) at NNLO, 50% at NNNLO.
• For each matching, we randomly choose Moreover, for each matching we randomly pick one of the jackknife samples (α, σ) used in step (iii) of Sec. 2.3 to generate our final "continuum/combined" results for α and σ, listed in the last line of Table 2. This ensures that the statistical uncertainties of these parameters and their correlation are properly taken into account.
The impact of finite volume effects on Λ (n f =2) MS was studied in [1] and found to be negligible in comparison to other errors. Similarly, the effects of non-vanishing light quark masses on Λ to be stable and constant within tiny statistical errors of ≈ ±1 MeV. Hence, we do not expect the non-vanishing light quark masses in the lattice QCD computation to induce any significant deviations from the limit of massless dynamical quark flavors assumed in the derivation of the perturbative static potential in Sec. 3. For these reasons, in the present study we do not account for any potential errors arising from finite volumes and non-vanishing light quark masses on the lattice.
Performing the matching procedure 20 000 times, we obtain 20 000 samples for Λ = 331(21) MeV, obtained from a momentum space analysis of the static potential in [2], our new result is slightly smaller, but still compatible within errors. Moreover, the error of our new result (26) is reduced by roughly 25%, which is why we consider the procedure presented here as superior. The main differences to our previous determination of Λ • a more straightforward technique: we first parameterize the discrete lattice potential by a continuous function (5), which is then straightforwardly transformed to momentum space. For instance, a cylinder cut as employed for the discrete Fourier transform in Sec 2.3 of [2] is no longer required.
• meanwhile, with b 4 an additional term in the perturbative expansion of the QCD βfunction (12) has been determined [41] and is included in our analysis.
For completeness, also note the result of our initial determination of Λ MS from a position space analysis of the quark-antiquark static potential based on the same lattice QCD data [1], yielding Λ (n f =2) MS [1] = 315(30) MeV, which is also consistent, but exhibits a larger error. This is related to the observation that our momentum space determinations of Λ The mean value of this result, is 11 MeV smaller than that in Eq. (26), based on the combination of both NNLO and NNNLO input. At the same time, the error is reduced by roughly 25%, which is consistent with the findings of our previous momentum space analysis [2], yielding Λ

Complete analytic parameterization of the static potential
The analytic parameterization of the lattice potential V lat (r) derived in Eq. (5), and the perturbative static potential V pert (r) as defined in Eqs. (16)- (20), with the strong coupling given by Eq. (23) and µ = 1/r, can be combined to provide a complete analytic parameterization of the quark-antiquark static potential V (r) for n f = 2 valid up to the string breaking distance. This parameterization will be useful for various applications, such as heavy-quark phenomenology.

Construction of the analytic parameterization of V (r)
The analytic parameterization (5) of the lattice potential V lat (r), with parameters α and σ fixed to the "continuum/combined" results provided in the last line of Table 2, accurately describes the quark-antiquark potential for heavy-quark separations r > ∼ 3a (cf. Table 2). For our smallest lattice spacing, the latter condition corresponds to r > ∼ 0.12 fm. This expression is valid up to the string breaking distance r sb ≈ 1 fm. On the other hand, the perturbative static potential V pert (r) at NNNLO, Eqs. (16)- (20), with the strong coupling given by Eq. (23), µ = 1/r, fixed to the value provided in Eq. (26), is expected to exhibit small systematic errors for r < ∼ 0.12 fm [77].
To provide an analytic parameterization of V (r) for r < ∼ r sb , we use these two results and connect them in a smooth way, allowing for a constant offset V 0 between the perturbative and lattice static potentials. More specifically, we define where r 1 , r 2 ≈ 0.1 fm . . . 0.2 fm, fulfilling r 1 < r 2 , and is a second degree polynomial. The coefficients A, B, C and the constant shift V 0 are chosen such that V (r) is a continuous and smooth function, by solving the following system of four simple linear equations, with the primes denoting differentiations with respect to r.

Error analysis, when using the analytic parameterization of V (r)
In Secs. 2-4 possible uncertainties and errors associated with the lattice QCD computation and the perturbative calculation were discussed and quantified in detail. For example we did not only provide individual errors for α and σ, but also accounted for their correlation. In this section, we propose the following procedure, allowing for the proper inclusion and propagation of these uncertainties to observables whose determination is based on our complete analytic parameterization of the static potential (28) (an explicit example is presented in Sec. 6): • Let X denote the observable. For example X could be a specific difference of two bottomonium masses, e.g. X ≡ m η b (2S) − m η b (1S) (cf. Sec. 6).
• Repeat the calculation of X very often, namely N 100 times, by randomly sampling the parameters of V (r). The results are X j , with j = 1, . . . , N . To this end we choose: α and σ according to the following 2-dimensional Gaussian probability distribution, characterized by the parameters α, ∆α, σ, ∆σ and corr(α, σ) of the "continuum/combined" results listed in the last line of Table 2, according to a Gaussian probability distribution parameterized by our result in Eq. (26). For completeness, note that in principle, Λ (n f =2) MS , α and σ are also correlated. As this correlation is rather small, we neglect it in the error analysis.
• The mean value of the obtained results X j serves as our estimate for X. Its error is defined via the standard deviation. More specifically, • Check that both X and ∆X are essentially independent of N . If not, increase N .

The bottomonium spectrum in the Born-Oppenheimer approximation
In the following, we use the parameterization of the static potential (28) to compute the bottomonium spectrum. To this end, we employ the Born-Oppenheimer approximation [80], which consists of two steps. Its first step amounts to the computation of the static potential, assuming the light quarks and gluons as dynamical degrees of freedom and the b quark and its antiquark b as static. In the second step, this constraint is relaxed, and the Schrödinger equation for the relative coordinate of thebb pair is solved with the potential computed in the first step, assuming a finite b quark mass m b .
Since the static potential manifestly neglects 1/m b corrections, encoding e.g. spin effects of the heavy quarks, this approach certainly does not allow us to obtain very accurate results for the bottomonium spectrum. We rather intend at performing an exemplarily calculation, utilizing the parameterization (28) of the static potential and including a full error propagation along the lines of Sec. 5.2. Moreover, it can be considered as a preparatory step for a more refined computation, accounting for such 1/m b corrections, to be determined within potential non-relativistic QCD and lattice QCD (cf. e.g. [81][82][83][84][85]). Also note, that the Born-Oppenheimer approach without 1/m b corrections has recently been used for the study of heavy tetraquarks (cf. e.g. [86][87][88][89][90][91][92]), where no experimental data is available. A comparison of the theoretical predictions for the bottomonium masses with corresponding experimental results, might provide us with an estimate of the systematic error associated with this approach.
To this end, we solve the Schrödinger equation where r is the relative coordinate of the bb pair, V (r) is the parameterization (28) of the static potential and µ = m b /2 the reduced mass of the b quark. For its explicit value, we employ either m b = m b,MS = 4.18 GeV determined in the MS scheme [93], or m b = m b,qm = 4.977 GeV from quark models [94]. Since the potential is radially symmetric, Eq. (33) can be separated in a radial and an angular equation, with the latter being trivial to solve. The radial equation reads − 1 2µ where we made use of the ansatz ψ( r) = (u n,l (r)/r)Y l,m (ϑ, ϕ). The solutions u n,l (r) are labeled by the principal quantum number n and the azimuthal quantum number l, corresponding to the orbital angular momentum of thebb pair. This equation can be solved numerically using standard techniques. Here we use a 4th order Runge-Kutta shooting method combined with Newton's method for root finding.
The resulting energy eigenvalues E n,l can be related to bottomonium masses M n,l via where m η b (1S) = 9399 MeV is fixed by experimental input [93]. Clearly, experimental input is needed to calibrate the unknown constant shift between the energy eigenvalues E n,l and the bottomonium masses M n,l , which has its origin in the self energy of the static quarks and the lattice cutoff. In other words, using the Born-Oppenheimer approach one can only predict mass differences, but not absolute masses of bottomonium states. Moreover, the results are independent of the heavy quark spins, because the potential V (r) has been computed in the static limit. For example, for orbital angular momentum l = 0 there are degenerate J = 0 and J = 1 states, and for l = 1 there are degenerate J = l − 1, J = l and J = l + 1 states, with J denoting the total angular momentum of the bottomonium system.
In the following we present results for three different cases: where r sb = 1.13 fm is the string breaking distance determined with lattice QCD in [60].
Cases (A) and (B) allow us to infer, how strong bottomonium mass differences depend on the value of m b . 4 For cases (B) and (C) we keep the mass of the b quark fixed, and only alter the long-distance behavior of the static potential. This allows us to check whether string breaking effects have a sizable effect on the bottomonium spectrum. Of course, with the potential (36) and the approach adopted here, we can only determine bottomonium states below the bb threshold. States above are unstable resonances. In principle, the determination of such states is possible along the lines, but requires techniques from scattering theory (cf. e.g. [92]).
Our results for the mass differences ∆E n,l = E n,l − E 1,0 with l ∈ {0, 1} are collected in Table 3.
The statistical analysis has been performed with the method detailed in Sec. 5 Table 3: Mass differences ∆E n,l = E n,l − E 1,0 in units of GeV for the three cases (A), (B) and (C) discussed in the main text.
The resulting bottomonium masses are confronted with experimental data in Table 4 and Figure 6, where the common notation S for L = 0 and P for L = 1 orbital angular momentum is used. As discussed above, our computations do not account for the heavy quark spins, such that, e.g. the bottomonium states η b (1S) and Υ(1S) with quantum numbers J P = 0 − and J P = 1 − , respectively, are mass degenerate. Here, the upper label P refers to the parity of the state. In turn, the experimentally measured mass difference of the order of 50 MeV for these two states can serve as an estimate of the systematic error associated with our results.  For the low-lying states 1S, 2S, 1P and 2P our theoretical predictions are in good agreement with experiment within the expected systematic error of the order of 50 MeV. Higher states, in particular those above the bb threshold, should be treated with caution. Noteworthily, all masses below threshold, including those in its vicinity, as e.g. the 3S state, are essentially identical for cases (B) and (C), which implies that they are not affected by the flattening of the potential for heavy quark separations above the string breaking distance.
Let us also note, that the static bottomonium results of [76] are compatible with ours. Reference [76] does not account for string breaking effects, but models contributions inversely proportional to the heavy quark masses in the potential, encoding spin effects of the heavy quarks and hyperfine splitting. A possible future direction could be to compute 1/m b and 1/m 2 b corrections 9.5  Table 4. Our theoretical predictions for the three cases (A), (B) and (C) discussed in the main text are confronted with experimental data [93]. In our numerical computations, the mass of the 1S state is fixed to the state η b (1S) observed in experiment.

Conclusions and outlook
In this article, we have determined the parameter Λ MS for QCD with n f = 2 dynamical quark flavors by fitting the perturbative result for the static potential to lattice data in momentum space. Building on insights from our previous determinations [1,2], we have substantially improved and streamlined our strategy to extract the value of Λ MS , resulting in Λ (n f =2) MS = 302(16) MeV.
One of the main improvements devised in the present work is the use of an analytic parameterization of the discrete simulation data of the lattice static potential in position space. This renders complicated and time-consuming numerical techniques employed in our previous work [2] superfluous, such as the discrete Fourier transform combined with a cylinder cut, which possibly introduces large systematic errors. Besides, it immediately provides an analytical expression for the static quark-antiquark potential in the manifestly non-perturbative regime.
In a second step, we have used Λ MS as input parameter in the perturbative static potential.
Utilizing an approximate analytical expression for the strong coupling α s (µ) in terms of the dimensionless ratio µ/Λ MS valid for large values of µ, upon identification of µ = 1/r, the perturbative static potential in position space becomes an analytic function of the quark-antiquark separation r. This function accurately describes the small distance behavior of the static potential. Connecting it with the analytic parameterization of the lattice potential by means of an adequately chosen interpolating function, we have constructed a complete analytic parameterization of the static quark-antiquark potential in position space, valid up to the string breaking distance. If desired, the effect of string breaking can also be phenomenologically accounted for by letting the potential become constant beyond the string breaking distance or by using first principles lattice QCD input, e.g. from [60]. This all-distance potential encoding both perturbative and manifestly non-perturbative information has ample phenomenological applications.
As an immediate and exemplary phenomenological application, we have used this potential to determine the bottomonium spectrum in the static limit, based on the Born-Oppenheimer approximation. Fixing the lowest bound state with data provided by the Particle Data Group [93], all bound states below the bb threshold are in reasonable agreement with experiment. Note that small deviations are not surprising as in particular spin effects have been completely ignored in our current analysis.
Let us finally emphasize, that the strategy devised in the present work can be readily adopted to the determination of Λ MS , and for the construction of analytic all-distance potentials from other lattice configurations with, e.g., n f = 0, n f = 2 + 1 and n f = 2 + 1 + 1 dynamical quark flavors.