format: text

text:
# Perturbative Wilson loops with massive sea quarks on the lattice

## I Introduction

The calculation of Wilson loops in lattice perturbation theory is useful in a number of ways. An important application is the prediction of a strong coupling constant $`\alpha _{\overline{MS}}^{(5)}(m_Z)`$ from low energy hadronic phenomenology by means of non-perturbative lattice simulations El-Khadra:1992vn ; davi1 ; davi2 ; Spitz:1999tu ; Booth:2001qp ; davi3 . Perturbative calculations of Wilson loops are also employed in the context of mean field improvement programmes of the lattice action and operators Parisi:1980pe ; lepage .

In the limit of infinite Euclidean time separation, $`T\mathrm{}`$, Wilson loops give access to the perturbative lattice potential. Our results improve the understanding of the violations of rotational symmetry at short distances of this quantity, as well as of the short distance effects of including various flavours of sea quarks. Furthermore, in the limit of large distances, $`R\mathrm{}`$, the self-energy of static sources can be obtained from the potential, enabling the calculation of $`\overline{m_b}(\overline{m_b})`$ from non-perturbative simulations of heavy-light mesons in the static limit Martinelli:1999vt .

Perturbative computations of Wilson loops on the lattice now reach back as far as two decades. For this reason, much of the groundwork for the calculation covered in this article is well known, and has been reproduced by us. In addition to improving the precision of many previously existing results and presenting them in a more comprehensive way, we extend the calculations by incorporating the actions, sea quark masses and observables that are relevant for present day lattice simulations. We also perform the first calculation of non-planar Wilson loops. Because of the degree of replication we present a brief survey of the existing literature here.

In 1981 Müller and Rühl MullerRuhl  performed a highly rigorous calculation of the $`𝒪(g^4)`$ contribution to the Wilson loop for pure $`U(1)`$ and $`SU(2)`$ gauge theories and presented their result in terms of lattice integrals. In particular, they obtained a closed form for the $`T\mathrm{}`$ limit and hence the static potential. The plaquette in $`SU(N)`$ pure gauge theory was first calculated to this order by Di Giacomo and Rossi DiGiacomo:1981wt , again in 1981. In the same year Hattori and Kawai Hattori:1981ac  estimated $`𝒪(g^4)`$ contributions to larger Wilson loops, by modelling some of the relevant lattice contributions by their continuum counterparts. Rigorous results to this order were subsequently independently obtained by three groups of authors: Weisz, Wohlert and Wetzel Weisz:1984bn , Curci, Paffuti and Tripiccione Curci:1984wh  and Heller and Karsch Heller:1985hx . The first reference is the most general and also applies to Symanzik-improved pure Yang-Mills actions. However, numerical results on Wilson loops are only available for the Wilson action. The $`𝒪(g^4)`$ fermionic contribution to the plaquette was obtained by Hamber and Wu Hamber:1983ft  for Wilson quarks and for Kogut-Susskind (KS) quarks by Heller and Karsch Heller:1985eq . More recently, Panagopoulos et al. calculated the $`𝒪(g^6)`$ $`SU(N)`$ gluonic contribution Alles:1994dn  as well as the $`𝒪(g^4)`$ and $`𝒪(g^6)`$ corrections Alles:1998is  in presence of massive Wilson fermions. The only numerical calculation to-date of small Wilson loops for Symanzik improved gluonic actions was performed by Iso and Sakai Iso:1987dj .

More recently, the Heller and Karsch calculations Heller:1985hx  have been extended by Martinelli and Sachrajda Martinelli:1999vt  in a manner similar to that of Müller and Rühl MullerRuhl : the $`R,T\mathrm{}`$ limit of the perimeter term of Wilson loops has been taken to obtain the static self energy for $`SU(3)`$ gauge theory, with the addition of massless Wilson and Sheikholeslami-Wohlert (SW) sea quarks.

The structure of this paper is as follows. In Sec. II we describe the procedure used to calculate Wilson loops in perturbation theory and present the perturbative expansions of small Wilson loops, both planar and non-planar, to $`𝒪(g^4)`$ using a variety of actions. Although partially addressed by previous authors Weisz:1984bn ; Curci:1984wh ; Heller:1985hx ; Altevogt:1995cj  the static potential has not been well treated in literature, which we attempt to rectify in Sec. III, where we also obtain the static source self energy from the asymptotic $`R\mathrm{}`$ potential. This sets the stage for Sec. IV where these perturbative results are confronted with state-of-the-art lattice data from simulations with Wilson, SW and KS fermions: we determine the “$`\beta `$-shift” and the running coupling $`\alpha _{\overline{MS}}`$, before comparing perturbative and non-perturbative potentials at short distances. In Sec. V we conclude and summarize. Our article is augmented by two Appendices on the relations between different lattice and continuum renormalization schemes and on the perturbative $`\beta `$-shift.

## II Wilson loops in lattice perturbation theory

In this section we explain our method and notations, before displaying results on small Wilson loops with massive fermions to $`𝒪(g^4)`$. We also include some findings that apply to the Iwasaki and Symanzik-Weisz gauge actions.

### II.1 The Method

We consider $`SU(N)`$ lattice gauge theory on an infinite four dimensional isotropic hyper-cubic lattice with lattice spacing $`a`$. An oriented closed curve $`𝒞`$ of length $`na`$ touches a sequence of sites, $`\{x_i:x_i/a^4\}`$ such that $`x_{i+1}x_i\{\pm a\widehat{\mu }:\mu =1,\mathrm{},4\}`$, $`x_{n+1}=x_1`$. The Wilson loop $`W(𝒞)`$ around such a curve is the expectation value of the path ordered product of gauge links,

$$W(𝒞)=\frac{1}{N}\text{tr}\left[\underset{x_i𝒞}{}U_{x_i,\mu _i}\right],$$
(1)
where $`\mu _i\{\pm 1,\mathrm{},\pm 4\}`$ denotes the direction indicated by $`x_{i+1}x_i`$ and $`U_{x,\mu }=U_{x\widehat{\mu },\mu }^{}`$.

We write $`W(𝐑,T)`$ for a “rectangular” Wilson loop where the curve $`𝒞`$ contains two opposite lines with an extent of $`T`$ lattice units pointing into the time ($`\widehat{4}`$) direction, that are separated by a spatial distance $`𝐑a`$. The smallest such Wilson loop is an $`a\times a`$ square, the plaquette $`\mathrm{}=W(1,1)`$ which, appropriately normalized, is the expectation value of the Wilson gauge action,

$$S=\beta \underset{x,\mu >\nu }{}\frac{1}{N}\text{Re}\text{tr}U_{x,\mu \nu },$$
(2)
where $`\beta =2N/g^2`$ and $`U_{x,\mu \nu }=U_{x,\mu }U_{x+a\widehat{\mu },\nu }U_{x+a\widehat{\nu },\mu }^{}U_{x,\nu }^{}`$. Unless stated otherwise this is the gauge action which we combine either with the Wilson, SW or the KS fermionic actions.

Writing $`U_l=e^{igaA_l}`$, where $`l=(x,\mu )`$ denotes the link connecting $`x`$ with $`x+a\widehat{\mu }`$ and $`A_l`$ is a short hand notation for $`A_\mu (x+\frac{a}{2}\widehat{\mu })`$, and expanding the exponentials within Eq. (1) one obtains,

$$W(𝒞)=1g^2W_2g^3W_3g^4W_4\mathrm{},$$
(3)
with Weisz:1984bn ,

$$W_2=\frac{a^2}{2N}\text{tr}T^aT^b\underset{l_1,l_2}{}A_{l_1}^aA_{l_2}^b,$$
(4)
$`W_3`$ $`=`$ $`{\displaystyle \frac{ia^3}{6N}}\text{tr}T^aT^bT^c{\displaystyle \underset{l_1,l_2,l_3}{}}A_{l_1}^aA_{l_2}^bA_{l_3}^c`$
$`+3{\displaystyle \underset{l_1<l_2<l_3}{}}(A_{l_1}^aA_{l_2}^bA_{l_3}^cA_{l_3}^aA_{l_2}^bA_{l_1}^c),`$
$`W_4`$ $`=`$ $`{\displaystyle \frac{a^4}{24N}}\text{tr}T^aT^bT^cT^d24{\displaystyle \underset{l_1<l_2<l_3<l_4}{}}A_{l_1}^aA_{l_2}^bA_{l_3}^cA_{l_4}^d`$
$`+12{\displaystyle \underset{l_1<l_2<l_3}{}}(A_{l_1}^aA_{l_1}^bA_{l_2}^cA_{l_3}^d+A_{l_1}^aA_{l_2}^bA_{l_2}^cA_{l_3}^d`$
$`+A_{l_1}^aA_{l_2}^bA_{l_3}^cA_{l_3}^d)`$
$`+{\displaystyle \underset{l_1<l_2}{}}(6A_{l_1}^aA_{l_1}^bA_{l_2}^cA_{l_2}^d+4A_{l_1}^aA_{l_1}^bA_{l_1}^cA_{l_2}^d`$
$`+4A_{l_1}^aA_{l_2}^bA_{l_2}^cA_{l_2}^d)`$
$`+{\displaystyle \underset{l_1}{}}A_{l_1}^aA_{l_1}^bA_{l_1}^cA_{l_1}^d.`$
The expectation values in the above formulae depend on $`g^2`$ and are sums of gluonic position space $`n`$-point functions. The leading order \[$`𝒪(g^2)`$\] contribution to the Wilson loop is the sum of tree-level two-point functions depicted in Fig. 1. This is easily expressed in momentum space as a sum of phase factors multiplying the momentum space propagator:

$`W_2`$ $`=`$ $`{\displaystyle \underset{l_1,l_2}{}}{\displaystyle _{\pi /a}^{\pi /a}}{\displaystyle \frac{d^4k}{(2\pi )^4}}{\displaystyle \frac{a^2}{2N}}\text{tr}T^aT^be^{ik[x_2x_1+\frac{a}{2}(\widehat{\mu }_2\widehat{\mu }_1)]}`$ (7)
$`\times A_{\mu _1}^a(k)A_{\mu _2}^b(k)`$
$`=`$ $`W_2^{(0)}+g^2W_2^{(2)}+𝒪(g^4).`$
For the Wilson gauge action the gluon propagator reads,

$$A_{\mu _1}^a(k)A_{\mu _2}^b(k)=\frac{\delta ^{ab}(\delta _{\mu _1,\mu _2}\delta _{\mu _1,\mu _2})}{\widehat{k}^2}+\mathrm{\Pi }_{\mu _1\mu _2}^{ab}(k),$$
(8)
where $`\widehat{k}_\mu =2a^1\mathrm{sin}\left(k_\mu a/2\right)`$ and $`\mathrm{\Pi }`$ is of order $`g^2`$. The colour trace then results in the Casimir factor of the fundamental representation, $`C_F=\frac{1}{N}\text{tr}T^aT^a=(N^21)/(2N)`$. All formulae in this Section as well as in Sec. III below apply to the general $`SU(N)`$ case unless stated otherwise. In addition the constants $`C_F`$ are factorized out in such a way that the results apply to external sources in any representation $`D`$ of the gauge group, under the replacement $`C_FC_D`$.

In explicitly summing the above phase factors one easily obtains the well known tree-level expressions for planar Wilson loops. The vacuum polarization tensor $`\mathrm{\Pi }`$ has to be considered to leading order for an $`𝒪(g^4)`$ evaluation of Wilson loops which requires $`W_2^{(2)}`$. The relevant graphs for the gluonic sector are depicted in Fig. 2. The fourth diagram involving a two-gluon-two-ghost vertex (which is lattice specific) as well as a lattice contribution to the four-gluon vertex of the third diagram, result in a lattice tadpole correction to $`W_2^{(2)}`$ that is exactly $`(2N^23)/(24N)`$ times the tree-level value $`W_2^{(0)}`$. The two additional diagrams of Fig. 3 have to be considered in the case of dynamical fermions.

To the order at which we are working throughout this paper, $`W_3`$ only contributes through its contraction with the three gluon vertex, depicted in Fig. 4, and we write

$$W_3=gW_3^{(1)}+𝒪(g^3).$$
(9)
In momentum space the phase factor associated with a typical entry in the sum Eq. (II.1) is

$`{\displaystyle _{\pi /a}^{\pi /a}}{\displaystyle \frac{d^4kd^4k^{}}{(2\pi )^8}}{\displaystyle \frac{i}{N}}\text{tr}T^aT^bT^ce^{ik[x_2x_1+\frac{a}{2}(\widehat{\mu }_2\widehat{\mu }_1)]}`$ (10)
$`\times e^{ik^{}[x_3x_1+\frac{a}{2}(\widehat{\mu }_3\widehat{\mu }_1)]}A_{\mu _1}^a(kk^{})A_{\mu _2}^b(k)A_{\mu _3}^c(k^{}).`$
The term

$$W_4=W_4^{(0)}+𝒪(g^2)$$
(11)
only contributes to the Wilson loop at $`𝒪(g^4)`$ through its contraction with double gluon exchanges, as depicted in Fig. 5. The generic phase factor associated with this diagram is given below, where due to the lack of an appropriately symmetrized contracting vertex we perform the Wick contraction manually, preserving the ordering of colour factors:

$$_{\pi /a}^{\pi /a}\frac{d^4k}{(2\pi )^4}\frac{d^4k^{}}{(2\pi )^4}\frac{i}{N}\text{tr}T^aT^bT^cT^d\left\{\begin{array}{c}e^{ik[x_2x_1+\frac{a}{2}(\widehat{\mu }_2\widehat{\mu }_1)]}e^{ik^{}[x_4x_3+\frac{a}{2}(\widehat{\mu }_4\widehat{\mu }_3)]}A_{\mu _1}^a(k)A_{\mu _2}^b(k)A_{\mu _3}^c(k^{})A_{\mu _4}^d(k^{})\\ e^{ik[x_3x_1+\frac{a}{2}(\widehat{\mu }_3\widehat{\mu }_1)]}e^{ik^{}[x_4x_2+\frac{a}{2}(\widehat{\mu }_4\widehat{\mu }_2)]}A_{\mu _1}^a(k)A_{\mu _2}^b(k^{})A_{\mu _3}^c(k)A_{\mu _4}^d(k^{})\\ e^{ik[x_3x_2+\frac{a}{2}(\widehat{\mu }_3\widehat{\mu }_2)]}e^{ik^{}[x_4x_1+\frac{a}{2}(\widehat{\mu }_4\widehat{\mu }_1)]}A_{\mu _1}^a(k^{})A_{\mu _2}^b(k)A_{\mu _3}^c(k)A_{\mu _4}^d(k^{})\end{array}\right\}.$$
(12)
After evaluating colour traces, this term can be split into an “Abelian” and a “non-Abelian” part. The Abelian part is proportional to $`C_F^2/2`$ and can be obtained by expanding $`\mathrm{exp}(g^2W_2)`$. The remainder is specific to non-Abelian gauge theories and proportional to $`C_FN/2`$, the same colour factor that also accompanies the vacuum polarization diagrams.

The integrands were coded with the help of a C++ package written by one of the authors, and then integrated using VEGAS Lepage:1978sw . For finite Wilson loops to $`𝒪(g^2)`$ and $`𝒪(g^4)`$ we compute 4 and 8 dimensional momentum integrals, respectively. In the case of the potential one can analytically perform the limit $`T\mathrm{}`$. This allows us to compute the momentum integrals in one less dimension, leaving 3 and 7 dimensional integrals.

The fermionic Wilson-SW and KS vacuum polarization graphs have been computed after a further analytical integration over the 4-component of the internal loop momentum, leaving a 6 dimensional integration to be performed in the cases of the potential and the static self energy.

It is relatively easy to automate the generation of the phase factors that are required for arbitrary curves $`𝒞`$. We generate the perturbative expansion for both small rectangular Wilson loops and for the non-planar loops shown in Fig. 6. We refer to these non-planar loops as the chair and the (three-dimensional) “parallelogram”; both are contained within the unit cube and are used within certain improved gauge actions.

### II.2 Definitions

We write the tree-level contribution to the expansion of the Wilson loop with the Wilson gauge action as,

$$W_2^{(0)}=C_FW_T.$$
(13)
For Wilson and SW fermions with the Wilson gauge action we write the vacuum polarization insertions as,

$$W_2^{(2)}=C_F\left(\frac{N}{2}W_\mathrm{\Pi }+\frac{2N^23}{24N}W_T+n_fX_f\right).$$
(14)
For Wilson-SW fermions we define,

$$X_f=X_f^{(0)}+X_f^{(1)}c_{SW}+X_f^{(2)}c_{SW}^2,$$
(15)
where $`c_{SW}=0`$ in the Wilson case and $`c_{SW}=1+𝒪(g^2)`$ in the SW case. The $`X_f`$ will depend on the quark mass in lattice units $`ma`$. In the case of Wilson-SW fermions and the order of perturbation theory we are working at we have, $`ma=m_0a=(\kappa ^1\kappa _c^1)/2`$. As anticipated above, we have factorized out a lattice tadpole term that is proportional to $`W_T`$. For KS fermions we define analogously,

$$X_f=X_f^{(KS)}.$$
(16)
In this case $`n_f`$ has to be a multiple of four. For each of the above quantities we may write a similar expansion in terms of $`\alpha _L=g^2/(4\pi )`$, rather than $`g^2`$, with lower case coefficients $`w`$ and $`x`$, where $`w_T=4\pi W_T`$, $`w_\mathrm{\Pi }=(4\pi )^2W_\mathrm{\Pi }`$, $`x=(4\pi )^2X`$ and so forth.

We denote the spider contribution by,

$$W_3^{(1)}=\frac{C_FN}{2}W_S.$$
(17)
The double gluon exchange can be written as the sum of Abelian and non-Abelian pieces, where the Abelian piece is the next term in the exponentiation of the tree-level contribution,

$$W_4^{(0)}=\frac{C_F^2}{2}W_T^2+\frac{C_FN}{2}W_{NA}.$$
(18)
In combining the above Eqs. (13) – (18) with Eqs. (3), (7), (9) and (11), we obtain the expansion of the Wilson loop,

$$W=1W_{\text{LO}}g^2\left(W_{\text{NLO}}+C_Fn_fX_f\right)g^4,$$
(19)
where

$`W_{\text{LO}}`$ $`=`$ $`C_FW_T,`$ (20)
$`W_{\text{NLO}}`$ $`=`$ $`C_F\left({\displaystyle \frac{N}{2}}W_\mathrm{\Sigma }+{\displaystyle \frac{2N^23}{24N}}W_T{\displaystyle \frac{C_F}{2}}W_T^2\right),`$ (21)
$`W_\mathrm{\Sigma }`$ $`=`$ $`W_\mathrm{\Pi }+W_{NA}+W_S`$ (22)
The $`W_T^2`$ term within Eq. (21) originates from the exponentiation of the Abelian part of the one gluon exchange, Eq. (18). It is cancelled in the Taylor expansion of $`\mathrm{ln}W`$, which turns out to be proportional to $`C_F`$, at least to $`𝒪(g^4)`$. The latter observation, which is referred to as the “Casimir scaling” hypothesis in the literature (see e.g. Refs. Bali:2000un ; Deldar:2000vi ), implies that, for given $`N`$ and $`n_f`$, $`\mathrm{ln}W`$ only depends on the representation of the static colour sources, $`D`$, through its proportionality to the corresponding Casimir charge $`C_D=(\text{tr}T_D^aT_D^a)/\text{dim}_D`$. Note that in the case of the fundamental representation, by using the relation $`C_F=(N^21)/(2N)`$, Eq. (21) can be rearranged into the form,

$$W_{\text{NLO}}=(N^21)\left(W_a+\frac{1}{N^2}W_b\right),$$
(23)
where

$`W_a`$ $`=`$ $`{\displaystyle \frac{1}{24}}\left(6W_\mathrm{\Sigma }+W_T3W_T^2\right),`$ (24)
$`W_b`$ $`=`$ $`{\displaystyle \frac{1}{16}}\left(W_T+2W_T^2\right).`$ (25)
For improved gluonic actions Eq. (21) does not apply and hence $`W_a`$ and $`W_b`$ can take somewhat different forms. The factorization Eq. (23) is frequently employed throughout the literature, e.g. in Refs. DiGiacomo:1981wt ; Iso:1987dj ; Alles:1994dn .

### II.3 Small Wilson loop results

The pure gauge results on small Wilson loops for the Wilson gluonic action are displayed in Tab. 1. We reproduce the known plaquette results of Di Giacomo and Rossi DiGiacomo:1981wt  and Allés et al. Alles:1998is 111In fact the result of this reference is more precise than ours: they obtain $`W_{\text{NLO}}=0.033910993(1)`$ while we quote $`W_{\text{NLO}}=0.033911(1)`$ in the table. and the small planar Wilson loop results of Wohlert, Weisz and Wetzel Weisz:1984bn . Boldface values have been calculated by us for the first time. We also calculated tree-level Wilson loops for the Symanzik-Weisz (SyW) Weisz:1983zw  and Iwasaki (I) Iwasaki:1984cj  actions. These results are displayed in Tabs. 2 and 3, respectively, together with the $`𝒪(g^4)`$ values of Iso and Sakai Iso:1987dj  (italicized) which we have not validated ourselves.

We are able to reproduce the known contributions to the plaquette for massless KS Heller:1985eq  and Wilson Hamber:1983ft  quarks with increased precision. The calculation of the SW contribution222Our result has already been used in Refs. Marcantonio:2001fc ; Booth:2001qp ; Bali:2001yh ., as well as all results for massive quarks are new. We also calculate the one-loop $`n_f`$ piece of the plaquette with the Iwasaki gauge action which has been used, for instance, in the dynamical fermion simulations of the CP-PACS group Yoshie:2000wd . We display a selection of $`m=0`$ results in Tab. 4. In Fig. 7 we show the mass dependence of $`X_f`$ for the case of the Wilson gluonic action, combined with all three different fermionic actions. The open symbols correspond to the respective $`m=0`$ limits. The numbers that are plotted in the figure are also displayed in Tab. 5.

In the region $`ma<0.25`$, which is most relevant for lattice simulations, this dependence is relatively weak and can be parametrized by,

$$X_f(ma)=X_f(0)+bma+c(ma)^2,$$
(26)
with $`X_f(0)`$ as in the first row of Tab. 5 and $`b0`$, 0.0024 and 0.0024 and $`c0.0017`$, $`0.0023`$ and $`0.0017`$ for KS, Wilson and SW fermions, respectively. As expected the mass dependence remains weak for any Wilson loop with linear dimensions that are small compared to the inverse quark mass.

## III The perturbative static potential

After a historical survey of the existing literature we will briefly explain our method and the notations that we adopt, before presenting results on the $`𝒪(g^4)`$ potential. We will then discuss limiting cases, including the continuum limit and lattice artefacts as well as mass dependent terms in the conversion between lattice and $`\overline{MS}`$ schemes at finite lattice spacing. We conclude with paragraphs on the perturbative $`\beta `$-shift, boosted lattice perturbation theory and the static source self energy.

### III.1 The method and definitions

The one-loop potential has first been calculated for pure $`SU(2)`$ gauge theory by Müller and Rühl MullerRuhl . Subsequently, a closed formula for pure $`SU(N)`$ gauge theory with the standard Wilson as well as Symanzik improved gauge actions was derived by Weisz and Wohlert Weisz:1984bn . The coefficients of the perturbative expansion have first been evaluated numerically by Altevogt and Gutbrod Altevogt:1995cj  for $`SU(2)`$ pure gauge theory with Wilson glue on isotropic and anisotropic lattices.

Here we restrict our calculations to the case of the Wilson gauge action on an isotropic lattice but include massive Wilson, SW and KS fermions as well as $`N>2`$. We also improve the numerical precision with respect to earlier results and incorporate off-axis lattice separations $`𝐑`$ of the colour sources. We define the static potential,

$`aV(𝐑a)`$ $`=`$ $`\underset{T\mathrm{}}{lim}{\displaystyle \frac{d\mathrm{ln}W(𝐑,T)}{dT}}`$ (27)
$`=`$ $`V_1(𝐑)g^2+V_2(𝐑)g^4+\mathrm{}`$ (28)
$`=`$ $`v_1(𝐑)\alpha _L+v_2(𝐑)\alpha _L^2,`$ (29)
where $`\alpha _L=g^2/(4\pi )`$, $`v_i(𝐑)=(4\pi )^iV_i(𝐑)`$. $`W(𝐑,T)`$ denotes a (generalized) rectangular Wilson loop with the temporal extent of $`T`$ lattice units. The expansion of the potential in terms of $`g^2`$ is suitable for calculations in lattice perturbation theory and for the comparison with data from lattice simulations, while the expansion in terms of $`\alpha _L`$ turns out to be more convenient to relate our results to those obtained in a continuum scheme.

The only gluon exchanges that contribute to the $`T`$ dependence of $`W(𝐑,T)`$ in the limit $`T\mathrm{}`$ \[at least up to $`𝒪(g^6)`$\] are those between temporal lines of a Wilson loop. After exploiting translational invariance in time, Eqs. (7) and (8) result in,

$$V_1(𝐑)=C_F\frac{d^3q}{(2\pi )^3}\frac{2\mathrm{sin}^2\left(\frac{\mathrm{𝐪𝐑}a}{2}\right)}{\widehat{𝐪}^2}.$$
(30)
The $`𝒪(g^4)`$ coefficient $`V_2`$ consists of two parts: one contribution from $`W_2`$, i.e. the vacuum polarization tensor $`\mathrm{\Pi }`$, inserted into Eq. (30), and a second contribution that originates from the “non-Abelian” part $`W_{NA}`$ of $`W_4`$ \[Eq. (18)\]. The spider diagrams, $`W_3`$, do not contribute to $`V(𝐑)`$ at this order. A closed formula for $`V_2`$ can be found in Ref. Weisz:1984bn , and we extend this calculation by incorporating the fermionic contributions to the vacuum polarization tensor. In our calculation we produce $`V_2`$ automatically from the Feynman rules and find agreement with the published analytic form.

The potential can be factorized into an interaction part and a part that is due to the static source self energy, $`V_S=2\delta m_{\text{stat}}`$,

$$V(𝐑a)=V_{\text{int}}(𝐑a)+V_S.$$
(31)
Since in perturbation theory $`V_{\text{int}}0`$ as $`r\mathrm{}`$ we can identify the self energy $`V_S`$ with $`lim_r\mathrm{}V(r)`$ and write,

$`aV_S`$ $`=`$ $`V_1(\mathrm{})g^2+V_2(\mathrm{})g^4+\mathrm{},`$ (32)
$`aV_{\text{int}}(𝐑a)`$ $`=`$ $`V_{L,1}(𝐑)g^2+V_{L,2}(𝐑)g^4+\mathrm{},`$ (33)
where $`V_{L,i}(𝐑)=V_i(𝐑)V_i(\mathrm{})`$. The limit $`R\mathrm{}`$ can easily be realized by replacing the factors $`\mathrm{sin}^2[\mathrm{𝐪𝐑}a/2]=[1\mathrm{cos}(\mathrm{𝐪𝐑}a)]/2`$ within the external 3-dimensional Fourier transformations in the expressions for $`V_1`$ \[cf. Eq. (30)\] and $`V_2`$ by the constant value $`1/2`$.

We write,

$`V_1(𝐑)`$ $`=`$ $`C_FV_T(𝐑)`$ (34)
$`V_2(𝐑)`$ $`=`$ $`C_F\{{\displaystyle \frac{N}{2}}[V_\mathrm{\Pi }(𝐑)+V_{NA}(𝐑)]`$
$`+{\displaystyle \frac{2N^23}{24N}}V_T(𝐑)+n_fV_f(𝐑)\},`$
where $`V_\mathrm{\Pi }`$ denotes the gluonic vacuum polarization and $`V_{NA}`$ is a contribution, specific to non-Abelian gauge theories, where the ordering of gluon vertices along the Wilson loop has to be considered. In the case of KS quarks we denote the fermionic vacuum polarization contribution to the potential by $`V_f^{(KS)}`$ while for Wilson-SW quarks we split this contribution into three parts,

$$V_f(𝐑)=V_f^{(0)}(𝐑)+c_{SW}V_f^{(1)}(𝐑)+c_{SW}^2V_f^{(2)}(𝐑).$$
(36)
### III.2 Results on the potential and violations of rotational symmetry

In Tab. 6 we display $`V_T`$, $`V_\mathrm{\Pi }`$ and $`V_{NA}`$ as well as $`V_2`$ for $`N=3,n_f=0`$ for small on- and off-axis distances $`𝐑`$. At the origin all $`V`$ are zero and $`V_T(1,0,0)=1/6`$. In Tab. 7 we display the corresponding results on the different $`V_f`$’s for massless quarks while in Tabs. 89 the results for distances $`R3`$ as well as for $`R\mathrm{}`$ are shown for quarks of various masses.

In Fig. 8 we separately display the different $`𝒪(g^4)`$ contributions to the quantity $`RV_2(𝐑)`$ as a function of $`R`$. The factor $`R`$ results in a cancellation of the leading order Coulomb behaviour. Note that the lattice tadpole contribution to the vacuum polarization (solid circles) is numerically equally important as the sum of the remaining vacuum polarization diagrams and the non-Abelian contribution. As one would naïvely expect, the effect of fermions (open squares) goes into the opposite direction, relative to the pure gauge contributions. In a continuum calculation using dimensional regularization $`RV_2(𝐑)`$ contains only a logarithmic dependence on R, but in the lattice calculation the $`R`$ dependence of the points is dominated by a linear piece arising from the static source self energy contribution to $`V_2(𝐑)`$, $`V_{S,2}=V_2(\mathrm{})`$. In Fig. 9 we subtract this term, before multiplying the result with $`R`$. Note the logarithmic scale. The static source self-energy was previously known to order $`\alpha ^2`$ with massless Wilson and SW quarks Martinelli:1999vt  and to $`𝒪(\alpha ^3)`$ in the pure gauge theory DiRenzo:2001nd ; Trottier:2001vj . The mass-dependence, as well as the results for KS sea quarks, are new.

In the limit $`ra`$, i.e. $`R=r/a1`$ rotational symmetry should be restored such that we are able to compare our result to calculations performed in a continuum scheme (continuous curves in Figs. 8 and 9). The potential has been calculated in the $`\overline{MS}`$ scheme to order $`\alpha ^3`$ Schroder:1999vy ; Melles:2000dq . The conversion between the $`\overline{MS}`$ and the lattice scheme has been worked out to this order too for Wilson-SW fermions Christou:1998ws ; Bode:2001uz  (for KS fermions only to $`𝒪(\alpha ^2)`$ Weisz:1981pu ; Luscher:1995zz ), and $`V_S`$ vanishes by definition to all orders in dimensional regularization. The form of the large $`R`$ expectation is worked out in Appendix A.4 \[Eqs. (192) – (194)\]. Here we display the parametrization of the curves in the limit of massless sea quarks:

$`rV_{L,2}(𝐫/a)`$ $`\stackrel{r/a\mathrm{}}{}`$ $`rV_{c,2}(r)`$
$`=`$ $`{\displaystyle \frac{C_F}{4\pi }}\left[b_1+a_1^R+2\beta _0\mathrm{ln}(r/a)\right].`$
The constants $`b_1`$ and $`a_1^R=a_1+2\gamma \beta _0`$ are defined in Eqs. (134) and (154), respectively. $`b_1`$ originates from the conversion between the $`\overline{MS}`$ and the lattice scheme, Eq. (133). The logarithmic running is proportional to the $`\beta `$-function coefficient $`\beta _0=(11N2n_f)/(12\pi )`$.

We can define an effective Coulomb coupling $`C_F\alpha _R`$ from the potential:

$`r\left[V(r)V_S\right]`$ $`=`$ $`C_F\alpha _R\left(r^1\right)`$ (38)
$`=`$ $`C_F\alpha _LRV_{L,2}(R)g^4.`$ (39)
Note that in this limit, $`C_F\alpha _L=RV_{L,1}(R)`$. In Fig. 9 the enhancement of this effective coupling, $`C_F\alpha _R`$, is demonstrated, relative to the tree-level value as a function of $`R`$. The points are indeed consistent with the logarithmic running proportional to $`2\beta _0C_F\mathrm{ln}R/(4\pi )`$ that is expected from the QCD $`\beta `$-function, Eq. (121) \[as well as from Eq. (III.2) above\]. The lattice tadpole terms do not contribute to this running but renormalize the overall value of the coupling. Lattice simulations are usually performed around $`g^21`$, where a fit to quenched lattice data Bali:2000vr ; Bali:2000gf  yields $`eC_F\alpha _R[1/(3a)]0.3`$. The tree-level expectation in this case, however, is substantially smaller: $`C_F\alpha _R=C_F\alpha _L0.106`$. As one can read off from the figure, the one-loop value at $`R=3`$ adds about 0.08 to this, but still the non-perturbative result is underestimated by perturbation theory in terms of the bare coupling $`g^2`$ by more than 30 %. We will discuss the improvement resulting from the use of so-called boosted perturbation theory in Sec. IV.3 below.

At small $`R`$ values the lattice results scatter around the continuous curves. We shall investigate these violations of rotational symmetry in more detail: in the top half of Fig. 10 we plot ratios of lattice and continuum $`𝒪(g^2)`$ (open squares, $`V_{L,1}/V_{c,1}1`$) and $`𝒪(g^4)`$ (open circles, $`V_{L,2}/V_{c,2}1`$) contributions for on-axis as well as for some off-axis distances $`𝐑`$ for $`n_f=0`$; the relative violations of rotational symmetry become smaller with bigger $`R`$ and the one-loop violations are smaller than the tree-level ones. In the lower part of the figure the ratio of the two ratios is displayed, which amounts to replacing the continuum $`1/R`$ term that multiplies the logarithmic running of $`V_{L,2}`$ by a lattice “$`[1/𝐑]`$” function. In doing so we isolate those lattice artefacts that appear specifically at order $`g^4`$; for instance, the lattice tadpole contributions cancel from this combination. The sign is opposite to that of the tree-level differences, which explains the weaker relative violations at $`𝒪(g^4)`$ in the comparison with $`𝒪(g^2)`$ in the upper part of the figure.

### III.3 The continuum limit

The limit of large $`R=r/a`$ has two interpretations: taken at a fixed lattice spacing $`a`$ (sufficiently small for perturbation theory still to be reliable at a distance $`r`$) continuum perturbative predictions should be met. We have already discussed this scenario above and indeed demonstrated this agreement at large $`R`$ in Figs. 9 and 10. On the other hand taking the limit of large $`R`$ at fixed physical $`r`$ corresponds to the continuum limit $`a0`$. In Fig. 11 we investigate the approach to the continuum limit by displaying ratios of lattice and continuum perturbation theory results $`V_{L,i}/V_{c,i}`$ for $`i=1,2`$ versus $`1/R^2=a^2/r^2`$ for pure $`SU(3)`$ gauge theory and on-axis separations $`R2`$. For the Wilson gauge action we expect lattice artefacts to be a polynomial in $`a^2`$ and indeed no linear term is found. The solid curves represent fits that are quadratic plus quartic in the lattice spacing $`a`$. For off-axis separations such as $`𝐑(1,1,0)`$ we find a similar picture, but with different $`a^2`$ and $`a^4`$ coefficients.

The same comparison is displayed for the fermionic contributions $`V_{L,f}(𝐑)=V_f(𝐑)V_f(\mathrm{})`$ and $`V_{c,f}(r)=lim_{a0}V_{L,f}(r/a)`$ alone for massless quarks in Fig. 12. The fit curves are quadratic plus cubic for SW ($`c_{SW}=1`$ to this order in perturbation theory) and KS fermions while for Wilson fermions we indeed require the expected linear term and attempt a linear plus quadratic fit. The coefficient of the linear term turns out to be so small that the quadratic term already dominates for distances $`a/r>0.1`$. We see that in spite of the more favourable functional dependence on $`a`$ the perturbative lattice artefacts of SW fermions are numerically bigger over the whole displayed $`R`$ range than those of Wilson fermions. However, the KS action “outperforms” both Wilson and SW for this observable. As the quark mass is increased the violations of the continuum symmetry become more pronounced in the Wilson case while for the other two actions the change is only small, as a comparison of Fig. 12 with Fig. 13, that corresponds to a quark mass $`ma=0.1`$, reveals: “improvement” enables one to simulate on coarser lattices at the same physical quark mass.

From Fig. 9 it is obvious that for realistic values of $`n_f`$ the potential is dominated by the gluonic contributions and the additional violations of rotational symmetry due to sea quarks that we have just discussed are an interesting but numerically subleading effect.

We shall now return to the limit where continuum perturbation theory is met, i.e. we investigate the behaviour at fixed (small) $`a`$ and large $`r=Ra`$. In the massless case this limit and the limit of large $`R=r/a`$ at fixed $`r`$ discussed above are equivalent, up to non-perturbative effects. However, as soon as a second scale $`ma>0`$ is introduced, the situation $`r>m^1`$ becomes distinguishable from $`r<m^1`$. In Fig. 14 we again display the self energy subtracted fermionic contribution $`V_{L,f}(𝐑)`$, this time multiplied by $`R`$, in analogy to Fig. 9. This combination, multiplied by $`C_Fn_f`$ isolates the shift that is induced onto an effective coupling $`C_F\alpha _R`$, due to sea quarks. The logarithmic slope of the three curves that are displayed in the figure is determined by the fermionic contribution to $`\beta _0`$ and is therefore universal. The inclusion of sea quarks not only reduces $`\beta _0`$, relative to the pure gauge case, and hence slows down the running of $`\alpha _R`$ but it also decreases the absolute value of the $`𝒪(g^4)`$ coefficient. In the case of massive quarks, however, this effect is (over)compensated by an increase in the coupling $`g^2`$, relative to the pure gauge case, if we require the same physics at scales $`rm^1`$ where the sea quarks decouple. This will be discussed in Sec. IV.3.3 below.

While for massless flavours the coupling runs logarithmically, for massive quarks the running gradually switches off around333The constant $`C_05.2`$ Brodsky:1999fr  is introduced in Eq. (156). $`R0.3/(ma)1/(\sqrt{C_0}ma)`$: physics at length scales $`Ram^1`$ is insensitive to the presence of massive sea quarks, at least in perturbation theory. In Fig. 15 we again compare the contributions to the effective Coulomb couplings resulting from the three different fermionic actions but for $`ma=0.1`$. At this mass value Wilson and KS results accidentally happen to be very close to each other. In Fig. 16 we investigate the mass dependence of the fermionic term somewhat more systematically by studying the example of Wilson quarks alone. In order to relate the lattice to a continuum scheme at large $`R`$, we not only have to subtract the (sea quark mass-dependent) static source self energy contribution to eliminate terms proportional to $`R`$ from the figure but also the vertical offset, $`K_1(ma)`$, at small $`R`$ changes due to a mass dependent term that appears in the conversion factor $`b_1`$ between the $`\overline{MS}`$ and the lattice scheme at finite $`a`$. We will discuss this effect in some detail in Sec. III.4 below.

In the limit $`rm^1`$ all curves indeed approach a constant value: the running of the effective coupling is not affected by the presence of massive sea quark flavours anymore. The resulting potential is the same as that in the pure gauge case, at least in perturbation theory, albeit with a different overall normalization of the effective Coulomb coupling. We shall discuss this shift of the QCD coupling constant at a given scale and the possibility of matching to the quenched theory in Sec. III.5 below.

While at large $`R`$ the lattice data and the continuum curve nicely coincide with each other \[once the offset $`K_1(ma)`$ has been determined from the large $`R`$ data and subtracted\], at small $`R`$ we observe large $`𝒪(am)`$ lattice artefacts, in addition to mass independent $`𝒪(a)`$ effects: even at the relatively small value $`ma=0.03`$ the lattice points lie systematically below the large $`R`$ expectation. Fig. 15 confirms this qualitative pattern for the other two fermionic actions.

### III.4 $`\mathrm{\Delta }K_1`$

As we have seen above the choice of lattice action not only affects the quality of rotational symmetry at small $`R`$ but also the overall normalization. The $`\overline{MS}`$ scheme is related to the lattice scheme via

$$\alpha _{\overline{MS}}(a^1)=\alpha _L+b_1\alpha ^2+\mathrm{},$$
(40)
with the conversion factor

$$b_1=\pi /(2N)+k_1N+K_1(ma)n_f.$$
(41)
The numerical constant $`k_1`$ is known for a variety of gluonic actions and $`K_1(0)`$ is known for Wilson, SW and KS quarks and independent of the gluonic action, cf. Eqs. (A.2) – (137).

The coefficients of the $`\beta `$-function in the $`\overline{MS}`$-scheme as well as in the lattice scheme in the continuum limit do not depend on the quark mass: both are “mass-independent” renormalization schemes. In lattice simulations it is often worthwhile to analyse quantities prior to an extrapolation to the continuum limit. One such example is determinations of the QCD running coupling from small Wilson loops El-Khadra:1992vn ; davi1 ; davi2 ; Spitz:1999tu ; Booth:2001qp ; davi3  that are ill-defined in this limit. We will discuss this technique in Sec. IV.2 below. Another example is simulations within the framework of an effective field theory like NRQCD Thacker:1991bm  that requires a finite momentum cut-off.

At finite $`a`$ and $`m>0`$ the lattice scheme becomes mass-dependent, as indicated by the argument of the function,

$$K_1(ma)=K_1(0)+\mathrm{\Delta }K_1(ma),$$
(42)
within Eq. (41). In this case universality is lost and the coefficients of the perturbative $`\beta `$-function acquire additional contributions that, like $`\mathrm{\Delta }K_1`$, will in general depend on the dimensionful observable that is studied:

$$\beta _0^{\overline{MS}}\beta _0^L=n_f\frac{d\mathrm{\Delta }K_1(ma)}{d\mathrm{ln}a^2}.$$
(43)
By definition, $`\mathrm{\Delta }K_1(0)=0`$. However, the behaviour of $`\mathrm{\Delta }K_1`$ is also constrained in the limit $`ma\mathrm{}`$ where the sea quarks decouple and therefore $`\beta _0^{\overline{MS}}=(11N)/(12\pi )`$ while $`\beta _0^L=\beta _0^{\overline{MS}}n_f/(6\pi )`$, i.e.,

$$\mathrm{\Delta }K_1\stackrel{ma\mathrm{}}{}\frac{1}{3\pi }\mathrm{ln}(Dma),$$
(44)
with a dimensionless constant $`D`$ that we calculate in Sec. III.5 below.

We have seen above that in the limit $`rm^1`$ the behaviour of the massless theory is emulated while the running of the potential for $`rm^1`$ is effectively quenched: the only scale that is relevant in perturbation theory in these limits is the distance $`r`$ (or, equivalently, momentum $`qr^1`$) and therefore lattice artefacts disappear like $`(a/r)^\nu `$ \[or $`(aq)^\nu `$\] with some positive integer power $`\nu `$, whose value depends on the lattice action, as $`r\mathrm{}`$. However, in the intermediate range of quark masses lattice corrections $`(am)^\nu `$ become relevant and the universality of the $`\beta `$-function is lost.

The function $`K_1(ma)`$ of Eq. (134) can in principle be read off from figures such as Figs. 15 and 16 up to short distance lattice artefacts: $`K_1(ma)=C_F/(4\pi )V_{L,2}(1)\sqrt{C_0}ma+𝒪[(ma)^2,a^\nu ]`$ \[cf. Eqs. (III.2), (171) and (196)\]. For the (improved) SW action not only the lattice artefacts are much more pronounced and qualitatively different from the Wilson and KS cases but also this overall offset $`K_1`$ is enhanced (although not its mass dependence).

We determine $`K_1(ma)`$ from two parameter fits,

$`{\displaystyle \frac{4\pi R}{C_F}}V_{L,2}`$ $`=`$ $`A_0+A_1\mathrm{ln}R+{\displaystyle \frac{n_f}{3\pi }}\text{Ein}\left(C_0^{1/2}maR\right)`$ (45)
$`+`$ $`n_f\mathrm{\Delta }K_1(ma)+{\displaystyle \frac{c}{R^\nu }},`$
to the $`R>9`$ on-axis data points, with one redundant parameter $`c`$ whose role is to parametrize any residual lattice artefacts. The above functional form with the constants $`A_0`$ and $`A_1`$ of Eqs. (A.4) – (195) is motivated by Eq. (192) of Appendix A.4. The special function $`\text{Ein}(x)`$ denotes the normalized exponential integral defined in Eq. (172).

All fits turned out to be stable within statistical errors under variations of the fit range as well as consistent with the off-axis data. The results on $`\mathrm{\Delta }K_1`$ are displayed in Tab. 10. In Fig. 17 these shifts are plotted relative to $`K_1(0)`$, together with cubic splines. Since $`K_1(0)<0`$, at large $`ma`$ the curves will diverge towards the negative direction. Note that the slope of $`\mathrm{\Delta }K_1(ma)`$ at $`ma=0`$ has also been determined in Ref. Booth:2001qp  for SW fermions, however, by matching the vacuum polarization in momentum space, rather than the position space potential and $`\mathrm{\Delta }K_1`$ is not universal.

For the two $`𝒪(a)`$ improved actions $`\mathrm{\Delta }K_1(ma)`$ is small (a few %), relative to $`K_1(0)`$, at least for $`ma1/4`$. This is very different in the case of Wilson fermions, where even at $`ma=0.05`$ $`\mathrm{\Delta }K_1`$ corrects $`K_1(0)`$ by more than 50 %. Using the static position space potential as the prescription that defines $`\mathrm{\Delta }K_1`$, the $`ma`$ dependence is well parametrized by $`\mathrm{\Delta }K_1(x)=ax+bx^2`$ with $`a=1.20\pm 0.02`$ and $`b=5.7\pm 0.3`$ for $`x=ma0.1`$ in the Wilson case. For SW fermions the fit parameters read $`a=0.44\pm 0.08`$ and $`b=3.5\pm 0.8`$ while the KS data are compatible with zero within our accuracy over the whole range $`x0.25`$ and can be fitted by a straight line with $`a=0.04\pm 0.04`$.

### III.5 The “$`\beta `$-shift”

We have seen above that at distances $`rm^1`$ the running of the coupling is not affected by the presence of sea quarks anymore: at large distances, at least in perturbation theory, the effect of massive quarks can be integrated out into a shift of the coupling constant of the quenched theory. In contrast this is not possible for a theory with massless quarks, which completely decouples from the quenched case at all scales.

The theory with massive sea quarks and the quenched theory can be matched in the infra red by use of an intermediate mass-dependent scheme, such as the $`R`$-scheme, defined through the inter-quark potential in position space,

$`\alpha _R\left(r^1\right)`$ $`=`$ $`{\displaystyle \frac{r}{C_F}}V_c(r)`$ (46)
$`\stackrel{R=r/a\mathrm{}}{}`$ $`{\displaystyle \frac{Ra}{C_F}}V_{\text{int}}(R)`$ (47)
Inserting the perturbative expansion of $`V_{\text{int}}(R)`$ in terms of $`\alpha _L`$ one obtains,

$`\alpha _R(\mu )`$ $`=`$ $`\alpha _L^{(n_f)}+e_1^{(n_f)}(R)\alpha _{L}^{(n_f)}{}_{}{}^{2}+\mathrm{}`$ (48)
$`=`$ $`\alpha _L^{(0)}+e_1^{(0)}(R)\alpha _{L}^{(0)}{}_{}{}^{2}+\mathrm{},`$ (49)
where $`\mu =(Ra)^1`$, with known coefficients, $`e_1^{(n_f)}(R)`$. Note that these coefficients will depend on the function $`K_1(ma)`$ discussed above. The requirement of a unique (static source self-energy subtracted) potential $`V_{\text{int}}(R)`$ at $`R(ma)^1`$ then results in,

$$\alpha _L^{(0)}=\left\{1+\left[e_1^{(n_f)}(\mathrm{})e_1^{(0)}(\mathrm{})\right]\alpha \right\}\alpha _L^{(n_f)}+\mathrm{}.$$
(50)
One can now re-write the above equation in terms of the lattice parameter $`\beta =2N/g^2=N/(2\pi \alpha _L)`$. For $`n_f`$ degenerate quark flavours with mass $`m`$ and $`N=3`$ the result reads,

$`\beta ^{(n_f)}`$ $`=`$ $`\beta ^{(0)}+\mathrm{\Delta }\beta ,`$ (51)
$`\mathrm{\Delta }\beta `$ $`=`$ $`{\displaystyle \frac{n_f}{2\pi ^2}}\left[\mathrm{ln}(Dma)+3\pi \mathrm{\Delta }K_1(ma)\right],`$ (52)
with the numerical values for the constant,

$$D=\sqrt{C_0}\mathrm{exp}\left[3\pi K_1(0)\frac{5}{6}\right]:$$
(53)
$`D^W`$ $`=`$ $`0.448\pm 0.002,`$ (54)
$`D^{SW}`$ $`=`$ $`0.0238\pm 0.0001,`$ (55)
$`D^{KS}`$ $`=`$ $`0.726\pm 0.002,`$ (56)
for Wilson, SW and KS fermions, respectively. Since $`\mathrm{\Delta }\beta 0`$ as $`ma\mathrm{}`$, the $`\mathrm{ln}(Dma)`$ term has to be cancelled by $`3\pi \mathrm{\Delta }K_1(ma)`$ at large $`ma`$ in Eq. (52): the constant $`D`$ is identical with that appearing within Eq. (44). We discuss the matching procedure in some more detail in Appendix B.

Naïvely one would assume perturbation theory to be applicable as long as the relative $`\alpha `$-shift is small. To leading non-trivial order in perturbation theory the $`\beta `$-shift is linear in $`n_f`$. At the next order, the situation will be complicated by additional terms that are proportional to $`n_f^2`$ as can be seen from Eq. (B.2) of Appendix B.

### III.6 Boosted perturbation theory and $`q^{}`$

Lattice perturbation theory is well known for its bad convergence, partly due to large contributions from lattice tadpole diagrams. We recall that $`V_2(𝐑)`$ incorporates such a contribution: $`(2N^23)/(24N)V_1(𝐑)`$ \[cf. Eq. (III.1)\]. Hence reordering the series in terms of a better behaved expansion parameter like $`\alpha _V(q)`$, the coupling defined by the static QCD potential in momentum space,

$$\stackrel{~}{V}(q)=4\pi C_F\frac{\alpha _V(q)}{q^2},$$
(57)
is desirable in many cases. In some respect this is similar to the situation in continuum perturbation theory and resembles an expansion in terms of a renormalized, rather than a “bare”, lattice coupling parameter. We will refer to such techniques as “boosted perturbation theory”.

At $`𝒪(\alpha ^2)`$ we can write,

$$\alpha _V(q^{})=\alpha _L+\left[a_1+b_12\beta _0\mathrm{ln}(q^{}a)\right]\alpha ^2,$$
(58)
with $`a_1`$, $`b_1`$ and $`\beta _0`$ as defined in Eqs. (154), (156), Eqs. (134) – (A.2) and Eq. (122), respectively. The “optimal” scale $`q^{}`$ depends on the underlying process. It has been argued by Brodsky, Lepage and Mackenzie (BLM) Brodsky:1983gc ; lepage  that the logarithmic average of the momenta exchanged at tree-level is a particularly good choice for the scale $`q^{}`$ to be used within a one-loop calculation. The scale optimization procedure can also be generalized to higher order perturbative calculations Brodsky:1998mf ; Hornbostel:2002af . We illustrate the original recipe for the case of the position space potential. Here the tree-level calculation, Eq. (30), yields,

$$\mathrm{ln}\left[aq^{}(𝐑)\right]=\frac{d^3q\mathrm{ln}q\mathrm{sin}^2\left(\frac{\mathrm{𝐪𝐑}a}{2}\right)/\widehat{𝐪}^2}{d^3q\mathrm{sin}^2\left(\frac{\mathrm{𝐪𝐑}a}{2}\right)/\widehat{𝐪}^2}.$$
(59)
For $`R1`$ we find,

$$aq^{}(R)=f_1+\frac{f_2}{R}+f_3\frac{\mathrm{ln}R}{R}+\mathrm{},$$
(60)
with the numerical constants,

$$f_11.672,f_21.00,f_30.41.$$
(61)
The constants $`f_2`$ and $`f_3`$ have been obtained from a fit to $`R4`$ data while the $`r\mathrm{}`$ limit $`f_1`$ has been calculated directly. The fit curve and $`q^{}`$ values are displayed in Fig. 18. Note that violations of rotational invariance in $`q^{}(𝐑)`$ are remarkably small. Analogously we can obtain $`q^{}`$ values for the inter-quark force. We display results for the plaquette $`\mathrm{}`$, the chair, the “parallelogram”, the “rectangle” $`W(1,2)`$, the static self energy $`\delta m_{\text{stat}}=V_S/2`$, the potential and the force $`F_{R_1R_2}=a[V(R_2)V(R_1)]/(R_2R_1)`$ at selected distances in Tab. 11. In addition, the $`q^{}`$ values for the “tadpole improved” quantities $`m_{\text{stat}}+\mathrm{ln}\mathrm{}/(4a)`$ as well as $`V(1)+\mathrm{ln}\mathrm{}/(2a)`$ and $`V(2)+\mathrm{ln}\mathrm{}/(2a)`$ are included.

One can further convert between the $`V`$ and the $`\overline{MS}`$ scheme:

$$\alpha _{\overline{MS}}(\mu )=\alpha _V(q^{})\left[a_1+2\beta _0\mathrm{ln}(\mu /q^{})\right]\alpha ^2+\mathrm{}.$$
(62)
It has been argued Brodsky:1983gc  that the scale

$$\mu =e^{5/6}q^{},$$
(63)
at which the $`n_f`$ dependence of $`a_1`$ above is cancelled (in the massless case) was the optimal choice of scale for this conversion. Note that the ratio $`\mu /q^{}`$ above is independent of the number of colours $`N`$ or flavours $`n_f`$. The $`\mu `$ values are also included into the table.

The average plaquette turns out to be the most ultra violet quantity, with $`q^{}3.402/a`$, followed by the potential $`V(R=1)`$. Due to the self energy contribution $`V_S=V(\mathrm{})=2\delta m_{\text{stat}}`$, $`q^{}`$ is bounded from below by $`1.672/a`$ for the potential at all distances. However, in the case of the force $`q_F^{}0`$ as $`R\mathrm{}`$. As one would expect $`q^{}`$ also approaches zero at large distances for the potential if $`V_S`$ is subtracted. We find for instance $`q^{}a0.59`$ for $`V(1)V_S`$ as opposed to $`q^{}a2.86`$ for $`V(1)`$ alone.

### III.7 The static source self energy

We define,

$$V_S=\underset{R\mathrm{}}{lim}V(R)=2\delta m_{\text{stat}},$$
(64)
where $`\delta m_{\text{stat}}`$ is the lattice pole mass of a fundamental static colour source, e.g. a quark in the $`m\mathrm{}`$ limit. The above relation only holds in perturbation theory where the interaction energy vanishes as $`r\mathrm{}`$.

The static (or residual) mass $`\delta m_{\text{stat}}`$ has been calculated to $`𝒪(\alpha ^3)`$ for pure gauge theory DiRenzo:2001nd ; Trottier:2001vj  as well as to $`𝒪(\alpha ^2)`$ for massless Wilson and SW quarks Martinelli:1999vt . This has enabled a number of authors to obtain the $`b`$ quark mass $`\overline{m}_b(\overline{m_b})`$ from lattice simulations of static-light mesons. Since the sources are static the $`𝒪(\alpha )`$ result can be obtained from a tree-level calculation and we count this as leading order (LO), discarding the trivial $`𝒪(1)`$ value $`\delta m_{\text{stat}}=0`$. The one-loop $`𝒪(\alpha ^2)`$ result is then next to leading order (NLO) and the two-loop $`𝒪(\alpha ^3)`$ value NNLO. Note that the counting conventions employed in some of the literature differ from the one defined above, in which $`𝒪(\alpha ^n)`$ corresponds to an NmLO calculation with $`m=n1`$ for which diagrams involving $`m`$ loops have to be computed. For instance in Ref. Martinelli:1999vt  the value $`m=n`$ is used, creating the (wrong) impression that a two-loop calculation is required to obtain the $`𝒪(\alpha ^2)`$ result.

We have computed the static quark mass shift as a function of the mass of three species of sea quarks: KS, Wilson and SW. The results, most of which can also be read off from the last seven lines of Tab. 9, are displayed in Tab. 12. The labelling conventions are identical to those of Eqs. (32) – (36):

$`a\delta m_{\text{stat}}`$ $`=`$ $`\delta M_1g^2+\delta M_2g^4,`$ (65)
$`\delta M_2`$ $`=`$ $`\delta M_{pg}+C_Fn_f\delta M_f,`$ (66)
where

$`\delta M_1`$ $`=`$ $`{\displaystyle \frac{C_F}{2}}V_1(\mathrm{}),`$ (67)
$`\delta M_{pg}`$ $`=`$ $`{\displaystyle \frac{C_F}{2}}\{{\displaystyle \frac{N}{2}}[V_\mathrm{\Pi }(\mathrm{})+V_{NA}(\mathrm{})]`$ (68)
$`+`$ $`{\displaystyle \frac{2N^23}{24N}}V_T(\mathrm{})\}`$
and $`\delta M_f=\delta M_f^{(KS)}`$ for KS quarks and

$$\delta M_f=\delta M_f^{(0)}+c_{SW}\delta M_f^{(1)}+c_{SW}^2\delta M_f^{(2)}$$
(69)
for Wilson-SW fermions. We find the numerical values,

$`\delta M_1`$ $``$ $`0.126365504C_F`$ (70)
$`\delta M_{pg}`$ $``$ $`0.05292C_F,`$ (71)
where the latter number applies to $`SU(3)`$ gauge theory (cf. the last line of Tab. 6). Note that this result differs in the least significant digit from the value obtained in Ref. Martinelli:1999vt  that can be translated into $`\delta M_{pg}0.05297C_F`$. The mass dependence of $`\delta M_f`$ \[which in Eq. (66) is accompanied by a factor $`C_Fn_f`$\] is visualized in Fig. 19 for Wilson and KS quarks. $`\delta M_f`$ approaches zero like $`1/(ma)`$ as $`ma\mathrm{}`$.

Finally, in Fig. 20 we compare the results on $`\delta M_f(ma)`$ for quark masses $`ma0.25`$, normalized to $`\delta M_f(0)`$, from all three quark actions with each other. While in absolute terms $`\delta M_f`$ for SW quarks turns out to be much bigger than for the other two actions (cf. Tab. 12) we find the relative variation with the quark mass to be much weaker for both $`𝒪(a)`$ improved actions than for Wilson fermions. This is consistent with our observations for small Wilson loops and $`\mathrm{\Delta }K_1`$ above. For $`ma0.1`$ $`\delta M_f(ma)`$ is well parametrized by a quadratic function,

$$\delta M_f(ma)=\delta M_f(0)+bma+c(ma)^2,$$
(72)
with $`\delta M_f(0)=0.0011667(3),0.0013404(2)`$ and $`0.0030351(11)`$, $`b=0.00101(8)`$, 0.00350(1) and 0.00245(5) and $`c=0.0014(8)`$, $`0.00543(5)`$ and $`0.0020(4)`$ for KS, Wilson and SW quarks, respectively.

Starting from the expansion of the plaquette,

$$\mathrm{}=1c_1\alpha _Lc_2\alpha ^2+\mathrm{},$$
(73)
we can rearrange the above expansion of the self energy in the following way,

$`a\delta m_{\text{stat}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}\mathrm{}}{4}}+m_1^t\alpha _L+m_2^t\alpha _L^2`$ (74)
$`=`$ $`{\displaystyle \frac{\mathrm{ln}\mathrm{}}{4}}+m_1^t\alpha _{\overline{MS}}(\mu )+\overline{m}_2^t\alpha ^2,`$ (75)
where $`\mu =q^{}e^{5/6}=0.363/a`$, $`m_1^t=(4\pi )\delta M_1c_1/4`$, $`m_2^t=(4\pi )^2\delta M_2c_2/4c_1^2/8`$. We will refer to practices such as adding the $`\mathrm{ln}\mathrm{}`$ term non-perturbatively and subtracting it perturbatively as “tadpole improvement” of the observable while in the second step, Eq. (75) we “boost” the perturbation theory. Here we choose to convert everything into the $`\overline{MS}`$ scheme. We prefer this to plaquette based schemes or the $`V`$ scheme as now all mass dependence is made explicit in the coefficient function $`\overline{m}_2^t(ma)`$. A conversion of the result into another scheme of choice including mass dependent schemes is easily possible with the help of Appendix A. We will obtain the $`\overline{MS}`$ coupling from the average plaquette and other quantities in Sec. IV.2 below.

For $`SU(3)`$ gauge theory with $`n_f`$ quark flavours of masses $`ma0.1`$ we find,

$`m_1^t`$ $`=`$ $`1.0700768,`$ (76)
$`m_2^t`$ $`=`$ $`7.6104+n_f\left[a_f+b_fma+c_f(ma)^2\right],`$ (77)
$`\overline{m}_2^t`$ $`=`$ $`0.5839+n_f\left[\overline{a}_f+\overline{b}_fma+\overline{c}_f(ma)^2\right],`$ (78)
with the constants $`a_f`$, $`b_f`$, $`c_f,\mathrm{}`$ of Tab. 13. Tadpole improvement reduces the pure gauge value $`\delta m_{pg}11.14`$ to $`m_{pg}^t7.61`$ and boosting reduces this further down to $`\overline{m}_{pg}^t0.584`$. The fermionic coefficients $`\delta m_fa_f\overline{a}_f`$ also undergo a reduction in each step. The NNLO coefficient $`\delta m_386m_3^t67.5\overline{m}_3^t2.85`$ is only known for $`n_f=0`$.

When applying the above result to extract say $`m_b`$ from lattice simulations in the static limit one can readily ignore the mass dependent terms for Wilson fermions, which are $`𝒪(a)`$ effects. However, it does not harm to include them either. In the cases of massive KS and SW sea quarks which are $`𝒪(a)`$ improved at least the $`b_f`$s have to be included.

## IV Comparison with non-perturbative data

We shall compare and apply our perturbative calculations to non-perturbative results obtained in lattice simulations. For this purpose we will use several data sets that have been obtained by different collaborations. All quenched reference data are from simulations of one of the co-authors of this article and collaborators (Ref. Bali:2000gf  and references therein) while $`n_f=2`$ data have been provided by the SESAM/T$`\chi `$L Collaboration Bali:2000vr ; Bali:2002  (Wilson fermions), the UKQCD Collaboration Allton:2001sk ; Booth:2001qp  (SW fermions) and the MILC Collaboration Tamhankar:2000ce  (KS fermions444While KS quarks are only defined for multiples of four mass-degenerate fermion flavours, many authors have attempted to emulate $`n_f=2`$ (or even $`n_f=3`$) by using the (positive) square root of the $`n_f=4`$ fermionic determinant in their simulation. The MILC Collaboration adopt this strategy. However, it is not clear that the resulting lattice action corresponds to a local field theory like QCD. In our analysis we shall set this problem aside. We remark that the perturbation theory generated by the use of this action indeed corresponds to replacing $`n_f=4`$ by $`n_f=2`$ to all orders, at least as long as no external fermion lines are encountered. The comparison with simulation results, however, might be meaningless should no universal continuum limit exist.).

We will apply our perturbative results to the “$`\beta `$-shift” encountered when including sea quarks — relative to the quenched approximation. Subsequently we determine the QCD running coupling “constant” for $`n_f=0`$ and $`n_f=2`$ from lattice data of the correlation length Sommer:1994ce  $`R_0=r_0/a`$, the average plaquette and the short distance static potential. In this context we will also make use of the CP-PACS Yoshie:2000wd  ensemble obtained with $`n_f=2`$ SW fermions and the Iwasaki gluonic action Iwasaki:1984cj . Finally we shall also compare perturbative and non-perturbative lattice potentials with the aim of parametrizing lattice artefacts and to resolve the differences in the running of the Coulomb couplings between quenched and un-quenched data sets.

### IV.1 The non-perturbative “$`\beta `$-shift”

We study the situation of $`n_f`$ mass-degenerate flavours of sea quarks. Two parameters can be varied: the coupling $`\beta =6/g^2`$ and the lattice quark mass $`ma`$ which in the case of Wilson-SW fermions is related to the parameter $`\kappa `$. Each $`\kappa `$-$`\beta `$ (or $`ma`$-$`\beta `$) combination can be translated into a pion mass $`m_\pi r_0`$ and a lattice spacing $`a/r_0`$ where $`r_0`$ is a correlation length defined below. As $`\beta \mathrm{}`$ we reach the continuum limit, $`a/r_00`$. Sending the quark mass $`mr_0`$ to zero corresponds to a vanishing pion mass, $`m_\pi r_00`$ (at finite $`\beta `$: only up to violations of chiral symmetry). However, in general the two limits will “mix”: varying $`\beta `$ at fixed $`ma`$ will not only affect $`a/r_0`$ but also to some extent the ratio $`m_\pi r_0`$ while a variation of $`ma`$, keeping the coupling fixed, will result in a change of $`a/r_0`$ as well.

In view of the computational cost of lattice simulations incorporating sea quarks, predicting by what amount the coupling has to be shifted in order to compensate for the change in $`a/r_0`$ that is induced by varying the quark mass $`ma`$ is certainly desirable. We will explore this possibility by comparing the expectation Eq. (52) on the perturbative $`\beta `$-shift against results from lattice simulations with $`n_f=2`$ sea quarks of masses $`0.0098ma0.3`$.

At large distances the potentials will be dominated by non-perturbative effects: the quenched potential will linearly rise ad infinitum. However, as soon as sea quarks are introduced the $`Z_3`$ symmetry of the action is broken and at some distance Bali:2000vr  $`r_c(2.10+1.56mr_0)r_0`$ string breaking will set in; obviously the linear quenched behaviour at large distances can never be emulated by a theory with a completely flat potential at large $`r>r_c`$. Sea quarks will also affect the running of the coupling at short distances $`r<0.3m^1`$. Hence, the best we can hope for is that sea quarks decouple within a window of distances $`m^1r<r_c`$.

We shall define a

$$\mathrm{\Delta }\beta _{r_0}(n_f,r_0/a,ma)=\beta ^{(n_f)}(r_0/a,ma)\beta ^{(0)}(r_0/a),$$
(79)
by non-perturbatively matching un-quenched and quenched $`\beta `$ values that correspond to the same scale Sommer:1994ce  $`r_00.5`$ fm, implicitly defined through the relation,

$$r_0^2\frac{dV(r)}{r}|_{r=r_0}=1.65.$$
(80)
If non-perturbative effects cancel from Eq. (79) and leading order perturbation theory applies then $`\mathrm{\Delta }\beta _{r_0}`$ will only depend on $`ma`$ (and $`n_f`$) but not on the lattice spacing $`a/r_0`$. Perturbation theory will obviously break down for $`ma1`$ in which case the coefficients of the perturbative expansion explode \[cf. Eq. (52)\]. This is a reflection of the fact that the running of the coupling in the theory with massless sea quarks differs from the $`n_f=0`$ case at all scales.

On the other hand, unless $`ma1/D`$ the relative importance of the non-universal term $`\mathrm{\Delta }K_1(ma)`$ will increase and in general the matching of the running of the inter-quark force at a scale $`r_0`$ is not equivalent to the matching of the running of the potential at larger distances anymore. In principle one can work out the matching condition for $`r_0`$ in perturbation theory. It turns out, however, that at distances $`r_00.5`$ fm the non-perturbative contribution to the lattice potential is already substantial. For instance a perturbative $`n_f=0`$ calculation yields $`r_0/a=5.33`$ at $`\beta 1.49`$ while non-perturbatively this result is obtained at $`\beta 6.0`$ and $`\beta `$ depends only logarithmically on $`a`$!

Given the above situation, it is difficult to identify a sensible way of combining perturbative and non-perturbative results at large distances. Consequently, we refrain from attempting to do this but separately employ purely perturbative and purely non-perturbative matching strategies: on the lattice simulation side of things we match $`r_0/a`$. In perturbation theory we match a quantity that is suitable for perturbative treatment, namely the inter-quark potential at large distances, rather than a quantity that is inspired by (non-perturbative) phenomenology, like $`r_0`$. We then hope that some, ideally most, non-perturbative effects cancel each other on the right hand side of Eq. (79) and leave the $`\beta `$-shift within the available window of $`r_0/a`$ and $`ma`$ untouched. Clearly, the perturbative matching will break down if either $`ma`$ or $`\beta `$ are too small. In addition the strategy that we adopt requires $`ma<1/D`$ to limit violations of universality as well as $`mr_0>0.3`$ to ensure that the sea quarks do not affect the running of the coupling at distances around $`r_0`$.

The quenched $`\beta `$ value that corresponds to a given $`r_0/a`$ is determined by use of the interpolation,

$$r_0/a=\mathrm{exp}\left(d_0+d_1x+d_2x^2+d_3x^3\right),x=\beta 5.9,$$
(81)
with $`d_01.489,d_11.982,d_20.630,d_31.522`$, obtained from a fit to lattice results within the window Bali:2000gf  $`5.5\beta 6.2`$. In principle we can also use the more recent precision data of Refs. Guagnelli:1998ud ; Necco:2002xg . However, the difference is insignificant in view of the present level of accuracy of un-quenched data and the matching of some of the MILC data requires an interpolation that extends to lattice spacings coarser than those investigated in the latter two references. In the un-quenched simulations, different groups used different procedures to extract $`r_0/a`$ which partly explains why different collaborations can obtain very different error estimates with similar computational efforts. Aiming only at a qualitative comparison we shall not attempt to reanalyse all data in one and the same way but prefer to cite the published values and errors.

In the case of KS fermions Tamhankar:2000ce ; Heller:1994rz  there is no additive mass renormalization. For Wilson fermions Bali:2000vr ; Bali:2002  we define the lattice quark mass via,

$$ma=\frac{1}{2}\left(\frac{1}{\kappa }\frac{1}{\kappa _c}\right),$$
(82)
where $`\kappa _c`$ has been obtained from a chiral extrapolation of the non-perturbative $`m_\pi ^2`$ at fixed555 In principle other extrapolations are possible like keeping $`r_0/a`$ fixed Allton:2001sk . However, the whole matching idea is based on a semi-quenched philosophy: what value of the $`n_f=0`$ coupling will produce the same infra red physics, e.g. $`r_0/a`$? Since at $`\kappa _c`$ the low energy physics will be very different anyway, the “natural” choice in this case seems to be keeping the coupling constant fixed. Having said this, to the order of perturbation theory at which we work both approaches are equivalent anyway. $`\beta `$666To all orders in perturbation theory $`ma`$ is defined via Eq. (82). To the order that we work at, $`\kappa _{c,\text{pert}}=1/8`$, which means that the $`\kappa `$ values employed in the lattice simulation all correspond to negative masses, $`ma`$, when directly plugged into the perturbative expansion. Obviously this is not a sensible choice, and hence we formulate our perturbation theory in terms of the mass $`ma`$ rather than $`\kappa `$. Subsequently, we determine the $`\kappa `$ that corresponds to a given $`ma`$ value non-perturbatively. In the case of SW fermions we proceed in an analogous way: on the perturbation theory side we use $`ma`$ and $`c_{SW}=1`$, which is the value consistent with the order at which we work. However, in the simulation we use the $`\kappa `$ and $`c_{SW}`$ values respectively which result in the same $`ma`$ value and that eliminate $`𝒪(a)`$ lattice artefacts non-perturbatively.. In the case of the SW data Allton:2001sk ; Booth:2001qp  we do not know $`\kappa _c`$ but we can use the published values of the PCAC quark mass Allton:2001sk , which is equivalent to any other definition to leading order perturbation theory.

We compile results from the three collaborations in Tabs. 1416. We have sorted Tab. 16 (with the exception of $`\beta =5.6`$) with respect to the quark mass and then $`\beta `$ while the other tables are sorted the other way around. Unfortunately, no PCAC mass value is available for the most critical UKQCD data set \[$`(\beta ,\kappa )=(5.2,0.13565)`$\], however, we arrived at the estimate $`ma=0.0179(9)`$ by extrapolating $`ma`$ as a linear function of $`(m_\pi r_0)^2`$ to the value Hepburn:2002fh  $`(m_\pi r_0)^2=0.132(2)`$. The perturbative prediction Eq. (52) is included in the last column of the tables where the error is due to the uncertainty in $`\mathrm{\Delta }K_1(ma)`$.

In Figs. 2123 we compare the non-perturbative results with the expectations. In all cases $`ma`$ is much smaller than $`1/D`$ as desired. However, $`mr_0`$, varying from 0.06 to 0.24 in the Wilson case, from 0.09 to 0.46 in the SW case and from 0.04 to 0.6 in the KS case, is not exactly big relative to $`r_0=0.3m^1`$ (the distance below which the perturbative running of the coupling changes and sea quark effects cannot completely be compensated for by a scale redefinition alone). Nonetheless, the first two figures reveal an excellent agreement with the expectation within 10 % for Wilson-type quarks, even for such light quark masses. Moreover, no significant $`\beta `$-dependence is observed, again in agreement with Eq. (52), indicating higher order corrections to be small. This is very different for the KS fermions depicted in Fig. 23: here the agreement with the prediction only improves as the lattice spacing is reduced.

In Fig. 24 we finally compare the three formulations within the window $`0.1ma0.009`$ for $`4<r_0/a<6.5`$ data. The Wilson and SW results seem to fall onto almost universal curves that differ by less than 10 % from the prediction while the KS results deviate by much more and some $`\beta `$-dependence is evident, Whether this is due to a slower convergence of the perturbative series, due to the inexact updating algorithm employed or due to $`n_f`$ not being a multiple of four is an open question. Around $`m0.05r_0^120`$ MeV the reliability of the matching for Wilson-SW quarks finally appears to break down (left-most data point in Fig. 24): the behaviour becomes “truly un-quenched”.

We remind the reader that $`\beta `$-shifts are in general not independent of the quantity that is used to match quenched and un-quenched theories. This ambiguity exists non-perturbatively as well as in perturbation theory. The qualitative agreement between prediction and simulation for Wilson-type quarks indicates that, at least at present masses, physics at hadronic scales is not yet strongly affected by quark loops, which is consistent with the phenomenological success of the quenched approximation.

In simulations with non-perturbatively improved SW quarks the lines of constant $`a/r_0`$ and constant $`m_\pi r_0`$ are significantly tilted with respect to both axes of the $`\beta `$-$`ma`$ plane Irving:2001hs . This observation is consistent with the small value, $`D^{SW}0.0238`$, but causes practical problems as going to lighter quark masses at sensible $`a/r_0`$ values requires simulations at small $`\beta `$s for which the non-perturbative determination of the improvement coefficient $`c_{SW}`$ causes problems Jansen:1998mx . In the worst case it might even be conceivable that the slope of the variation of $`a/r_0`$ with $`\kappa \kappa _c`$ at fixed $`\beta `$ eventually diverges and that the continuum limit $`\beta \mathrm{}`$ does not at all exist for light quarks Aoki:1983qi ; Aoki:2001xq . The variation of $`r_0`$ with the quark mass, however, is reduced when actions with only moderately negative values of $`K_1(0)`$ \[and therefore $`D`$ values that are of order one\] like the Wilson or KS action are employed.

Our result might be taken as an indication that the perturbative matching of $`n_f1`$ to $`n_f`$ QCD couplings in the $`\overline{MS}`$ scheme via the intermediate mass-dependent $`V`$-scheme Brodsky:1998mf  (or by calculating and matching a physical amplitude Bernreuther:1983zp ) is likely to be quite reliable, even at masses as light as or lighter than that of the charm quark. In this case the matching condition reads,

$$\frac{1}{\alpha _{\overline{MS}}^{(n_f)}(\mu )}=\frac{1}{\alpha _{\overline{MS}}^{(n_f1)}(\mu )}+\frac{1}{3\pi }\left[\mathrm{ln}\left(\sqrt{C_0}\frac{m}{\mu }\right)\frac{5}{6}\right]\alpha ,$$
(83)
which may be rewritten as,

$$\alpha _{\overline{MS}}^{(n_f)}(\mu )=\alpha _{\overline{MS}}^{(n_f1)}(\mu ),$$
(84)
where,

$$\mu =\sqrt{C_0}\mathrm{exp}\left(\frac{5}{6}\right)m=0.990(3)m,$$
(85)
to this order in perturbation theory777 Eqs. (84) – (85) can also be cast into $`\alpha ^{(n_f)}(m)=\alpha ^{(n_f1)}(m)+c\alpha ^2`$ where the numerical value $`c0.0033`$ is of the same magnitude as Bernreuther:1983zp  $`c=7/(72\pi ^2)0.0099`$. Note that the difference between pole and $`\overline{MS}`$ masses is irrelevant to this order in perturbation theory.. Clearly, the applicability of the matching method in this case is ultimately limited by the reliability of the perturbative running of the $`\overline{MS}`$ coupling at small scales.

### IV.2 Determining $`\alpha _{\overline{MS}}`$

Lattice simulations yield correlation lengths and masses that are functions of a set of input parameters, namely the inverse lattice coupling $`\beta `$ and lattice quark masses $`ma`$ (or in the case of Wilson-type fermions $`\kappa `$-values that can be related to the quark masses). One can then for instance obtain the strong coupling constant in the $`\overline{MS}`$ scheme,

$$\alpha _{\overline{MS}}(a^1)=\alpha _L+b_1\alpha _L^2+\mathrm{},$$
(86)
once a physical scale has been assigned to the lattice spacing $`a`$. One such possibility would be to “measure” $`R_0=r_0/a`$ on the lattice and to equate $`r_00.5`$ fm. Other input scales with a more direct connection to experiment are possible, for instance the proton mass $`m_N`$ or the pion decay constant $`f_\pi `$. If QCD with the right number of quark flavours and masses is simulated the resulting $`a`$ should become independent of the choice of the experimental input quantity. In this way the strong coupling constant can be determined from low energy hadron phenomenology.

There are of course higher order perturbative as well as non-perturbative corrections to Eq. (86) which, however, will vanish as $`a0`$. In practice these corrections are big in the range of lattice spacings that will ever be realistically accessible Bali:1993ru ; lepage . One way to improve the convergence is to convert between the $`\overline{MS}`$ and the lattice schemes at an optimized BLM scale Brodsky:1983gc ; lepage ; Brodsky:1998mf  $`\mu =e^{5/6}q^{}`$, rather than at $`a^1`$. We refer to this reorganization of the perturbative series as boosted perturbation theory. Another possibility (that can be combined with the BLM scheme) is to “measure” the coupling on the lattice. This can either be done from a quantity that depends on the lattice spacing like the average plaquette Parisi:1980pe ; lepage ; El-Khadra:1992vn ; davi1 ; davi2  or at a scale $`\mu a^1`$ from quantities that have a well defined continuum limit Booth:1992bm ; Bali:1993ru ; Necco:2001gh ; Garden:2000fg . Here we shall follow the former strategy.

#### IV.2.1 The method

We collect the average plaquette, the potential at $`R=1`$ and at $`R=2`$ obtained in quenched simulations as well as from simulations with sea quarks in Tabs. 1721. In addition we include the “force”,

$`F_{12}`$ $`=`$ $`a\left[V(2)V(1)\right],`$ (87)
$`F_{23}`$ $`=`$ $`a\left[V(3)V(2)\right].`$ (88)
We will restrict our discussion to a one-loop determination of the running coupling. We shall also calculate two-loop corrections for the pure gauge case. Two-loop results are also known for the case of the plaquette with massive Wilson quarks888Unfortunately, these results are of limited use since they have been obtained at $`\kappa ^1`$ values that correspond to negative quark masses, after subtracting $`\kappa _c^1`$ to the same order in perturbation theory. Alles:1998is .

Obviously many ways, some better than others, to extract $`\alpha _{\overline{MS}}`$ exist that are consistent with perturbation theory at a given order. For convenience we adopt the procedure detailed in Refs. davi1 ; davi2  but point out that many alternative ways are equally justified (e.g. that of Ref. Booth:2001qp ). We only differ from Refs. davi1 ; davi2  in as far as, once the coupling has been converted into the $`\overline{MS}`$ scheme, we evolve the scale and extract the $`\mathrm{\Lambda }`$-parameter by numerically integrating the four-loop $`\beta `$ function, rather than using a perturbatively truncated formula.

We can write the plaquette as,

$$\mathrm{}=1c_1\alpha _Lc_2\alpha _L^2c_3\alpha _L^3,$$
(89)
where $`c_1`$ and $`c_2`$ can be read off from Tab. 1 with the help of Eqs. (19) – (22) for the pure gauge case and the fermionic contributions can be identified in Tabs. 4 and 5. In particular one finds $`c_1=4\pi C_F/4`$. $`c_3`$ for the pure gauge case and sources in the fundamental representation reads Alles:1994dn ; Alles:1998is :

$`c_3`$ $`{\displaystyle \frac{(N^21)N}{8}}[0.0063538{\displaystyle \frac{0.0181239}{N^2}}`$ (90)
$`+{\displaystyle \frac{0.0185221}{N^4}}](4\pi )^3.`$
We now define the quantity,

$`P`$ $`=`$ $`{\displaystyle \frac{\mathrm{ln}\mathrm{}}{c_1}}`$ (91)
$`=`$ $`\alpha _L+\stackrel{~}{c}_1\alpha _L^2+\stackrel{~}{c}_2\alpha _L^3,`$ (92)
$`\stackrel{~}{c}_1`$ $`=`$ $`{\displaystyle \frac{c_2}{c_1}}+{\displaystyle \frac{c_1}{2}},`$ (93)
$`\stackrel{~}{c}_2`$ $`=`$ $`{\displaystyle \frac{c_3}{c_1}}+c_2+{\displaystyle \frac{c_1^2}{3}}.`$ (94)
The expansion of $`P`$ in terms of $`\alpha _L`$ suffers from big coefficients. For instance one obtains $`\stackrel{~}{c}_13.37`$, $`\stackrel{~}{c}_217.69`$ in $`SU(3)`$ pure gauge theory. The origin of these big numbers can be traced to the lattice tadpole diagrams lepage . Large coefficients are also encountered when $`\alpha _L`$ is converted into the more “physical” coupling $`\alpha _V`$ \[or $`\alpha _{\overline{MS}}`$ which is “close” to $`\alpha _V`$, cf. Eq. (175)\]. By re-expressing the series Eq. (92) in terms of $`\alpha _V`$, taken at a suitable scale $`q^{}`$ Brodsky:1983gc , one might hope that the two effects cancel in part lepage :

$$P=\alpha _V(q^{})+X_1(q^{}a)\alpha _V^2(q^{})+X_2(q^{}a)\alpha _V^3+\mathrm{}.$$
(95)
We obtain,

$`X_1(x)`$ $`=`$ $`\stackrel{~}{c}_1+B_1(x)`$ (96)
$`X_2(x)`$ $`=`$ $`\stackrel{~}{c}_2+B_2(x)+2\stackrel{~}{c}_1B_1(x),`$ (97)
where

$`B_1(x)`$ $`=`$ $`b_1a_1+2\beta _0\mathrm{ln}x`$ (98)
$`B_2(x)`$ $`=`$ $`b_2a_2+2(a_1b_1+b_1^2+a_1^2)`$ (99)
$`+`$ $`2\beta _1\mathrm{ln}x+B_1^2(x)B_1^2(1).`$
The constant $`b_1`$ is defined in Eqs. (134) – (137), $`b_2`$ in Eqs. (140) – (A.2) and $`\beta _i`$ in Eqs. (122) – (123). The $`a_i`$ are defined in Eqs. (154) and (A.3) and have to be taken at the scale $`q^{}`$ for massive sea quarks. The numerical values for $`n_f=0`$ and $`q^{}a3.402`$ read: $`X_11.191`$, $`X_21.688`$. We can now truncate Eq. (95) at $`𝒪(\alpha ^2)`$ and obtain,

$`\alpha _V(q^{})`$ $`=`$ $`\alpha _PX_2\alpha _P^3+\mathrm{},`$ (100)
$`\alpha _P`$ $`=`$ $`{\displaystyle \frac{1}{2X_1}}\left(\sqrt{4PX_1+1}1\right).`$ (101)
Note that the coefficients $`X_i`$ inherit a mass dependence from $`a_i`$, $`b_i`$ (through $`\mathrm{\Delta }K_i`$) and $`\stackrel{~}{c}_i`$. With

$`\alpha _{\overline{MS}}(\mu )`$ $`=`$ $`\alpha _V(q^{})+Y_1\alpha _V^2(q^{})+Z_2\alpha _V^3(q^{}),`$ (102)
$`\mu `$ $`=`$ $`e^{5/6}q^{}1.478/a,`$ (103)
$`Y_1`$ $`=`$ $`\left(a_1{\displaystyle \frac{5}{3}}\beta _0\right),`$ (104)
$`Z_2`$ $`=`$ $`\left[a_2a_1^2{\displaystyle \frac{5}{3}}\beta _1\left(a_1{\displaystyle \frac{5}{3}}\beta _0\right)^2\right],`$ (105)
we arrive at,

$$\alpha _{\overline{MS}}(\mu )=\alpha _P+Y_1\alpha _P^2+Y_2\alpha _P^3,$$
(106)
where

$`Y_1`$ $`=`$ $`{\displaystyle \frac{2N}{3\pi }}{\displaystyle \frac{n_f}{6\pi }}\mathrm{ln}\left(1+{\displaystyle \frac{C_0m^2}{\mu ^2}}\right),`$ (107)
$`Y_2`$ $`=`$ $`Z_2X_2.`$ (108)
For $`m=0`$ $`Y_1`$ is independent999 For $`n_f`$ degenerate massive flavours one could in principle maintain the mass independence of this coefficient: $`Y_1Y_111N/(12\pi )\mathrm{ln}(1+C_0m^2/\mu ^2)`$, $`\mu \mu \times \sqrt{1+C_0m^2/\mu ^2}`$. In the most interesting case, QCD with non-degenerate flavours, this is however not possible. of $`n_f`$. For $`n_f=0`$, $`N=3`$ we find the numerical value, $`Y_20.9538`$: the NNLO correction is small.

In analogy to defining a coupling from the measured plaquette, other couplings can be computed from force and potential. We illustrate this procedure at NLO: in a first step we can write,

$`{\displaystyle \frac{aV(𝐑a)}{v_1(𝐑)}}`$ $`=`$ $`\alpha _L+{\displaystyle \frac{v_2(𝐑)}{v_1(𝐑)}}\alpha ^2`$ (109)
$`=`$ $`\alpha _V(q^{})+\left[{\displaystyle \frac{v_2(𝐑)}{v_1(𝐑)}}+B_1(q^{}a)\right]\alpha ^2,`$
where the function $`B_1(x)`$ is defined in Eq. (98). A similar expression can easily be written down for the force. The respective $`q^{}a`$ values can be found in Tab. 11. Consequently, $`\alpha _V(q^{})`$ can be obtained by solving the quadratic equation Eq. (109) and converted into $`\alpha _{\overline{MS}}(\mu )`$ via Eqs. (106) and (107), where we set $`Y_2=0`$ in our $`𝒪(\alpha ^2)`$ calculation.

Finally, we can run $`\alpha _{\overline{MS}}`$ numerically to arbitrarily high scales using the perturbative four-loop $`\beta `$ function vanRitbergen:1997va ; Tarasov:1980au  and then determine the $`\mathrm{\Lambda }_{\overline{MS}}`$-parameter, defined as in Eq. (125).

#### IV.2.2 Pure gauge theory

We display one- and two-loop results on $`\alpha _{\overline{MS}}(\mu )`$ as obtained from the logarithm of the average plaquette following the boosted perturbation theory procedure detailed above (as well as in Refs. davi1 ; davi2 ), in Tab. 22. This is then numerically converted into estimates of the QCD $`\mathrm{\Lambda }_{\overline{MS}}`$-parameter which are displayed in the last two columns of the table. In Fig. 25 we compare various methods of determining the $`\overline{MS}`$ $`\mathrm{\Lambda }`$-parameter as a function of the lattice spacing. The horizontal error band corresponds to the continuum limit result as obtained by the ALPHA Collaboration Garden:2000fg  by use of finite volume techniques. $`\alpha _L`$ refers to a conversion from the bare lattice coupling $`\alpha _L=3/[2\pi \beta (a)]`$ and $`\alpha _P`$ to a coupling, “measured” from the average plaquette. NLO and NNLO refer to results from a naïve perturbative conversion of $`\alpha _L`$ into $`\alpha _{\overline{MS}}(a^1)`$ at order $`\alpha ^2`$ and $`\alpha ^3`$ respectively while the abbreviations bNLO and bNNLO correspond to the boosted perturbation theory procedure detailed above.

We find that the ratios between NLO and NNLO results can be brought closer to unity, both by “measuring” the coupling and by “boosting” the perturbation theory. The combination of both methods indeed brings the result in agreement with the known number (error band). Interestingly, almost all bNLO and bNNLO results are within the expected range, even at rather coarse lattice spacings.

In Tab. 23 we display values $`\alpha _{\overline{MS}}(\mu )`$ obtained from the short distance lattice potential and force101010 One could also imagine to repeat this procedure for a “tadpole improved” potential, $`V(R)+\mathrm{ln}\mathrm{}/(2a)`$ at the corresponding $`q^{}`$s. Since in this case the result is obviously correlated with that which we obtained from the plaquette we leave this exercise to the interested reader and remark that any such additive corrections cancel in the force.. These are then converted into estimates of $`\mathrm{\Lambda }`$-parameters and compared to the corresponding bNLO estimates from the average plaquette in Tab. 23. We also display these numbers versus the inverse momentum scale $`1/\mu ^2`$ in Fig. 26. The curves represent results from phenomenological111111Since $`\mathrm{\Lambda }`$ cannot directly be obtained from position space Greens functions there is no theoretically well founded reason to assume the leading order lattice corrections to be quadratic in the lattice spacing. In fact we know that perturbative corrections to $`\mathrm{\Lambda }`$ exist which should be of order $`\mathrm{\Lambda }/\mathrm{ln}[(a\mathrm{\Lambda })^2]`$. quadratic plus quartic fits in $`1/(\mu r_0)`$. The $`\mathrm{\Lambda }_{\overline{MS}}r_0`$ values resulting from these continuum limit extrapolations are displayed in the last row of the table. The results from the high $`q^{}`$ quantities $`V(1)`$ and $`\mathrm{}`$ turn out to be in reasonable agreement with that from the ALPHA Collaboration, $`\mathrm{\Lambda }_{\overline{MS}}r_0=0.602(48)`$, however, this is not the case for $`F_{12}`$ or $`V(2)`$.

Unfortunately, there exists no first principles way to estimate systematic uncertainties in any determination of $`\mathrm{\Lambda }_{\overline{MS}}`$ that is partially based on perturbation theory at momenta as low as a few GeV. For instance the $`F_{12}`$ data exhibit a plateau for $`\beta 5.9`$ that is in no way inferior to that obtained from the average plaquette. However, the two extrapolated values differ by almost three standard deviations of the ALPHA Collaboration result. This suggests that the good agreement obtained for the plaquette might be partly accidental. We are unaware of any convincing argument why an estimate of $`\alpha _{\overline{MS}}(\mu )`$ extracted from $`F_{12}`$ at $`\beta =6.6`$, i.e. $`\mu r_05.65`$ should be less reliable than the value obtained from the plaquette at say $`\beta =5.8`$ ($`\mu r_05.4`$). The same can be said about the combination $`aV(2)`$ at $`\beta =6.6`$, $`\mu r_012.8`$ and the plaquette at $`\beta =6.3`$, $`\mu r_012.4`$ and yet the former two $`\beta =6.6`$ estimates $`\mathrm{\Lambda }_{\overline{MS}}^{F_{12}}r_0=0.717(13)`$ and $`\mathrm{\Lambda }_{\overline{MS}}^{V_2}r_0=0.764(24)`$ happen to lie significantly higher than the known continuum limit result.

We conclude that without the comparison between predictions from different non-perturbatively obtained observables we might easily have underestimated the systematic uncertainties of the approach. Having different observables as well as the ALPHA result at hand it becomes obvious that the 4 % difference between bNLO and bNNLO estimates from the plaquette do not necessarily reflect the whole truth, but that a conservative approach would instead argue in favour of a 15 % error on the $`\mathrm{\Lambda }`$-parameter in an $`𝒪(\alpha ^2)`$ determination like the above. Quantities with higher $`q^{}`$ values are certainly likely to behave better, but incorporating a number of measurements to estimate systematic errors is prudent when such a scale setting prescription has been used.

The present article contains the NLO perturbation theory that is necessary for the use of additional short distance quantities such as $`V(1)`$ and the chair and parallelogram Wilson loops. The observed scattering between extrapolations from different input quantities will hopefully shrink both as finer lattices are simulated, and once higher order corrections are known for potential and force. In the meantime we quote the value,

$$\mathrm{\Lambda }_{\overline{MS}}^{(0)}r_0=0.609(4)(90),$$
(110)
where the second error is systematic.

#### IV.2.3 Results with sea quarks

We shall limit our discussion of the fermionic case to $`𝒪(\alpha ^2)`$ since only the plaquette is known to NNLO and only for Wilson quarks at some negative mass values Alles:1998is . The fermionic correction to $`c_2`$ within Eqs. (89) – (92) can be calculated from Tab. 5 with the help of Eqs. (19) and the definitions Eqs. (14) – (16). The corrections to the potential are available in Tabs. 89, with the conventions of Eqs. (III.1) and (36) and $`v_i=(4\pi )^iV_i`$. We also include data that has been obtained by the CP-PACS Collaboration Yoshie:2000wd  by combining SW fermions with the Iwasaki gluon action into the analysis. The corresponding perturbative expansion of the plaquette can be read off from Tabs. 3 and 4, with the definitions Eqs. (19) – (21). In this case $`b_11.2466+n_fK_1^{SW}(ma)`$ Aoki:1998ar  and $`K_1`$ is independent of the gauge action. We find $`q^{}a3.213`$ for the plaquette. We display CP-PACS simulation parameters and results on $`r_0/a`$ and the plaquette in Tab. 25. The quark masses are obtained via Eq. (82). We set $`c_{SW}=1`$ in the analysis of SW-Wilson and SW-Iwasaki (SW-I) results which is consistent at NLO.

We display our results on $`\alpha _{\overline{MS}}(\mu )`$ as well as on $`\mathrm{\Lambda }_{\overline{MS}}r_0`$ in Tabs. 2629. $`\delta \mathrm{\Lambda }_m^{\mathrm{}}`$ refers to the contribution to the $`\mathrm{\Lambda }_{\overline{MS}}`$-parameters due to the mass dependence of the coefficients $`c_2`$, $`b_1`$ (through $`\mathrm{\Delta }K_1`$) and is included into the $`\mathrm{\Lambda }_{\overline{MS}}`$ estimates. This mass dependence has not been considered prior to the present study. For SW and KS fermions the effect turns out to be of the order of the statistical uncertainty as long as $`ma<0.1`$. However, for Wilson fermions this is a substantial effect, ranging from 10 % at $`ma0.06`$ to 2 % at $`ma0.01`$. In this case the main effect is due to the term $`\mathrm{\Delta }K_1(ma)`$ in the matching between lattice and $`\overline{MS}`$ scheme and we hence conclude that un-improved Wilson fermions are not particularly well suited for unambiguous determinations of the running coupling, unless data at rather small quark masses are available. We have neglected the mass dependence in our analysis of the SW-I data (Tab. 29), which is equally small as in the SW-Wilson case.

We use our perturbative result on the potential to estimate the ratio $`\mathrm{\Lambda }_{\overline{MS}}r_0`$ from the potential $`aV(1)`$, $`aV(2)`$ and force $`F_{12}`$ for the three $`n_f=2`$ data sets with Wilson glue. The resulting numbers are compared to the results from the plaquette in Tabs. 3032. Like in the quenched case the estimates from $`V(2)`$ consistently turn out to be bigger by about 15 %, in comparison to those obtained from $`V(1)`$ or the plaquette, however, in the case of KS quarks the estimates from $`F_{12}`$ are somewhat more in line with those from the plaquette. At present we do not have data at sufficiently many different lattice spacings for a detailed comparison like the one presented in Fig. 26 for the $`n_f=0`$ case. Hence, we will only use the most ultra violet quantity, the average plaquette, to extract the $`n_f=2`$ $`\mathrm{\Lambda }`$-parameter and assume a systematic uncertainty of 15 % due to higher order perturbative corrections, based on our study of the pure gauge case as well as on the numbers in Tabs. 3032.

In Fig. 27 we compare our $`n_f=2`$ estimates on $`\alpha _{\overline{MS}}(\mu )`$ from the plaquette to the $`n_f=0`$ results. Due to the variation in the sea quark masses, action and lattice spacings used, the $`n_f=2`$ points scatter around quite a bit but yield consistently larger values of $`\alpha `$, compared to the quenched case as already observed in lattice simulations of the static potential, for instance in Ref. Bali:2000vr  as well as in Sec. IV.3 below. It is also clear from the figure that at present the $`\mu `$-window covered by $`n_f=2`$ results is much smaller than the one available in pure gauge studies.

Since chiral symmetry is explicitly broken for the fermionic actions used we should first extrapolate our $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}`$ estimates obtained at similar quark masses in physical units $`mr_0`$ to the continuum limit, before attempting a chiral extrapolation: $`mm_u,m_d0`$. Due to the small range of lattice spacings covered by the presently available data, however, we are forced to interchange the ordering of these two limits: after discarding lattices with $`2.5a>r_0`$, i.e. the KS data with $`\beta <5.415`$ or $`ma>0.05`$, we extrapolate in $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0`$ at fixed $`\beta `$-values to small quark masses and subsequently send the lattice spacing to zero.

The chiral extrapolations are illustrated in Fig. 28 for Wilson glue and in Fig. 29 for the SW-I action. We include the $`\beta =1.8`$ data points into the latter figure but we will discard these coarse lattice results from our continuum limit extrapolation. The last column of Tab. 30 reveals that the slopes of $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0`$ as a function of $`mr_0`$ would have been much steeper for Wilson fermions had we neglected the $`ma`$ dependence of the matching between lattice and $`\overline{MS}`$ schemes. In the course of each extrapolation the lattice spacing in physical units $`a/r_0`$ changes. However, $`\mathrm{\Lambda }`$ is no spectral quantity and there is no theoretical handle on the functional form of such an extrapolation anyway. We find linear fits,

$$\left[\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0\right](mr_0)=\left[\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0\right](0)+cmr_0,$$
(111)
with parameters $`[\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0](0)`$ and $`c`$, to be consistent with the $`mr_0<0.25`$ ($`mr_0<0.35`$ for the SW-I action) data. The resulting parameter values as well as $`a/r_0(m=0)`$ can be found in Tab. 33. For the SW data at $`\beta =5.20`$ no fit is possible since $`N_{DF}=0`$ and at $`\beta =5.29`$ we included a rather heavy mass, $`mr_00.45`$.

Finally we attempt a continuum limit extrapolation in Fig. 30, where we have generously estimated error bars for the SW quarks at $`\beta =5.20`$. The lines represent extrapolations, linear in $`a/r_0`$ for Wilson fermions and quadratic in $`a/r_0`$ for KS and SW quarks. The extrapolated values are $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0=0.638(3)`$ with $`\chi ^2/N_{DF}=0.75/1`$ for KS quarks, $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_00.764`$ and $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_00.575`$ for Wilson and SW fermions, respectively, in which cases $`N_{DF}=0`$. For the SW-I action we obtain $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0=0.747(2)`$ with $`\chi ^2/N_{DF}=0.05/1`$. As a final result we quote the (unweighted) average of the KS and SW-I data sets,

$$\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0=0.69\pm 0.15,$$
(112)
where the error estimate reflects the uncertainties due to higher order perturbative corrections as well as in the extrapolations, both chiral and to the continuum limit.

The above number should be related to the recent more optimistic estimate $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0=0.553(34)`$ by the QCDSF and UKQCD collaborations Booth:2001qp  which is based on the UKQCD SW data alone. Note that their result agrees with our above value, $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_00.575`$, obtained from the same SW data set by use of somewhat different methods in the conversion between schemes and extrapolations. For instance in Ref. Booth:2001qp  $`\alpha _{\overline{MS}}(\mu )`$ was extracted at NNLO, estimating the as yet unknown NNLO fermionic contribution to the plaquette, while we consistently worked at NLO.

Fig. 30 paints a superficially disappointing picture for universality of the continuum limit. Universality absolutely must hold for lattice field theory to be a worthwhile exercise, so one might be concerned about the large difference in particular between the two different results obtained with SW fermions. While we try to estimate systematics by incorporating this difference into our error, it is interesting to consider how one could expect the figure to mature as improved data from future simulations become available.

Firstly, we note that the KS and SW-I continuum limit extrapolations are much better controlled than those for Wilson or SW quarks with Wilson glue, due to the availability of an additional data point. There are sizable $`𝒪(a^2)`$ corrections to the current KS and SW data. The Wilson continuum limit result has been obtained by drawing a straight line through two data points but there is no reason to believe $`𝒪(a^2)`$ corrections to be any smaller (nor indeed any larger) in this case than in the SW and KS cases. Likewise, it is not clear whether part of the large slope of the SW-I data set (which is only perturbatively improved such that $`𝒪(\alpha _La)`$ corrections are still present) might be due to a residual linear contribution. One might hope that the scatter of the results will be reduced by future data at different lattice spacings that would allow us to discriminate between linear and quadratic terms.

Secondly, we note that perturbative coefficients for the effects of the clover term in the SW action are particularly large throughout the literature, and refer the reader to Eq. (A.2) for an example of an $`𝒪(\alpha ^2)`$ coefficient and Eqs. (A.2) and (A.2) for $`𝒪(\alpha ^3)`$. There is no automatism that guarantees $`𝒪(\alpha ^3)`$ corrections to the plaquette in determinations of $`\alpha _{\overline{MS}}`$ to be small, just because this happened to be the case in the quenched theory with Wilson glue. Given the experience of the large $`c_{SW}^\nu `$ corrections we find it at least plausible that perturbative corrections to the SW determinations are larger than for the two other formulations. A calculation of these might result in a better agreement between the different continuum extrapolations too.

We remain confident that a combination of more lattice data and higher order perturbative calculations will greatly reduce the systematic spread shown in Fig. 30. Of course eventually the ordering of chiral and continuum limit extrapolations should be interchanged, unless data with chirally symmetric lattice fermions become available.

We finally note that efforts are on the way of calculating the $`n_f=2`$ $`\mathrm{\Lambda }`$-parameter with the help of finite size scaling techniques in the Schrödinger functional scheme. At present the error from this approach DellaMorte:2002vm  is of $`𝒪(50\text{\%})`$ and the scale has still to be related to a physical quantity like $`r_0`$ or $`f_\pi `$, however, some improvement can be expected in the near future.

#### IV.2.4 Implications for the “real” world

Running the estimate $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0=0.69(15)`$ to the mass of the $`Z`$ vector boson $`m_Z91.2`$ GeV results in the value,

$$\alpha _{\overline{MS}}^{(2)}(m_Z)=0.0945(30)(7),$$
(113)
where the second error that is smaller than 1 % is due to the 5 % scale ambiguity Bali:1997am  in the phenomenological value $`r_00.5`$ fm. This subdominant uncertainty can be reduced by using a quantity directly accessible in experiment like $`f_\pi `$ or $`m_N`$ to set the scale. Nature of course does not contain two but five quark flavours lighter than the $`Z`$ boson and we have to address this problem if we intend to produce numbers that have phenomenological relevance.

The main difficulty of a QCD prediction is the effect of the strange quark. Given that we cannot resolve the shift between $`\mathrm{\Lambda }_{\overline{MS}}^{(0)}r_0=0.60(5)`$ and $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0=0.69(15)`$ within the errors, for the moment being we assume that the strange quark which is heavier than the up and the down quarks anyway will not have a large impact on this ratio either and guess $`\mathrm{\Lambda }_{\overline{MS}}^{(3)}r_0\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0`$. However, we try to incorporate the uncertainty of this assumption (as well as the scale error of $`r_0`$) by somewhat inflating the systematic error:

$$\mathrm{\Lambda }_{\overline{MS}}^{(3)}=270(70)\text{MeV}.$$
(114)
Our result of Sec. III.5 indicates that a perturbative matching of the $`n_f=3`$ and $`n_f=4`$ couplings around the charm threshold should be reliable. Further, the difference between a three and a four-loop $`\overline{MS}`$ running from a scale of 1.1 GeV upwards corresponds to a relative shift in $`\alpha _{\overline{MS}}`$ at the $`Z`$ resonance of only 0.26 %. Hence, we perturbatively match the coupling at the charm and bottom thresholds. We allow pole masses $`1.1`$ GeV$`<m_c<1.5`$ GeV and 4.4 GeV $`<m_b<5.0`$ GeV and arrive at the estimate,

$$\alpha _{\overline{MS}}^{(5)}(m_Z)=0.1133(59),$$
(115)
using the four-loop $`\beta `$-function. The effect of varying the respective quark masses $`m_c`$ and $`m_b`$ as well as differences between a three and four-loop running of the coupling have been added linearly in the above error.

The result is compatible with the average quoted in the “Review of Particle Physics” Hagiwara:pw : $`\alpha _{\overline{MS}}^{(5)}(m_Z)=0.1172(20)`$ but lower than results from D0/CDF or LEP experiments, a tendency that is consistent with previous lattice estimates El-Khadra:1992vn ; davi1 ; davi2 ; Spitz:1999tu ; Booth:2001qp . In Fig. 31 we relate various lattice studies in chronological order to ours. We find that, although all lattice results on this quantity El-Khadra:1992vn ; davi1 ; davi2 ; Spitz:1999tu ; Booth:2001qp ; davi3  (as well as the PDG world average Hagiwara:pw ) lie within the error bars that we obtain, many of the previous numbers turn out to be incompatible with each other, within their quoted uncertainties. For instance we obtain a $`\chi ^2/N_{DF}=17.2/5`$ when averaging the six lattice results, and this assuming that the quoted errors, that are dominated by systematics, only represent one standard deviation!

The least well controlled step in our prediction Eq. (115) is certainly the extrapolation $`n_f=0,2n_f=3`$. A first lattice study with “$`n_f=3`$” exists davi3 , however, it is too early to investigate quark mass effects and to attempt a continuum limit extrapolation in this case. Once lattice results with 2+1 quark flavours of a similar quality as those analysed here become available a determination of $`\alpha _{\overline{MS}}(m_Z)`$ with about 4 % accuracy using the methods explored above should become realistic, and by both, a full fledged $`𝒪(\alpha ^3)`$ calculation in lattice perturbation theory and more simulation points that would allow for a better controlled continuum limit extrapolation this error should be reducible to 1 – 2 % within the next few years.

We experienced that systematic errors are easily underestimated in this kind of study, in particular since neither the size of higher order perturbative corrections nor the functional form of the chiral and continuum limit extrapolations are fully controlled. However, using a set of different lattice observables and actions certainly helps in arriving at realistic error estimates. Perturbation theory is naturally most reliable when the scales involved are high and it is very expensive to increase the energy $`a^1`$ in lattice simulations, at least if one keeps the simulation volume fixed in physical units. This makes the average plaquette a particularly valuable quantity in studies like the present one. Future predictions might benefit if lattice data on the plaquette, not only in the fundamental but also in the adjoint and higher representations, were available as well as data on other short distance quantities like the “chair”, the “parallelogram” and the “rectangle”.

### IV.3 Comparison to the non-perturbative static potential

We shall compare our perturbation theory results with data from simulations of QCD without and with sea quarks. First we explain what we mean by boosted and tadpole improved boosted perturbation theory in the present context. Next we investigate different possibilities to parametrize the pure gauge static potential at short distances, based on our results of Sec. III. We conclude with an attempt to resolve the different running of the coupling when including sea quarks from non-perturbative data.

#### IV.3.1 Parametrizations of the potential

In what follows we will refer to the parametrization

$$aV(𝐑a)=g^2V_1(𝐑)$$
(116)
with $`V_1`$ defined in Eq. (30) \[and Eq. (34)\] as leading order (LO). Apart from a different value of $`g^2=6/\beta `$ this result does not depend on $`n_f`$. The NLO result can be found in Eqs. (31) – (36). The LO result above corresponds to

$$aV(Ra)aV_S(a)C_F\frac{g^2}{4\pi R}(R\mathrm{}),$$
(117)
where $`aV_S=C_Fg^2\times 0.252731\mathrm{}`$ while the parametrizations of the NLO and NNLO curves in the $`R1`$ limit with and without massive sea quarks can be found in Appendix A.4. In the pure gauge case we will also compare against NNLO expectations for which the small $`R`$ lattice effects have not yet been calculated. When including fermions this is not possible since even in the massless case the contribution to the self energy is only known to NLO.

In addition to expansions in terms of the bare coupling we shall also compare against “boosted” perturbation theory expectations. In this case we use the $`\mathrm{\Lambda }_{\overline{MS}}`$ values of Tabs. 24 and 26, calculated from the respective plaquette at NLO as input. Since we shall compute the coupling $`\alpha _V(q)`$ at NNLO this is somewhat inconsistent. However, we also aim at a comparison with the fermionic case where such a determination of the $`\mathrm{\Lambda }`$-parameter is at present only possible at NLO and prefer to compare like with like. As we have seen above, the numerical difference between NLO and NNLO turns out to be small anyway for this quantity. Using $`\mathrm{\Lambda }_{\overline{MS}}r_0`$ as input we numerically run $`\alpha _{\overline{MS}}(\mu )`$ to the required scales $`\mu =e^{5/6}q^{}`$, according to the four-loop $`\beta `$-function. Some of the respective $`aq^{}(𝐑)`$ values can be read off Tab. 11 and for large $`R`$ we employ the parametrization Eq. (60). Slightly deviating from Sec. IV.2, the resulting $`\alpha _{\overline{MS}}(\mu )`$ is then always converted into $`\alpha _V(q^{})`$ at the two loop level \[$`𝒪(g^6)`$\], at bLO as well as at bNLO and bNNLO, via Eqs. (153) – (A.3) and Eq. (167). This means that what we call bLO includes some NLO and NNLO contributions and bNLO includes a conversion from the $`\overline{MS}`$ to the $`V`$ scheme at NNLO. However, the conversion from the momentum space continuum to the position space lattice potential is done at the order indicated.

Finally, we introduce “tadpole improved” boosted perturbation theory. We start from the observation that each link that appears within any lattice observable contains classes of tadpole diagrams that are lattice specific. We try to cancel these ultra violet contributions in part by a non-perturbative prescription Parisi:1980pe ; lepage : the second of the two Wilson loops that appear within the difference $`aV(a𝐑)=lim_T\mathrm{}\{\mathrm{ln}[W(𝐑,T)]\mathrm{ln}[W(𝐑,T+1)]\}`$ contains two more links than the first one. We hence prefer to perturbatively expand the combination $`aV(a𝐑)+\mathrm{ln}\mathrm{}/2`$ instead, within which the net number of links is zero. We have already made use of this trick in Eq. (74). Finally, we reexpress the result in terms of $`\alpha _V(q_N^{})`$, boosting our perturbative series, where the new $`q_N^{}(𝐑)`$ values can easily be calculated from those of the potential $`q^{}(𝐑)`$ and the plaquette $`q_{\mathrm{}}^{}`$:

$$\mathrm{ln}(aq_N^{}(𝐑))=\frac{\mathrm{ln}[aq^{}(𝐑)]V_1(𝐑)\mathrm{ln}(aq_{\mathrm{}}^{})c_1/2}{V_1(𝐑)c_1/2}.$$
(118)
The $`\alpha _V(q_N^{})`$ are then obtained as described above. We shall refer to this third method as tbLO, tbNLO and tbNNLO, depending on the order of the conversion from $`\alpha _V[q_N^{}(𝐑)]`$ to $`aV(𝐑)`$.

#### IV.3.2 Pure gauge theory

In Fig. 32 we compare LO, NLO, NNLO and bLO, bNLO and bNNLO expectations against the static potential at $`\beta =6.0`$. At this lattice spacing we use $`\mathrm{\Lambda }_{\overline{MS}}r_0=0.562(3)`$ and there is no free parameter. The un-boosted results come closer to the non-perturbative data as the order of the calculation increases but even at NNLO only about half of the potential can be explained, even at distances as short as $`R=1`$. The boosted results are more in line with the data but still only on a very qualitative level at best. It is also obvious that if one tuned the $`\mathrm{\Lambda }`$-parameter such that $`V(1)`$ is reproduced, $`V(2)`$ would be undershot, in agreement with our experience from Sec. IV.2 above. We can define an effective Coulomb coupling,

$$e(R)=R[aV(R)aV_S],$$
(119)
and find $`e0.106`$ at LO while $`e(3)0.19`$ and $`e(3)0.26`$ at NLO and NNLO, respectively. Finally, bNLO results in $`e(3)0.49`$.

Next in Fig. 33 we try out the tbLO, tbNLO and tbNNLO parametrizations and indeed tadpole improvement vastly reduces the difference with the NP potential. To exclude that this is just accidental we also display a comparison with the potential at $`\beta =6.4`$ in Fig. 34 where the plaquette yields $`\mathrm{\Lambda }_{\overline{MS}}r_0=0.603(10)`$, and with the smaller gauge coupling the situation becomes even more convincing. It turns out to be absolutely essential to compare like with like. Contrary to previous claims Bali:1999ai ; Necco:2001gh  not only the force but also the short range potential can be understood in terms of perturbation theory, provided that the self energy contribution is dealt with in a consistent way: when comparing with NP lattice results the self energy should not be subtracted from the perturbative expansion while when matching with the $`\overline{MS}`$ scheme in which $`V_S^{\overline{MS}}=0`$ by definition it has to be subtracted at the same scales $`\mu (r)`$ at which the interaction energy $`V_{\text{int}}(r)=V(\mu ;r)V_S(\mu )`$ is calculated, to avoid renormalon ambiguities. A similar observation has been made in a recent comparison with continuum perturbation theory in a renormalon based approach Pineda:2002se . In Ref. Bali:1999ai  a comparison was made between the NP lattice potential at short distances and the perturbative $`\overline{MS}`$ scheme prediction. In this study one and the same perturbative self energy had been used at all distances, which was obtained from a fit to the data, and consequently, no agreement was found.

In Fig. 35 we display the difference between non-perturbative and perturbative potentials, $`\mathrm{\Delta }V(r)=V(r)V^{\text{tbNLO}}(r)`$ for the example of $`\beta =6.4`$ as a function of $`R=r/a`$ in physical units $`r_0/a=9.89(16)`$. Within the statistical errors of the non-perturbative potential all points lie on a smooth curve indicating that the lattice effects are accounted for by the perturbative result. However, there are deviations between prediction and lattice data: at short distances we overestimate the NP potential by 30 MeV and around 0.3 fm we underestimate the result by 110 MeV. We fit the difference for $`r<0.61r_00.3`$ fm by the phenomenological ansatz,

$$\mathrm{\Delta }V(r)=b+cr^2$$
(120)
with $`b=0.068(3)r_0^1`$, $`c=0.99(6)r_0^3`$ and $`\chi ^2/N_{DF}=18.2/17`$. Allowing for an additional linear term does not improve the quality of the fit and the term turns out to be in agreement with zero: unlike in Ref. Bali:1999ai  no linear term is required at short distances to understand the data and this without any free parameter in the computation of $`\mathrm{\Delta }V(r)`$! We notice that the coefficient $`c(1.6\mathrm{\Lambda }_{\overline{MS}})^3`$ is of order one in units of the QCD $`\mathrm{\Lambda }`$-parameter.

In view of our plan to resolve the differences in the running of the coupling between quenched and un-quenched simulations, in Fig. 36 we display the dimensionless combination $`e(r/a)`$ of Eq. (119) for $`\beta =6.6`$ where $`5.1a0.4r_01\text{GeV}^1`$. At this $`\beta `$ value the average plaquette yields $`\mathrm{\Lambda }_{\overline{MS}}r_0=0.601(7)`$. We have also included bare LO, NLO and NNLO results for comparison. The self energies have been subtracted from all perturbative potentials, LO, NLO, NNLO as well as tbNLO at the respective orders in perturbation theory. As what we are doing is a linear transformation we do not wish to subtract different $`V_S`$ values from the NP and tbNLO potentials which would spoil the qualitative agreement between the results found above. Ideally one would subtract a non-perturbatively determined $`V_S`$ which then guarantees scaling between different NP data sets. The difference between NP and tbNLO however will then result in the tbNLO coupling to diverge linearly at large $`R`$ where perturbation theory will become unreliable anyway. In the absence of a non-perturbative determination of $`V_S`$ we will take the $`\beta =6.6`$ tbNLO value $`aV_S^{\text{tbNLO}}0.650`$ as an estimate which is indeed close at least to the tbNNLO result $`aV_S^{\text{tbNNLO}}0.670`$. Doing this means that the tbNLO Coulomb coupling will increase logarithmically with $`R`$ and no linear term is present.

In the figure we see a gap opening up between NP and tbNLO potentials for $`Ra>4a0.15`$ fm. The NP data starts diverging towards negative values, which is expected as $`e(r)\sigma r^2`$ as $`r\mathrm{}`$ where $`\sigma `$ denotes the string tension. In addition, our estimate of an NP self energy might be wrong by a few per cent which will result in an unwanted linear contribution to $`e`$. We find relatively large values of the Coulomb coupling, $`e0.4`$, that naturally disagree with $`e0.3`$ as obtained from phenomenological fits to a Cornell type potential Bali:2000vr ; Bali:2000gf , $`aV(R)=V_0e/R+a^2\sigma R`$, in which the linear confining term contributes at all distances: the parameter $`e`$ within this parametrization can only be interpreted as a Coulomb coupling at distances $`r\sqrt{\sigma }`$. We observe that the running of the coupling is quite steep at short distances which provides us with hope that differences between the quenched and un-quenched $`\beta `$-functions can be resolved from simulation data.

In Fig. 37 we display the same comparison for a similar $`R`$ range in physical units at $`\beta =6.0`$ where $`2.1a0.4r_01\text{GeV}^1`$. In this case the difference between tbNLO and tbNNLO self energies is 15 % such that we cannot regard these values as reliable anymore. Instead we compute $`a^{6.0}V_S^{6.0}=a^{6.0}V^{6.0}(r_0)(R_0^{6.6}/R_0^{6.0})a^{6.6}[V^{6.6}(r_0)V_S^{6.6}]`$, where the superscripts refer to the respective $`\beta `$ values, and add this non-perturbative difference to the $`\beta =6.6`$ $`aV_S`$ estimate above. In doing this we arrive at $`aV_S^{\text{NP}}0.83`$ which we subtract from both NP and tbNLO data sets. Now at large $`R`$ the tbNLO coupling will diverge linearly towards negative values as our $`V_S^{\text{NP}}`$ estimate differs from $`V_S^{\text{tbNLO}}`$ by as much as $`0.35a^1`$. It is obvious when comparing for instance the dimensionless NP couplings at $`R3.5`$ with $`R8.5`$ that both NP data sets scale nicely. We see however that at short distances the agreement with the perturbative expectation at $`\beta =6.0`$ is somewhat worse than that at $`\beta =6.6`$ above. The gap between NP and tbNLO curves opens up in the same region in lattice units, around $`R4`$, which corresponds to a larger distance in physical units. Since the perturbative expectation decreases at large $`r`$ only because of the difference between $`aV_S^{\text{NP}}`$ and $`aV_S^{\text{tbNLO}}`$ we should not expect agreement at $`R>\sqrt{5}`$ anymore and interpret this behaviour as accidental. This view is confirmed by our $`\beta =6.2`$ and $`\beta =6.4`$ data sets: the improved agreement at small $`r`$ goes along with the gap between perturbative and NP potentials opening up earlier.

#### IV.3.3 Comparison with $`n_f=2`$

We compare the quenched potential with a potential obtained with sea quarks for the example of two flavours of Wilson fermions at $`\beta =5.6`$ and $`\kappa =0.1565`$. These parameter values have been chosen to correspond to a quenched lattice coupling of $`\beta =6.0`$: in the $`n_f=2`$ case one finds Bali:2000vr  $`r_0=5.28(5)a`$ while a quenched $`\beta =6.0`$ yields $`r_0=5.33(3)a`$. The effective string tensions, determined by three parameter Cornell fits agree within errors too: $`\sigma a^2=0.0466(14)`$ vs. $`\sigma a^2=0.0479(7)`$, such that the linear term vanishes when subtracting the two potentials from each other. The $`\kappa `$ value corresponds to a quark mass $`ma0.040`$ and we will use this value to obtain the perturbative expectations.

In Fig. 38 we compare the NP $`n_f=2`$ potential with NLO and tbNLO perturbation theory. In the latter case we used the value $`\mathrm{\Lambda }_{\overline{MS}}r_0=0.570(6)`$, as obtained from the plaquette $`\mathrm{}=0.57073`$. The situation is very similar to that of the quenched potential at $`\beta =6.0`$ depicted in Figs. 32 and 33 which is why we refrain from including LO, bLO, bNLO and tbLO curves. The respective NNLO results are not known at present, even at large $`R`$.

In a next step we calculate the difference between the $`n_f=2`$ and the $`n_f=0`$ potentials. This is depicted in Fig. 39, together with the respective differences calculated in tbNLO and NLO. The NP linear confining term cancels in this combination, such that there is some reason to believe that this difference should be more perturbative than the potentials themselves, at least as long as $`R<12`$, the distance at which the QCD string will eventually “break” in the $`n_f=2`$ case. In the NP as well as perturbative data sets lattice artefacts are clearly visible that can be ascribed to the effect of the sea quarks. In the case of NLO we have depicted the $`n_f=0`$ expectation using the slightly smaller one loop matched value $`\beta =5.969`$ rather than $`\beta =6.0`$ for internal consistency.

As $`R\mathrm{}`$ all three curves approach constant values as the massive fermions start to decouple from the running. Although the shape of all curves is similar, even the tbNLO prediction underestimates the NP data by almost 30 %. This is not surprising as the vertical scale has been vastly inflated after subtracting two quantities of similar size from each other. For instance we can almost completely eliminate the gap between the NP and tbNLO results by increasing the $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0`$ value (or decreasing the $`\mathrm{\Lambda }_{\overline{MS}}^{(0)}r_0`$ value) by only 3 %. We find it quite encouraging that at least on a qualitative level we are able to resolve sea quark effects that only appear from $`𝒪(\alpha ^2)`$ onwards in perturbation theory. Unsurprisingly, a quantitative understanding is not possible without knowledge of the difference to second non-trivial order, i.e. $`𝒪(\alpha ^3)`$ in this case.

We now turn to the difference in the running of the two respective Coulomb couplings, between $`n_f=2`$ and $`n_f=0`$ and depict the quantity $`\mathrm{\Delta }e(R)=Ra[V^{(2)}(Ra)V^{(0)}(Ra)\mathrm{\Delta }V_S]`$ in Fig. 40. This combination further inflates scale and errors, in particular at large $`R`$. In the tbNLO case we can now identify deviations between the curve that constitutes the large $`R`$ limit and the points calculated in lattice perturbation theory at short distances. This bias, which also appeared in Fig. 16 and which turns out to be an $`𝒪(ma)`$ lattice artefact, has been discussed in Sec. III.3 above. Since at this lattice spacing, and within the scale of the figure, perturbative predictions and NP data cannot be expected to agree with each other anymore (cf. Fig. 39) we subtract different $`\mathrm{\Delta }V_S`$ values from the three data sets. While we compute $`a\mathrm{\Delta }V_S^{\text{NLO}}`$ and $`a\mathrm{\Delta }V_S^{\text{tbNLO}}`$ in perturbation theory, we fit $`a\mathrm{\Delta }V_S^{\text{NP}}0.062`$ from the NP data depicted in Fig. 39.

The NP data indicates an enhancement of the Coulomb coupling in the $`n_f=2`$ case at the given mass by almost 0.03 over the quenched coupling at short distances. This enhancement then reduces to zero at large distances where the two data sets are matched. Contrary to this we find the bare NLO expectation to be rather flat as a function of $`R`$: the running of the coupling with $`R`$ is faster when fermions are included, however, this effect is compensated for by a larger value of $`g^2`$, in comparison with the quenched case. This larger bare coupling can be translated into a deflation of the distance scale which in turn slows down the running again. In the end one would expect an increase of $`e`$ by 0.012 – 0.013 to NLO. The tbNLO prediction includes some slope but falls short of the NP data by a factor of about 3. After multiplying the points by this factor (open circles) the agreement with the NP $`\mathrm{\Delta }e(R)`$ is at best qualitative. We should not forget however that in the figure we have vastly inflated tiny differences. Unfortunately, $`n_f=2`$ lattice simulations at finer lattice spacings are not yet available such that we cannot test to what extent the quality of the parametrization improves as the continuum limit is approached.

## V Summary

We have calculated Wilson loops, the static potential and the static source self energy to $`𝒪(\alpha ^2)`$ in lattice perturbation theory with the Wilson gluonic action as well as with massive Wilson-Sheikholoslami-Wohlert and Kogut-Susskind fermions for arbitrary representation of the external sources. Some numbers for “improved” gluonic actions are provided too. The results are useful for understanding violations of rotational symmetry and the effect of including sea quarks on the lattice. Lattice results can be related to other regularization schemes like the perturbatively defined $`\overline{MS}`$ scheme, for instance to allow for an extraction of the QCD coupling at high energies from low energy hadronic phenomenology.

We find that at presently employed quark masses lattice spacing effects on the short distance potential are not smaller in absolute terms for the two non-perturbatively $`𝒪(a)`$ improved fermionic actions, relative to the Wilson action. However, improvement results in a reduced mass dependence around $`ma=0`$ of quantities like the static source self energy, that is required for determinations of the $`b`$ quark mass in an effective field theory framework, small Wilson loops and in the matching between lattice and $`\overline{MS}`$ schemes at finite $`a`$.

We perturbatively determined $`\beta =2N/g^2`$ shifts between the pure gauge theory and QCD with massive sea quarks. The results appear to be reliable within 10 % for quark masses down to 20 MeV for Wilson type fermions (which translates into less than 1 % uncertainty in the value of the predicted matched $`\beta `$) and within 25 % for “$`n_f=2`$” KS fermions within a window of lattice spacings that correspond to quenched $`\beta `$ values $`5.8<\beta <6.15`$. This method might turn out to be useful in predicting the impact of a shift in the sea quark mass onto correlation lengths like $`r_0/a`$ in future simulations. In particular we confirm that for SW fermions the lines of constant $`r_0/a`$ and $`m_\pi a`$ are significantly tilted with respect to those of constant $`\beta `$ and $`ma`$, respectively, a fact that complicates continuum limit extrapolations.

The success of the perturbative matching of the running coupling at flavour thresholds via an intermediate mass-dependent scheme, in our case the potential scheme, in lattice simulations (which allow us to compare with the correct non-perturbative result) means that the standard procedures employed in perturbative calculations Bernreuther:1983zp ; Brodsky:1998mf  should be quite reliable even at the charm quark mass.

We investigated so-called boosted perturbation theory methods Parisi:1980pe ; lepage  and found them to work well in many cases. However, such methods have to be digested with some caution: as boosted $`𝒪(\alpha ^2)`$ perturbation theory is tuned to reduce higher order and in particular $`𝒪(\alpha ^3)`$ corrections, differences between results calculated at these two orders are not necessarily indicative of the systematic uncertainties involved. In the absence of even higher order calculations, the perturbative uncertainty in for instance predictions of the running coupling can be estimated by varying both, the quantity from which the coupling is determined and the lattice spacing $`a`$.

In analysing quenched lattice data on the static potential and plaquette we obtain the value $`\mathrm{\Lambda }_{\overline{MS}}^{(0)}r_0=0.61(9)`$ which is in agreement with $`\mathrm{\Lambda }_{\overline{MS}}^{(0)}r_0=0.60(5)`$ obtained by use of finite size techniques by the ALPHA Collaboration Bode:1999hd . The error which is dominated by the uncertainty of higher order perturbative corrections can be systematically reduced by an $`𝒪(\alpha ^3)`$ calculation of the static potential. Using Wilson, KS and SW $`n_f=2`$ lattice data, the latter with Wilson as well as with Iwasaki glue, we estimate $`\mathrm{\Lambda }_{\overline{MS}}^{(2)}r_0=0.69(15)`$ in the continuum and chiral limit, where the error is again dominated by systematics. We find this value to depend strongly on both, the quark mass and the lattice spacing. Chiral and in particular continuum limit extrapolations are not fully under control but the use of different input quantities and lattice actions allows us to estimate the systematic uncertainties.

Based on the $`n_f=0`$ and $`n_f=2`$ estimates we guess $`\mathrm{\Lambda }_{\overline{MS}}^{(3)}=270(70)`$ MeV. This value corresponds to $`\alpha _{\overline{MS}}^{(5)}(m_Z)=0.1133(59)`$. In spite of the fact that more lattice data entered our analysis than has been used in any of the previous attempts to determine the $`n_f=2`$ running coupling, our error estimate is not smaller but larger than those that have been quoted by various lattice groups in the past El-Khadra:1992vn ; davi1 ; davi2 ; Spitz:1999tu ; Marcantonio:2001fc ; Booth:2001qp ; davi3 . The main reason is that the use of different quark actions and observables for the first time enabled us to realistically assess the uncertainties connected to higher order perturbative corrections and chiral and continuum extrapolations.

We estimate that the error on the QCD coupling at the $`Z`$ boson mass can be reduced to about 4 %, once a set of $`n_f=2+1`$ lattice simulations, similar to those that are already available for $`n_f=2`$, becomes available. By improving the perturbative calculation to NNLO as well as by incorporating a wider range of lattice spacings into the numerical simulations, errors of 1 – 2 % appear to be realistic within the next few years. It might be possible to further improve on this quality by a dedicated lattice simulation of running quark masses and coupling DellaMorte:2002vm .

We found that the lattice potential at distances smaller than about 1 GeV-1 is quite well described by NLO boosted tadpole improved perturbation theory without any free parameter, as long as no attempt is made to subtract the power divergent static self energy. This observation is somewhat in contrast to previous beliefs Bali:1999ai ; Necco:2001gh  and supports the view that perturbation theory for this quantity (and not only for the force) might be quite reliable in the continuum too. This is so as long as renormalons (in the $`\overline{MS}`$ scheme) and power divergences (in the lattice scheme) are properly dealt with when rearranging the perturbative series Pineda:2002se ; sumino . At distances larger than 1 GeV-1 the linear confining term sets in rather rapidly and it is doubtful that this behaviour can be emulated completely by higher order perturbative corrections. At distances shorter than 1 GeV-1 the difference between non-perturbative data and perturbative expectation is dominated by a quadratic term with a coefficient of $`𝒪(1)`$ in units of $`\mathrm{\Lambda }_{\overline{MS}}`$. The dependence of this coefficient on the way in which the perturbative series is organized, on the $`\mathrm{\Lambda }_{\overline{MS}}r_0`$ value that is used and on the lattice spacing deserves further study.

It turned out to be possible to resolve the difference in the logarithmic running of the coupling between two flavour QCD and the quenched approximation from non-perturbative data on the static potential. This difference is in rough qualitative agreement with the NLO expectation, however, about three times larger than expected. Unfortunately, at present we are limited to inverse lattice spacings of 2 GeV such that it is not possible to judge whether the discrepancy reduces as the scale is increased. As the magnitude of the effect should also increase as the quark mass is reduced, next generation lattice data at smaller quark masses should turn out helpful too.

###### Acknowledgements.

This work has been supported by the EU network HPRN-CT-2000-00145. During this work G.B. has benefited from a Heisenberg Fellowship (DFG grant Ba 1564/4-1) and a PPARC Advanced Fellowship (grant PPA/A/S/2000/00271). He has also received funding from PPARC grant PPA/G/O/1998/00559. P.B. acknowledges support from PPARC grant PPA/J/S/1998/00756. We thank Christine Davies, Antonio Pineda and Rainer Sommer for helpful comments and Sonali Tamhankar, Urs Heller and Alan Irving for making available to us MILC, HEMCGC and UKQCD results, respectively. We express our gratitude to the SESAM/T$`\chi `$L Collaboration for the permission to include some unpublished results on plaquette values and short distance potentials.

## Appendix A Perturbative relations between lattice and continuum schemes

We review the known perturbative relations between the lattice $`\beta `$-function and the modified minimal subtraction $`(\overline{MS})`$ scheme for various lattice actions in $`SU(N)`$ gauge theory with $`n_f`$ fermion flavours. We also display continuum perturbation theory results for the (mass dependent) static potential in position space, which constitute the limit against which our calculations will converge for $`ra`$.

### A.1 Conventions

We define the QCD $`\beta `$-function as,

$$\beta (\alpha )=\frac{d\alpha }{d\mathrm{ln}\mu ^2}=\beta _0\alpha ^2\beta _1\alpha ^3\beta _2\alpha ^4\mathrm{}.$$
(121)
While in mass-independent schemes like lattice regularization121212Strictly speaking, in lattice regularization, $`\beta _i`$ only become scale independent in the continuum limit $`a0`$, i.e. when scales $`\mu a^1`$ and quark masses $`m_ia^1`$ are considered. In general $`\beta _i`$ will obtain additional non-universal contributions that are functions of $`am_i`$ and $`a\mu `$. While contributions of the latter type are genuinely non-perturbative in character, quark mass dependent contributions arise in perturbation theory too. or minimal subtraction the first two coefficients,

$`\beta _0`$ $`=`$ $`\left({\displaystyle \frac{11}{3}}N{\displaystyle \frac{2}{3}}n_f\right){\displaystyle \frac{1}{4\pi }},`$ (122)
$`\beta _1`$ $`=`$ $`\left[{\displaystyle \frac{34}{3}}N^2\left({\displaystyle \frac{13}{3}}N{\displaystyle \frac{1}{N}}\right)n_f\right]{\displaystyle \frac{1}{(4\pi )^2}},`$ (123)
are universal, $`\beta _2`$ will in general depend on the renormalization scheme. In the $`\overline{MS}`$ scheme $`\beta _2`$ reads Tarasov:1980au ,

$`\beta _2^{\overline{MS}}`$ $`=`$ $`[{\displaystyle \frac{2857}{54}}N^3({\displaystyle \frac{1709}{54}}N^2{\displaystyle \frac{187}{36}}{\displaystyle \frac{1}{4N^2}})n_f`$ (124)
$`+({\displaystyle \frac{56}{27}}N{\displaystyle \frac{11}{18N}})n_f^2]{\displaystyle \frac{1}{(4\pi )^3}}.`$
The QCD $`\mathrm{\Lambda }`$-parameter in a given scheme is defined as,

$$\mathrm{\Lambda }=\underset{\mu \mathrm{}}{lim}\mu \mathrm{exp}\left(\frac{1}{2\beta _0\alpha (\mu )}\right)\left[\beta _0\alpha (\mu )\right]^{\frac{\beta _1}{2\beta _0^2}}.$$
(125)
We discuss the conversion between two different renormalization schemes: let,

$$\alpha ^{}(\mu )=\alpha (\mu )+c_1\alpha ^2(\mu )+c_2\alpha ^3(\mu )+\mathrm{}.$$
(126)
From the definition of the $`\beta `$-function, Eq. (121), one infers,

$`\beta _0^{}`$ $`=`$ $`\beta _0,\beta _1^{}=\beta _1,`$ (127)
$`\beta _2^{}`$ $`=`$ $`\beta _2c_1\beta _1+(c_2c_1^2)\beta _0,`$ (128)
$`\mathrm{\Lambda }^{}`$ $`=`$ $`\mathrm{\Lambda }e^{c_1/(2\beta _0)}.`$ (129)
If not only the scheme but also the matching scale $`\mu `$ is shifted, one easily obtains the relation,

$$\alpha ^{}(\mu ^{})=\alpha (\mu )+c_1^{}\alpha ^2(\mu )+c_2^{}\alpha ^3(\mu )+\mathrm{},$$
(130)
with,

$`c_1^{}`$ $`=`$ $`c_12\beta _0\mathrm{ln}\left({\displaystyle \frac{\mu ^{}}{\mu }}\right),`$ (131)
$`c_2^{}`$ $`=`$ $`c_2\left(2\beta _1+4c_1\beta _0\right)\mathrm{ln}\left({\displaystyle \frac{\mu ^{}}{\mu }}\right)`$
$`+4\beta _0^2\mathrm{ln}^2\left({\displaystyle \frac{\mu ^{}}{\mu }}\right)`$
$`=`$ $`c_22\beta _1\mathrm{ln}\left({\displaystyle \frac{\mu ^{}}{\mu }}\right)+c_{1}^{}{}_{}{}^{2}c_1^2.`$
### A.2 Conversion between the $`\overline{MS}`$ and lattice schemes

The $`\overline{MS}`$ coupling is related to the lattice coupling $`\alpha _L=g^2/(4\pi )`$ via,

$$\alpha _{\overline{MS}}(a^1)=\alpha _L+b_1\alpha _L^2+b_2\alpha ^3+\mathrm{},$$
(133)
where for the Wilson and Sheikholeslami-Wohlert (SW) quark actions one obtains Hasenfratz:1980kn ; Weisz:1981pu ; Luscher:1995np ; Alles:1998is ; Christou:1998ws ,

$$b_1=\frac{\pi }{2N}+k_1N+K_1(ma)n_f,$$
(134)
where $`K_1(ma)=K_1(0)+\mathrm{\Delta }K_1(ma)`$, with the numerical values,

$`k_1`$ $`=`$ $`2.135730074078457(2),`$ (135)
$`K_1(0)`$ $`=`$ $`0.08414443(6)`$
$`+0.0634188788(1)c_{SW}`$
$`0.3750240693(1)c_{SW}^2,`$
for the Wilson parameter $`r=1`$. The tree-level value of the SW coefficient, $`c_{SW}=1+3.3414(9)\alpha _L+\mathrm{}=c_{SW}^{(0)}+c_{SW}^{(1)}g^2+\mathrm{}`$, is $`c_{SW}^{(0)}=1`$. In the case of Kogut-Susskind (KS) fermions $`K_1(0)`$ reads Weisz:1981pu ; Luscher:1995zz ,

$$K_1^{KS}(0)=0.03298341916(1).$$
(137)
Note that in our notation $`n_f`$ KS quark flavours correspond to $`n_f`$ (not $`4n_f`$) flavours of Dirac quarks in the continuum limit. Therefore, $`n_f`$ has to be a multiple of four. The value of $`k_1`$ applies to the Wilson gluonic action only. However, the contribution $`K_1(ma)`$ (whose quark mass dependence has been calculated for the first time in this publication and is displayed in Tab. 10) remains the same for Symanzik type improvements of the Wilson gluonic action.

At vanishing quark mass or lattice spacing, $`ma^1`$ we obtain the massless limit, i.e. $`\mathrm{\Delta }K_1(0)=0`$. The mass dependent term results in a change of the QCD $`\beta `$-function at finite $`a`$:

$$\beta _0^L=\beta _0^{\overline{MS}}\underset{i}{}\frac{a}{2}\frac{d\mathrm{\Delta }K_1(m_ia)}{da}.$$
(138)
From this we can infer, $`lim_{ma0}d\mathrm{\Delta }K_1(ma)/d\mathrm{ln}a=lim_{ma0}d\mathrm{\Delta }K_1(ma)/d\mathrm{ln}(am)=0`$: the dependence is regular in $`ma`$ near the origin. Since the $`\beta `$-function is no on-shell quantity, we cannot exclude terms proportional to $`ma`$, even for improved quark actions. By demanding that in the limit $`ma\mathrm{}`$ the quenched $`\beta `$-function is re-obtained, we find,

$$\underset{x\mathrm{}}{lim}\frac{d\mathrm{\Delta }K_1(x)}{d\mathrm{ln}x}=\frac{1}{3\pi },$$
(139)
i.e. $`\mathrm{\Delta }K_1`$ will grow logarithmically for large $`ma`$. This is expected since $`ma\mathrm{}`$ means that the coupling on the left hand side of Eq. (133) runs with $`n_f`$ active flavours while the running on the lattice side of this equation corresponds to that of pure Yang-Mills theory. The difference has to be compensated for by $`\mathrm{\Delta }K_1(ma)`$. We will determine the integration constant in the region $`ma\mathrm{}`$ in Eq. (B.2) below from matching the un-quenched and quenched static potentials in the infra red.

The coefficient $`b_2`$ is at present only known for massless Wilson-SW quarks and Wilson gluonic action:

$$b_2=b_1^2+\frac{3\pi ^2}{8N^2}+k_2+k_3N^2+\left(\frac{K_2}{N}+K_3N\right)n_f$$
(140)
with Christou:1998ws ; Bode:2001uz

$`k_2`$ $`=`$ $`2.8626215972(6),`$ (141)
$`k_3`$ $`=`$ $`1.24911585(3),`$ (142)
$`K_2`$ $`=`$ $`0.1890(2)`$
$`0.02492(5)c_{SW}0.83585(3)c_{SW}^2`$
$`0.79942(5)c_{SW}^30.012947(2)c_{SW}^4,`$
$`K_3`$ $`=`$ $`0.1579(3)`$
$`0.00540(6)c_{SW}+0.7684(1)c_{SW}^2`$
$`+0.033842(5)c_{SW}^3+0.006920(2)c_{SW}^4.`$
From Eqs. (128) and (129) we infer,

$`\beta _2^L`$ $`=`$ $`\beta _2^{\overline{MS}}+b_1\beta _1(b_2b_1^2)\beta _0,`$ (145)
$`\mathrm{\Lambda }_{\overline{MS}}`$ $`=`$ $`\mathrm{\Lambda }_L\mathrm{exp}\left({\displaystyle \frac{b_1}{2\beta _0}}\right).`$ (146)
This results for instance in the numerical values,

$`\mathrm{\Lambda }_{\overline{MS}}^{(0)}`$ $``$ $`28.809\mathrm{\Lambda }_L^{(0)},`$ (147)
$`\mathrm{\Lambda }_{\overline{MS}}^{(2)}`$ $``$ $`41.053\mathrm{\Lambda }_L^{W,(2)}27.379\mathrm{\Lambda }_L^{SW,(2)},`$ (148)
for $`N=3`$ colours.

### A.3 Conversion between the $`\overline{MS}`$ scheme and schemes based on the static potential

We define the coupling $`\alpha _V`$ through the static QCD potential in momentum space,

$$\stackrel{~}{V}(q)=4\pi C_F\frac{\alpha _V(q)}{q^2},$$
(149)
where $`q=|𝐪|`$ and $`C_F=(N^21)/(2N)`$. The momentum space potential is related to that in position space via a Fourier transformation,

$$\stackrel{~}{V}(q)=\frac{d^3x}{(2\pi )^3}[V(r)V_S]e^{i\mathrm{𝐪𝐱}},$$
(150)
where $`r=|𝐱|`$. $`V_S`$ denotes the self energy term, that vanishes in dimensional regularization but will in general be present, and the position space potential,

$$V(r)=C_F\frac{\alpha _R(r^1)}{r}+V_S,$$
(151)
defines the coupling $`\alpha _R`$. By Fourier transformation one obtains the relation between the two schemes,

$$\alpha _R(\mu )=\alpha _V\left(e^\gamma \mu \right)\left(1+\frac{\pi ^2\beta _0^2}{3}\alpha ^2\right),$$
(152)
which holds for fermion masses $`m_i\mu `$. $`\gamma =0.57721566\mathrm{}`$ is the Euler constant.

We write,

$$\alpha _V(\mu )=\alpha _{\overline{MS}}(\mu )+a_1\alpha _{\overline{MS}}^2(\mu )+a_2\alpha ^3+\mathrm{}.$$
(153)
The coefficients $`a_1`$ and $`a_2`$ are known131313 In fact the leading log contribution to $`a_3`$ Brambilla:1999qa  and, more recently, the complete NNLO result Pineda:2000gz  have been calculated too.. For $`n_f`$ massless quark flavours they read Schroder:1999vy :

$`a_1`$ $`=`$ $`\left({\displaystyle \frac{31}{9}}N{\displaystyle \frac{10}{9}}n_f\right){\displaystyle \frac{1}{4\pi }},`$ (154)
$`a_2`$ $`=`$ $`\{({\displaystyle \frac{4343}{162}}+4\pi ^2{\displaystyle \frac{\pi ^4}{4}}+{\displaystyle \frac{22}{3}}\zeta _3)N^2`$ (155)
$``$ $`\left[\left({\displaystyle \frac{1798}{81}}+{\displaystyle \frac{56}{3}}\zeta _3\right){\displaystyle \frac{N}{2}}+\left({\displaystyle \frac{55}{3}}16\zeta _3\right){\displaystyle \frac{C_F}{2}}\right]n_f`$
$`+`$ $`{\displaystyle \frac{100}{81}}n_f^2\}{\displaystyle \frac{1}{16\pi ^2}}.`$
$`\zeta _3`$ denotes the value of the $`\zeta `$-function $`\zeta (3)=1.20205690\mathrm{}`$. The corresponding result for massive quark flavours with masses $`m_i`$ has been obtained in Ref. Melles:2000dq :

$`a_1(\{m_i\})`$ $`=`$ $`a_1+{\displaystyle \frac{1}{6\pi }}{\displaystyle \underset{i=1}{\overset{n_f}{}}}\mathrm{ln}\left(1+{\displaystyle \frac{C_0m_i^2}{\mu ^2}}\right),`$ (156)
$`a_2(\{m_i\})`$ $`=`$ $`a_2+{\displaystyle \frac{19}{48\pi ^2}}{\displaystyle \underset{i=1}{\overset{n_f}{}}}\{\mathrm{ln}(1+{\displaystyle \frac{C_3}{C_2}}{\displaystyle \frac{m_i^2}{\mu ^2}}+{\displaystyle \frac{1}{C_2}}{\displaystyle \frac{m_i^4}{\mu ^4}})`$
$`+{\displaystyle \frac{C_32C_1}{\sqrt{C_3^24C_2}}}\{2\mathrm{ln}(1+{\displaystyle \frac{C_3+\sqrt{C_3^24C_2}}{2C_2}}{\displaystyle \frac{m_i^2}{\mu ^2}})\mathrm{ln}[1+(2+{\displaystyle \frac{C_3}{C_2}}){\displaystyle \frac{m_i^2}{\mu ^2}}]\}\}`$
$`+{\displaystyle \frac{1}{54\pi ^2}}{\displaystyle \underset{i=1}{\overset{n_f}{}}}\left({\displaystyle \frac{31}{2}}N5n_f\right)\mathrm{ln}\left(1+{\displaystyle \frac{C_0m_i^2}{\mu ^2}}\right)+{\displaystyle \frac{1}{36\pi ^2}}{\displaystyle \underset{i,j=1}{\overset{n_f}{}}}\mathrm{ln}\left(1+{\displaystyle \frac{C_0m_i^2}{\mu ^2}}\right)\mathrm{ln}\left(1+{\displaystyle \frac{C_0m_j^2}{\mu ^2}}\right),`$
with the numerical values Brodsky:1999fr ,

$`C_0`$ $`=`$ $`5.19(3),`$ (158)
$`C_1`$ $`=`$ $`0.571(34),`$ (159)
$`C_2`$ $`=`$ $`0.221(15),`$ (160)
$`C_3`$ $`=`$ $`1.33(12).`$ (161)
A change in the argument on the left hand side of Eq. (153), $`\alpha _V(\mu )\alpha _V(\mu ^{})`$, results in the substitutions,

$`a_1(\mu ^{},\{m_i\})`$ $`=`$ $`a_1\left(\left\{m_i{\displaystyle \frac{\mu }{\mu ^{}}}\right\}\right)2\beta _0\mathrm{ln}\left({\displaystyle \frac{\mu ^{}}{\mu }}\right),`$ (162)
$`a_2(\mu ^{},\{m_i\})`$ $`=`$ $`a_2\left(\left\{m_i{\displaystyle \frac{\mu }{\mu ^{}}}\right\}\right)+4\beta _0^2\mathrm{ln}^2\left({\displaystyle \frac{\mu ^{}}{\mu }}\right)`$
$`\left[2\beta _1+4\beta _0a_1\left(\left\{m_i{\displaystyle \frac{\mu }{\mu ^{}}}\right\}\right)\right]\mathrm{ln}\left({\displaystyle \frac{\mu ^{}}{\mu }}\right),`$
in analogy to Eqs. (131) – (A.1).

The masses $`m_i`$ in Eqs. (156) – (A.3) denote pole masses,

$$m=\left[1+\frac{C_F}{4\pi }\left(43\mathrm{ln}\frac{m^2}{\mu ^2}\right)\alpha \right]m_{\overline{MS}}(\mu ).$$
(164)
Re-expressing Eq. (A.3) in terms of $`\overline{MS}`$ quark masses $`m_i`$ at a scale $`\mu `$, therefore, results in the replacement,

$$a_2^{\overline{MS}}(\{m_i\})=a_2(\{m_i\})+\frac{C_F}{12\pi ^2}\underset{i=1}{\overset{n_f}{}}\frac{43\mathrm{ln}\frac{m_i^2}{\mu ^2}}{1+\frac{\mu ^2}{C_0m_i^2}},$$
(165)
while Eq. (156) remains unaffected. For Wilson fermions one obtains Capitani:2001xi ,

$$m_{\overline{MS}}(a^1)=\left(1+16.95241\frac{C_F}{4\pi }\alpha \right)m_0.$$
(166)
Hence, in terms of lattice quark masses $`m_i(a)`$, $`a_2`$ reads:

$$a_2^L(\{m_i\})=a_2(\{m_i\})+\frac{C_F}{12\pi ^2}\underset{i=1}{\overset{n_f}{}}\frac{16.952413\mathrm{ln}m_i^2a^2}{1+\left(C_0m_i^2a^2\right)^1}.$$
(167)
We now go to position space and write,

$$\alpha _R(\mu )=\alpha _{\overline{MS}}(\mu )+a_1^R\alpha _{\overline{MS}}^2(\mu )+a_2^R\alpha ^3+\mathrm{}.$$
(168)
In following Eqs. (130) – (A.1) we obtain from Eqs. (152) and (153),

$`a_1^R`$ $`=`$ $`a_1+2\gamma \beta _0,`$ (169)
$`a_2^R`$ $`=`$ $`a_2+4\gamma \beta _0a_1+2\gamma \beta _1+\left(4\gamma ^2+{\displaystyle \frac{\pi ^2}{3}}\right)\beta _0^2.`$ (170)
The one loop result for massive quark flavours reads,

$$a_1^R(\{m_i\})=a_1^R+\frac{1}{3\pi }\underset{i=1}{\overset{n_f}{}}\text{Ein}\left(\frac{\sqrt{C_0}m_i}{\mu }\right)$$
(171)
with

$$\text{Ein}(x)=\gamma +\mathrm{ln}(x)+E_1(x)=\underset{\nu =1}{\overset{\mathrm{}}{}}\frac{(x)^\nu }{\nu \nu !}.$$
(172)
For the corresponding two-loop mass dependence of the position space potential we refer to Ref. Melles:2000dq . The exponential integral $`E_1(x)`$ is defined by,

$$E_1(x)=_x^{\mathrm{}}𝑑t\frac{e^t}{t},$$
(173)
i.e.

$$\text{Ein}(x)=_0^x𝑑t\frac{1e^t}{t}.$$
(174)
From Eqs. (126), (129), (153) and (154) one obtains the numerical values,

$$\mathrm{\Lambda }_V^{(0)}1.5995\mathrm{\Lambda }_{\overline{MS}}^{(0)},\mathrm{\Lambda }_V^{(2)}1.5213\mathrm{\Lambda }_{\overline{MS}}^{(2)},$$
(175)
while Eq. (152) results in,

$$\mathrm{\Lambda }_R^{(n_f)}=e^\gamma \mathrm{\Lambda }_V^{(n_f)}1.7811\mathrm{\Lambda }_V^{(n_f)}.$$
(176)
Note that for massive quarks the coefficients of the $`V`$\- and $`R`$-scheme $`\beta `$-functions become explicitly scale dependent. However, this does not affect the respective QCD $`\mathrm{\Lambda }`$-parameters, that are defined in the limit $`\mu \mathrm{}`$ \[Eq. (125)\].

We obtain:

$`{\displaystyle \frac{d\alpha _V}{d\mathrm{ln}\mu ^2}}=\left[\beta _0^{(n_f)}+{\displaystyle \frac{1}{6\pi }}{\displaystyle \underset{i=1}{\overset{n_f}{}}}\left(1+{\displaystyle \frac{\mu ^2}{C_0m_i^2}}\right)^1\right]\alpha _V^2(\mu )+\mathrm{}`$ (177)
$`{\displaystyle \frac{d\alpha _R}{d\mathrm{ln}\mu ^2}}=\left\{\beta _0^{(n_f)}+{\displaystyle \frac{1}{6\pi }}{\displaystyle \underset{i=1}{\overset{n_f}{}}}\left[1\mathrm{exp}\left({\displaystyle \frac{\sqrt{C_0}m_i}{\mu }}\right)\right]\right\}\alpha _R^2(\mu )+\mathrm{}.`$ (178)
In the situation, $`m_1=m_2=\mathrm{}=m_{n_f1}=0,m_{n_f}>0`$, the above coefficient functions smoothly interpolate between $`\beta _0^{(n_f1)}`$ (for $`\mu /m_{n_f}0`$) and $`\beta _0^{(n_f)}`$ (for $`\mu /m_{n_f}\mathrm{}`$) as one would expect.

### A.4 The perturbative potential at large $`R`$.

We denote a dimensionless distance measured in lattice units by $`R`$ while $`r=Ra`$ refers to the distance in physical units. At separations $`R1`$, i.e. $`ra`$, rotational invariance should be restored and the static potential obtained in lattice perturbation theory should coincide with the one computed above by use of dimensional regularization. In the large $`R`$ limit our results should therefore approach the continuum expressions given in Appendix A.3 above.

We define,

$`aV(𝐑a)`$ $`=`$ $`v_1(𝐑)\alpha _L+v_2(𝐑)\alpha _L^2+\mathrm{},`$ (179)
$`aV_S`$ $`=`$ $`v_1(\mathrm{})\alpha _L+v_2(\mathrm{})\alpha _L^2+\mathrm{},`$ (180)
where the self-energy $`V_S=2\delta m_{\text{stat}}`$ is twice the lattice pole mass of a static colour source. The coefficient functions $`v_i(𝐑)`$ have been calculated in the present article. At large separations we expect $`v_i(𝐑)`$ only to depend on the modulus $`R=|𝐑|`$ and can write,

$$v_{c,i}(R)=v_i(R)v_i(\mathrm{}).$$
(181)
Eq. (151) now implies,

$$a\left[V(𝐑a)V_S\right]C_F\frac{\alpha _R\left[(Ra)^1\right]}{R}(R\mathrm{}).$$
(182)
It follows,

$$v_1(𝐑)v_1(\mathrm{})v_{c,1}(R)=\frac{C_F}{R}(R\mathrm{}),$$
(183)
where

$$v_1(\mathrm{})=C_F\times 3.1759115\mathrm{},$$
(184)
for the Wilson gluonic action. $`v_2(\mathrm{})`$ is also known. Here we state the result for massless Wilson-SW quark flavours and $`N=3`$:

$`v_2(\mathrm{})`$ $`=`$ $`C_F\times \{16.714(4)+[0.42333(6)`$
$`+0.0516(2)c_{SW}0.5870(2)c_{SW}^2]n_f\}.`$
In the case of massless KS quarks the numerical value of the $`n_f`$ coefficient is $`C_F\times 0.36846(6)`$. Results on the mass dependence of the fermionic coefficients can be found in Tabs. 12 and 13. The $`n_f=0`$ result has first been derived in Refs. Weisz:1984bn ; Heller:1985hx . The precision has subsequently been increased by Martinelli and Sachrajda Martinelli:1999vt . We do, however, slightly disagree with this reference and obtain $`\delta m_{\text{stat},2}=C_Fv_2(\mathrm{})/2=11.143(3)`$, rather than 11.152 Martinelli:1999vt . The latter reference includes the Wilson and SW contributions for massless fermions as well, which we were able to reproduce with increased precision. The mass dependence of $`v_2(\mathrm{})`$ as well as the KS result are new. The $`n_f=0`$ value for $`v_3(\mathrm{})`$ has recently been obtained by means of stochastic perturbation theory DiRenzo:2001nd :

$$v_3(\mathrm{})=C_F\times 129(2).$$
(186)
This result agrees with the one of Ref. Trottier:2001vj  $`v_3(\mathrm{})=C_F\times 130(1)`$, obtained from a fit to high $`\beta `$ Monte Carlo data.

From Eqs. (179) – (183) we obtain,

$$\alpha _R\left[(Ra)^1\right]=\alpha _L+e_1(R)\alpha _L^2+e_2(R)\alpha _L^3+\mathrm{}$$
(187)
with

$$e_i(R)=\frac{v_{c,i+1}(R)}{v_{c,1}(R)}=\frac{R}{C_F}v_{c,i+1}(R).$$
(188)
In converting the scheme and scale by use of Eqs. (168), (133) and (130) – (A.1), we obtain the coefficients $`e_1(R)`$ and $`e_2(R)`$:

$`\alpha _R\left[(Ra)^1\right]`$ $`=`$ $`\alpha _{\overline{MS}}[(Ra)^1]+a_1^R(\{Rm_i\})\alpha _{\overline{MS}}^2[(Ra)^1]+a_2^R(\{Rm_i\})\alpha ^3`$ (189)
$`=`$ $`\alpha _L+[b_1+a_1^R(\{Rm_i\})+2\beta _0\mathrm{ln}R]\alpha _L^2+\{b_2+a_2^R(\{Rm_i\})+2b_1a_1^R(\{Rm_i\})`$
$`+`$ $`[2\beta _1+4(b_1+a_1^R(\{Rm_i\}))\beta _0]\mathrm{ln}R+4\beta _0^2\mathrm{ln}^2R\}\alpha ^3.`$
Therefore,

$`v_{c,2}(R)`$ $`=`$ $`{\displaystyle \frac{C_F}{R}}\left[b_1+a_1^R(\{Rm_i\})+2\beta _0\mathrm{ln}R\right],`$ (190)
$`v_{c,3}(R)`$ $`=`$ $`{\displaystyle \frac{C_F}{R}}\left\{b_2+a_2^R(\{Rm_i\})+2b_1a_1^R(\{Rm_i\})+[2\beta _1+4(b_1+a_1^R(\{Rm_i\}))\beta _0]\mathrm{ln}R+4\beta _0^2\mathrm{ln}^2R\right\}.`$ (191)
This amounts to

$$v_{c,2}(R)=\frac{C_F}{R}\left\{A_0+A_1\mathrm{ln}R+\underset{i=1}{\overset{n_f}{}}\left[\mathrm{\Delta }K_1(m_ia)+\frac{1}{3\pi }\text{Ein}(\sqrt{C_0}m_iaR)\right]\right\}$$
(192)
with141414We have factorized out the Casimir constant $`C_F`$ such that the above results also apply to the potential between charges in representation $`D`$, replacing $`C_F`$ by $`C_D`$. For $`N3`$ the constant $`C_0`$ has not yet been calculated.151515Rather than comparing the perturbative lattice potential at large distances with known results, one can also use Eq. (190) to determine the parameters $`b_i`$ from the calculated $`v_{i+1}(R)v_{i+1}(\mathrm{})`$ at large $`R`$. While in precision this procedure can in no way compete with the values displayed in Eqs. (A.2) or (137), calculating big Wilson loops and extracting the lattice potential by means of stochastic perturbation theory Parisi:1981ys ; DiRenzo:2001nd  might turn out to be a feasible alternative to diagrammatic techniques for more involved quark and gluonic actions.161616$`v_{c,i}`$ are subject to $`(am)^\nu `$ and $`(a/r)^\nu `$ lattice corrections, where $`\nu =1`$ in the case of Wilson fermions and $`\nu =2`$ for the KS and SW actions. While the $`r`$ dependent corrections vanish at large $`r`$ where we match to the continuum scheme, the $`am`$ corrections will in general be present in this limit. For $`0<ma<\mathrm{}`$ $`\mathrm{\Delta }K_1`$ can be defined by matching the lattice scheme to a mass-dependent scheme like the $`R`$-scheme. We do this by imposing that $`\mathrm{\Delta }K_1`$ contains all $`am`$ corrections to $`v_{c,2}`$. One can use different quantities to define $`\mathrm{\Delta }K_1`$, for instance the potential in momentum space. Since in general the $`(ma)^\nu `$ corrections for different observables will differ, $`\mathrm{\Delta }K_1`$ can only be fixed up to an $`(ma)^\nu `$ ambiguity. This is a reflection of the loss of universality of the $`\beta `$-function, due to the mass dependence of $`\beta _0`$ (and $`\beta _1`$). In the continuum limit, $`a0`$ universality will be restored. the numerical values for $`n_f`$ flavours of Wilson-SW fermions in $`SU(3)`$ QCD,

$`A_0`$ $``$ $`7.7164259753`$
$``$ $`[0.2338080.063419(1)c_{SW}+0.37524(1)c_{SW}^2]n_f,`$
$`A_1`$ $``$ $`1.750704374010.10610329540n_f.`$ (194)
For KS fermions the fermionic contribution to $`A_0`$ differs:

$$A_0^{KS}7.71642597530.1826473162n_f.$$
(195)
The function $`\text{Ein}(x)`$ is defined in Eq. (172) and can be approximated for small arguments by,

$$\frac{1}{R}\text{Ein}\left(\sqrt{C_0}m_iaR\right)=\sqrt{C_0}m_ia+𝒪\left[(m_iaR)^2\right]:$$
(196)
to leading order the coefficient $`v_{c,2}(R)`$ is parallel shifted in proportion to $`ma`$ at short distances.

At infra red scales, $`Ram_i^1`$, $`E_1(x)e^x/x`$, such that the massive flavours decouple from the logarithmic running of the potential. For $`n_f`$ massive degenerate flavours with masses $`m_i=m`$, we find:

$`v_{c,2}^{(n_f)}(R)`$ $`\stackrel{R\mathrm{}}{}v_{c,2}^{(0)}(R)`$
$`{\displaystyle \frac{n_fC_F}{R}}\left\{K_1(ma)+{\displaystyle \frac{1}{3\pi }}\left[\mathrm{ln}\left(\sqrt{C_0}ma\right){\displaystyle \frac{5}{6}}\right]\right\}.`$
For massless fermions we can parametrize the next order term in the expansion of the lattice potential $`v_{c,3}`$ by \[Eq. (191)\],

$$v_{c,3}=\frac{C_F}{R}\left(B_0+B_1\mathrm{ln}R+B_2\mathrm{ln}^2R\right).$$
(198)
The fermionic part of $`b_2`$, which is required to calculate the contribution to $`B_0`$ that is proportional to $`n_f`$, is only known for massless Wilson-SW fermions. Restricting ourselves to massless Wilson quarks and $`SU(3)`$ we obtain the numerical values,

$`B_0`$ $`73.8166085803`$
$`3.667(1)n_f+0.0530952n_f^2,`$
$`B_1`$ $`28.3102065046`$
$`2.295710n_f+0.0496157n_f^2,`$
$`B_2`$ $`3.06496580518`$
$`0.3715110067n_f+0.01125790929n_f^2.`$
$`B_2`$ is known exactly while the accuracy of the term within $`B_0`$ that is proportional to $`n_f`$ is limited by the precision of the constants $`K_2`$ and $`K_3`$ \[Eqs. (A.2) – (A.2)\].

## Appendix B The perturbative $`\beta `$-shift

### B.1 Matching quenched and un-quenched

In Eq. (A.4) we have seen that at distances $`rm^1`$ the running of the coupling is not affected by the presence of sea quarks anymore: at large distances, at least in perturbation theory, the effect of massive quarks can be integrated out into a shift of the coupling constant of the quenched theory. The matching can be done within an intermediate mass-dependent scheme, such as the $`V`$\- or $`R`$\- schemes. For simplicity we assume $`n_f`$ mass-degenerate flavours. The discussion below can easily be generalized to the non-degenerate case. For $`Ram^1`$ we write,

$`{\displaystyle \frac{a[V(Ra)V_S]}{v_{c,1}(R)}}`$ $`=`$ $`\alpha _L^{(n_f)}+e_1^{(n_f)}(R)\alpha _L^2+\mathrm{}`$ (202)
$`=`$ $`\alpha _L^{(0)}+e_1^{(0)}(R)\alpha _L^2+\mathrm{},`$
where $`e_i^{(n_f)}(R)`$ are the expansion coefficients of Eq. (187) for $`n_f`$ massive flavours with lattice quark mass $`m=(\kappa ^1\kappa _c^1)/(2a)`$. The matched coupling that results in the same physics at $`R(ma)^1`$ is,

$$\alpha _L^{(0)}=Z(m)\alpha _L^{(n_f)}$$
(203)
with

$$Z(m)=1+z_1(m)\alpha _L^{(n_f)}+z_2(m)\alpha _L^2+\mathrm{},$$
(204)
$`z_1(m)=`$ $`\underset{R\mathrm{}}{lim}\left[e_1^{(n_f)}(R)e_1^{(0)}(R)\right],`$ (205)
$`z_2(m)=`$ $`\underset{R\mathrm{}}{lim}\{e_2^{(n_f)}(R)e_2^{(0)}(R)`$ (206)
$`2e_1^{(0)}(R)[e_1^{(n_f)}(R)e_1^{(0)}(R)]\}.`$
In fact it turns out to be easier to determine $`Z(m)`$ in momentum space in the limit $`q0`$. Since one and the same physical potential is matched, momentum and position space results have to agree \[modulo the $`(am)^\nu `$ ambiguity within $`\mathrm{\Delta }K_1`$\] and we explicitly checked this in our calculation.

Using the $`\overline{MS}`$ scheme in momentum space we demand for $`q0`$:

$`\alpha _V^{(n_f)}(q)`$ $`=`$ $`\alpha _{\overline{MS}}^{(n_f)}\left(a^1\right)+a_1^{(n_f)}(q,\{m\})\alpha ^2+\mathrm{}`$
$`=`$ $`\alpha _{\overline{MS}}^{(0)}\left(a^1\right)+a_1^{(0)}(q)\alpha ^2+\mathrm{}`$
$`=`$ $`\alpha _V^{(0)}(q),`$
where the $`a_i(q,\{m\})`$ are defined in Eqs. (162) – (A.3). This amounts to,

$$\alpha _{\overline{MS}}^{(0)}\left(a^1\right)=W(m)\alpha _{\overline{MS}}^{(n_f)}\left(a^1\right)$$
(208)
with

$$W(m)=1+w_1(m)\alpha _{\overline{MS}}^{(n_f)}\left(a^1\right)+w_2(m)\alpha _{\overline{MS}}^2\left(a^1\right)+\mathrm{},$$
(209)
$`w_1(m)`$ $`=`$ $`\underset{q0}{lim}\left[a_1^{(n_f)}(q,\{m\})a_1^{(0)}(q)\right]`$
$`=`$ $`a_1^{(n_f)}a_1^{(0)}+{\displaystyle \frac{n_f}{3\pi }}\mathrm{ln}\left(\sqrt{C_0}ma\right),`$
$`w_2(m)`$ $`=`$ $`\underset{q0}{lim}[a_2^{(n_f)}(q,\{m\})a_2^{(0)}(q)`$
$`2a_1^{(0)}(q)w_1(m)].`$
By use of Eq. (133) one can then relate the above $`\overline{MS}`$ scheme coefficients to the lattice coefficients of Eq. (204):

$$z_1=w_1+b_1^{(n_f)}b_1^{(0)}=w_1+K_1n_f,$$
(212)
$`z_2`$ $`=`$ $`w_2+b_2^{(n_f)}b_2^{(0)}+2\left(b_1^{(n_f)}b_1^{(0)}\right)\left(w_1b_1^{(0)}\right)`$ (213)
$`=`$ $`w_2+\left({\displaystyle \frac{K_2}{N}}+K_3N+2K_1w_1+K_1^2n_f\right)n_f.`$
### B.2 Results

In following the above considerations we obtain the one loop result,

$$z_1(m)=\left[K_1(ma)\frac{5}{18\pi }+\frac{1}{3\pi }\mathrm{ln}\left(\sqrt{C_0}ma\right)\right]n_f:$$
(214)
the infra red behaviour of the theory with massive fermions is the same as that in the absence of massive flavours, at a gauge coupling that decreases with the quark mass:

$`\alpha _L^{(0)}`$ $`=`$ $`\alpha _L^{(n_f)}`$
$`+`$ $`{\displaystyle \frac{n_f}{3\pi }}\left[\mathrm{ln}(Dma)+3\pi \mathrm{\Delta }K_1(ma)\right]\alpha ^2,`$
$`D`$ $`=`$ $`\sqrt{C_0}\mathrm{exp}[3\pi K_1(0)5/6].`$ (216)
The constant $`D`$ has the numerical values,

$`D^W`$ $`=`$ $`0.448(2),`$ (217)
$`D^{SW}`$ $`=`$ $`0.0238(1),`$ (218)
$`D^{KS}`$ $`=`$ $`0.726(2),`$ (219)
for Wilson, SW and KS fermions, respectively. From the values of $`D`$ we see that the coupling shift is much more pronounced for SW fermions than it is for Wilson or KS flavours. From Eq. (139) we infer that for large $`ma`$ $`\mathrm{\Delta }K_1(ma)=\mathrm{ln}(D^{}ma)/(3\pi )`$. Note that Eq. (B.2) fixes the coefficient $`D^{}=D`$ and therefore the leading order behaviour of $`\mathrm{\Delta }K_1(ma)`$:

$`\mathrm{\Delta }K_1(ma)`$ $`\stackrel{ma\mathrm{}}{}`$ $`{\displaystyle \frac{1}{3\pi }}\mathrm{ln}(Dma)`$
$`=`$ $`0.00105(30)K_1(0){\displaystyle \frac{1}{3\pi }}\mathrm{ln}ma.`$
In the above limit as well as in the limit $`am=0`$ $`\mathrm{\Delta }K_1(ma)`$ is uniquely determined. This is different in the intermediate $`ma`$ range where $`\mathrm{\Delta }K_1(ma)`$ will depend on the quantity used in the matching and is therefore subject to $`𝒪[(am)^\nu ]`$ ambiguities. This in turn will also affect the predicted $`\beta `$-shift.

The coefficient of the $`\alpha ^2`$ term diverges for $`ma^1`$: for small values of $`m`$, the physics at distance scales $`rm^1`$, at which the matching to the quenched theory has to be performed, will be dominated by non-perturbative effects. Eventually, at very small quark masses the quenched and un-quenched theories cannot be matched to each other anymore at any distance scale with relevance to hadronic physics and the un-quenched theory becomes genuinely different from the quenched one.

Finally we state the result for $`z_2(m)`$:

$`z_2(m)`$ $`=`$ $`\{[({\displaystyle \frac{62}{9}}+{\displaystyle \frac{56}{3}}\zeta _3){\displaystyle \frac{N}{2}}({\displaystyle \frac{55}{3}}16\zeta _3){\displaystyle \frac{C_F}{2}}+{\displaystyle \frac{19}{3}}{\displaystyle \frac{C_32C_1}{\sqrt{C_3^24C_2}}}\mathrm{ln}\left({\displaystyle \frac{C_3+\sqrt{C_3^24C_2}}{C_3\sqrt{C_3^24C_2}}}\right)`$
$`+`$ $`{\displaystyle \frac{38}{3}}\mathrm{ln}\left({\displaystyle \frac{m^2a^2}{\sqrt{C_2}}}\right)]{\displaystyle \frac{1}{16\pi ^2}}+(16.952413\mathrm{ln}m^2a^2){\displaystyle \frac{C_F}{12\pi ^2}}+{\displaystyle \frac{K_2}{N}}+K_3N\}n_f`$
$`+`$ $`\left\{\left[{\displaystyle \frac{100}{81}}+{\displaystyle \frac{4}{9}}\mathrm{ln}^2\left(C_0m^2a^2\right)\right]{\displaystyle \frac{1}{16\pi ^2}}+\left[{\displaystyle \frac{20K_1}{9}}+\left({\displaystyle \frac{4K_1}{3}}{\displaystyle \frac{95}{81\pi }}\right)\mathrm{ln}\left(C_0m^2a^2\right)\right]{\displaystyle \frac{1}{4\pi }}+K_1^2\right\}n_f^2,`$
where $`m`$ denotes the lattice mass $`(\kappa ^1\kappa _c^1)/(2a)`$, rather than the pole mass (hence the correction term containing the numerical factor 16.95241). Unfortunately, we cannot yet apply this result to lattice data. For even in the case of Wilson-SW fermions where $`K_2(0)`$ and $`K_3(0)`$ are known, the quark mass dependence of these coefficients still has to be computed171717In the case of Wilson fermions one might be tempted to argue that $`𝒪(ma)`$ terms are lattice artefacts anyway. However, in the determination of the $`\beta `$-shift it is exactly this $`ma`$ dependence of $`r_0/a`$ that we are looking for and this indeed depends on the lattice action. Hence the variation of $`K_2`$ and $`K_3`$ with $`ma`$ cannot be neglected..