QCD phenomenology of static sources and gluonic excitations at short distances

New lattice data for the Pi_u and Sigma_u^- potentials at short distances are presented. We compare perturbation theory to the lower static hybrid potentials and find good agreement at short distances, once the renormalon ambiguities are accounted for. We use the non-perturbatively determined continuum-limit static hybrid and ground state potentials at short distances to determine the gluelump energies. The result is consistent with an estimate obtained from the gluelump data at finite lattice spacings. For the lightest gluelump, we obtain Lambda_B^{RS}(nu_f=2.5/r_0)=[2.25\pm 0.10(latt.)\pm 0.21 (th.)\pm 0.08 (Lambda_{MS})]/r_0 in the quenched approximation with r_0^{-1} approx 400 MeV. We show that, to quote sensible numbers for the absolute values of the gluelump energies, it is necessary to handle the singularities of the singlet and octet potentials in the Borel plane. We propose to subtract the renormalons of the short-distance matching coefficients, the potentials in this case. For the singlet potential the leading renormalon is already known and related to that of the pole mass, for the octet potential a new renormalon appears, which we approximately evaluate. We also apply our methods to heavy light mesons in the static limit and from the lattice simulations available in the literature we obtain the quenched result \bar{Lambda}^{RS}(nu_f=2.5/r_0)=[1.17\pm 0.08(latt.)\pm 0.13 (th.) \pm 0.09 (Lambda_{MS})]/r_0. We calculate m_{b,MS}(m_{b,MS}) and apply our methods to gluinonia whose dynamics are governed by the singlet potential between adjoint sources. We can exclude non-standard linear short-distance contributions to the static potentials, with good accuracy.


0
in the quenched approximation with r −1 0 ≈ 400 MeV. We show that, to quote sensible numbers for the absolute values of the gluelump energies, it is necessary to handle the singularities of the singlet and octet potentials in the Borel plane. We propose to subtract the renormalons of the short-distance matching coefficients, the potentials in this case. For the singlet potential the leading renormalon is already known and related to that of the pole mass, for the octet potential a new renormalon appears, which we approximately evaluate. We also apply our methods to heavy light mesons in the static limit and from the lattice simulations available in the literature we obtain the quenched result Λ RS (ν f = 2.5 r −1 0 ) = 1.17 ± 0.08(latt.) ± 0.13(th.) ± 0.09(Λ MS ) r −1 0 . We calculate m b,MS (m b,MS ) and apply our methods to gluinonia whose dynamics are governed by the singlet potential between adjoint sources. We can exclude non-standard linear short-distance contributions to the static potentials, with good accuracy.

Introduction
In recent years, we have witnessed growing interest in the physics of gluelumps and static hybrid potentials. In many cases this has been driven by increasingly reliable lattice simulations of their properties [1,2,3,4,5,6,7]. These results expose models of low energy QCD to stringent tests and therefore enhance our understanding of the underlying dynamics. The short distance physics of the static hybrid potentials is of particular importance. In this region, hybrids and gluelumps are intimately related and well suited to investigate the interplay between perturbative and non-perturbative physics. At short distances, r, one is faced with widely separated scales: 1/r ≫ Λ QCD . In such situations, effective field theories (EFTs) are particularly useful since they enable the physics associated with the different scales to be factorized in a very efficient and model independent way. One EFT designed to deal with the kinematical case of interest to us corresponds to potential non-relativistic QCD (pNRQCD) [8], in the static limit [9].
In Ref. [9] the gluelumps and the short distance regime of the static hybrids were studied within this EFT framework and general features identified. Some results known from the past [10,5,11,12,13] were recovered within a unified framework and in some cases extended.
One can go beyond this analysis and use lattice data plus the knowledge of the (perturbative) octet potential to obtain numerical values for gluelump masses in a particular scheme. However, analogously to the situation with the static singlet potential, the convergence of the perturbative series of the octet potential does not appear very promising. This is a general problem when different scales are factorized, and in particular perturbative from non-perturbative ones. The bad convergence is also related to the problem of factorizing non-perturbative quantities, without defining their perturbative counterparts [14], and is usually believed to be due to the existence of singularities in the Borel transform of the perturbative quantity. These singularities appear to be due to scales of order e −n ×-(the typical scale of the perturbative quantity) in an n-loop calculation. In Ref. [15] one of the present authors proposed that, since these singularities are related to energy scales much lower than the ones that are supposedly included in the perturbative object, they should be subtracted from it and introduced in the matrix elements of the effective theory. This programme has been worked out for the pole mass and the static singlet potential [15,16]. Here we apply the same approach to the static octet potential. This will allow us to determine absolute values for the gluelump masses from the spectrum of the static hybrids, as well as to study up to which scale one can use perturbation theory to describe hybrid potentials. This paper is organized as follows. In Sec. 2 we will work out the rôle of gluelumps in pNRQCD, and how gluelumps and hybrid potentials are interrelated. In Sec. 3 we will then sketch how our lattice data have been obtained, before discussing and classifying renormalons and power corrections in the continuum MS scheme as well as in a lattice scheme in Sec. 4. In the same section we will also generalise the renormalon subtracted (RS) scheme of Ref. [15] to the case of the octet potential and discuss the scale dependence. In Sec. 5 we will obtain the gluelump masses in the RS scheme and relate these results to the lattice scheme. We will compare to previous literature and predict the gluelump spectrum. In Sec. 6 we will determine the binding energy of static-light mesons as well as the bottom mass, before we discuss generalisations to and relations with adjoint potentials, gluinonium and other objects with relevance to short-distance QCD in Sec. 7.

Hybrid potentials and gluelumps
We discuss the relationship between hybrid potentials and gluelumps at short distances. First we consider the EFT picture, before we discuss the symmetries that are relevant in the non-perturbative case. Finally we compare these expectations to lattice data.

pNRQCD and gluelumps
The pNRQCD Lagrangian at leading order in 1/m and in the multipole expansion reads [8,9], All the gauge fields in Eq. (1) are evaluated in R and t, in particular F µν a ≡ F µν a (R, t) , O]. The singlet and octet potentials V i , i = s, o are to be regarded as matching coefficients, which depend on the scale ν us separating soft gluons from ultrasoft ones. In the static limit "soft" energies are of O(1/r) and "ultrasoft" energies are of O(α s /r). Notice that the hard scale, m, plays no rôle in this limit. The only assumption made so far concerns the size of r, i.e. 1/r ≫ Λ QCD , such that the potentials can be computed in perturbation theory. Also note that throughout this paper we will adopt a Minkowski space-time notation. The spectrum of the singlet state reads, where m OS denotes an on-shell (OS) mass. One would normally apply pNRQCD to quarkonia and in this case m OS represents the heavy quark pole mass. For the static hybrids, the spectrum reads where φ(T /2, −T /2) ≡ φ(T /2, R, −T /2, R) = P exp −ig and H represents some gluonic field, for examples see Table 4 in Sec. 5.3. Eq. (3) allows us to relate the energies of the static hybrids E H to the energies of the gluelumps, This equation encapsulates one of the central ideas of this paper. The combination E H − E s is renormalon-free in perturbation theory [up to possible O(r 2 ) effects], and can be calculated unambiguously non-perturbatively: the ultraviolet (UV) renormalons related to the infrared (IR) renormalons of twice the pole mass cancel each other. However, Λ H contains an UV renormalon that corresponds to the leading IR renormalon of V o . The shapes (of some) of the E H (r) have been computed on the lattice for instance in Refs. [1,2,3,4,5,6]. On the other hand the values of (some) Λ H have also been computed within a variety of models as well as in lattice simulations [12,17]. Consistency would require that the values of Λ H obtained from E H − E s and the values of Λ H directly obtained from gluelump computations should agree. This will be checked in Sec. 5.2 below.
Gluelump states are created by a static source in the octet (adjoint) representation attached to some gluonic content (H) such that the state becomes a singlet under gauge transformations. This is what would happen to heavy gluinos in the static approximation. Hence sometimes gluelumps are also referred to as gluinoballs or glueballinos in the literature. Without further information, their energy is only fixed up to a global constant. Only the energy splittings between different states have well defined continuum limits in lattice simulations. In lattice regularization at a lattice spacing a the normalization ambiguity is reflected in a linear divergence ∝ a −1 while in dimensional regularization one encounters an UV renormalon. In the HQET picture of a heavy-light meson one faces a similar problem. In this situation one also has a static source (in the fundamental representation in this case), which has to be attached to some light-quark (and gluonic) content to become a singlet under gauge transformations. The binding energy Λ is again only defined up to a global constant [18] and only its sum with the pole mass is unambiguous: We will investigate this situation in Sec. 6 below.

Symmetries of hybrid potentials and gluelumps
The spectrum of open QCD string states can be completely classified by the quantum numbers associated with the underlying symmetry group, up to radial excitations. In this case, these are the distance between the endpoints, the gauge group representation under which these endpoints transform (in what follows we consider the fundamental representation), and the symmetry group of cylindrical rotations with reflections D ∞h . The irreducible representations of the latter group are conventionally labelled by the spin along the axis, Λ, where Σ, Π, ∆ refer to Λ = 0, 1, 2, respectively, with a subscript η = g for gerade (even) P C = + or η = u for ungerade (odd) P C = − transformation properties. All Λ ≥ 1 representations are two-dimensional. The one-dimensional Σ representations have, in addition to the η quantum number, a σ v parity with respect to reflections on a plane that includes the two endpoints. This is reflected in an additional ± superscript. The state associated with the static singlet potential transforms according to the representation Σ + g while the two lowest lying hybrid potentials are within the Π u and Σ − u representations, respectively. In contrast, point-like QCD states are characterised by the J P C of the usual O(3) ⊗ C rotation group as well as by the gauge group representation of the source. In the pure gauge sector, gauge invariance requires this representation to have vanishing triality, such that the source can be screened to a singlet by the glue. States created by operators in the singlet representation are known as glueballs, octet states as gluelumps. In contrast to gluelump states, where the octet source propagates through the gluonic background, the normalization of glueball states with respect to the vacuum energy is unambiguous.
Since D ∞h ⊂ O(3) ⊗ C, in the limit r → 0 certain hybrid levels must become degenerate. For instance, in this limit, the Σ − u state corresponds to a J P C = 1 +− state with J z = 0 while the Π u doublet corresponds to its J z = ±1 partners. The gauge transformation property of the hybrid potential creation operator will also change in this limit, 3 ⊗ 3 * = 1 ⊕ 8, such that hybrids will either approach gluelumps [cf. Eq. (3)] or glueballs, in an appropriate normalization. In the case of glueballs the correct normalization can be obtained by considering the difference E H (r)−E s (r) from which the pole mass cancels. We will discuss the situation with respect to gluelumps in detail in Sec. 4 below.
In perturbation theory, the ground state potential corresponds to the singlet potential while hybrid potentials will have the perturbative expansion of the octet potential.
Recently, Philipsen [19] suggested to non-perturbatively generalise the octet potential, employing a definition that resembles the perturbative one, after gauge fixing to Laplacian Coulomb gauge. He proved that this construction is equivalent to a gauge invariant correlation function whose eigenvalues will resemble masses of physical states. In the limit r → 0 the suggested operator will be an adjoint temporal Schwinger line, dressed with a non-local but symmetric gluon cloud, with the J P C quantum numbers of the vacuum. A similar construction is mentioned in paragraph 2 of Sec. 6 below, as a possible non-perturbative normalization point for gluelump energies. The static "octet" potential suggested in Ref. [19] will have the Σ + g symmetry and, up to a different non-perturbative off-set, the same perturbative expansion and power term/renormalon structure as the hybrid potentials discussed below. Due to the nature of its creation operator which is non-local, even in the r = 0 limit, at present it is not obvious to us how this non-perturbative state can be interpreted in terms of the local states we are considering in this paper, certainly an open question that should be addressed in the future.

Hybrid and gluelump mass splittings
We would like to establish if lattice data on hybrid potentials reproduces the degeneracies expected from the above discussion in the short distance region. In the limit r → 0, any given Λ ≥ 1 hybrid potential can be subduced from any J P C state with J ≥ Λ and P C = + for η = g or P C = − for η = u representations. For instance the Π u is embedded in 1 +− , 1 −+ , 2 +− , 2 −+ , · · ·. The situation is somewhat different for Λ = 0 states, which have the additional σ v parity: the Σ + g representation can be obtained from 0 ++ , 1 −− , 2 ++ , · · ·, Σ − g from 0 −− , 1 ++ , · · ·, Σ + u from 0 +− , 1 −+ , · · · and Σ − u from 0 −+ , 1 +− , · · ·. We list all combinations of interest to us in Table 1. The ordering of low lying gluelumps has been established in Ref. [12] and reads with increasing mass: with a 3 −− state in the region of the 4 −− and 1 −+ . The 2 +− and 3 +− as well as the 4 −− and 1 −+ states are degenerate within present statistical uncertainties 1 . The continuum limit gluelump masses are displayed as circles at the left of Fig. 1, where we have added the (arbitrary) overall constant 2.26/r 0 to the gluelump splittings to match the hybrid potentials. The similarity of this value to our estimate of the gluelump energy in Sec. 5.1 below is purely accidental.
Juge, Kuti and Morningstar [1] have, for the first time, comprehensively determined the spectrum of hybrid potentials. We convert their data, computed at their smallest lattice spacing a σ ≈ 0.2 fm, into units of r 0 ≈ 0.5 fm [20]. Since the results have been obtained Table 1: Expected degeneracies of hybrid potentials at short distance, based on the level ordering of the gluelump spectrum. Note that if the 3 +− gluelump turned out to be lighter than the 2 +− then the with an improved action and on anisotropic lattices with a τ ≈ a σ /4, one might expect lattice artifacts to be small 2 , at least for the lower lying potentials. Hence we compare these data, normalized to E Σ + g (r 0 ), with the continuum expectations of the gluelumps [12]. The full lines are cubic splines to guide the eye while the dashed lines indicate the gluelumps towards which we would expect the respective potentials to converge. The first 7 hybrid potentials are compatible with the degeneracies suggested by Table 1. The next state is trickier since it is not clear whether 2 +− or 3 +− is lighter. In the figure we depict the case for a light 2 +− . This would mean that (Σ + u , Π ′ u , ∆ u ) approach the 2 +− while (Σ −′ u , Π ′′ u , ∆ ′ u , Φ u ) approach the 3 +− . Note that of the latter four potentials only data for Π ′′ u and Φ u are available. Also note that the continuum states Π ′ u , Π ′′ u and Φ u are all obtained from the same E u lattice representation. For the purpose of the figure we make an arbitrary choice to distribute the former three states among the E ′ u , E ′′ u and E ′′′ u lattice potentials. To firmly establish their ordering one would have to investigate radial excitations in additional lattice hybrid channels and/or clarify the gluelump spectrum in more detail. Should the 2 +− and 3 +− hybrid levels be inverted then (Σ −′ u , Π ′ u , ∆ u , Φ u ) will converge to the 3 +− while (Σ + u , Π ′′ u , ∆ ′ u ) will approach the 2 +− . We note that the ordering of the hybrid potentials, with a low Σ + u , makes the first interpretation more suggestive. Finally the Σ +′′ g potential seems to head towards the 0 ++ gluelump but suddenly turns downward, approaching the (lighter) sum of ground state potential and scalar glueball [21,22] instead. The latter type of decay will eventually happen for all lattice potentials but only 2 On the lattice the relevant symmetry group is D 4h rather than D ∞h (see e.g. Ref. [23]). In the continuum limit the A 1η potentials will correspond to Σ + η , the A 2η potentials to Σ − η and the E η potentials to Π η , where η = u, g. The radial excitations could in principal correspond to higher spin potentials and in fact one of the three observed excitations of E u will correspond to the Φ u ground state. In all other cases, associating the lowest possible continuum spin to a given lattice potential seems to agree with the ordering suggested by the gluelump spectrum (as well as in the large distance string limit [1]). B 1η and B 2η both correspond to ∆ η . In either case (as well as for ∆ ′ g ), at the short distances displayed in the figure, the two lattice representations agree with each other, supporting the view that violations of rotational symmetry are small. In this case we only display the lattice representations with better statistical accuracy, i.e. the B (′) 1η s.
-+ Figure 1: Different hybrid potentials [1] at a lattice spacing a σ ≈ 0.2 fm ≈ 0.4 r 0 , where r 0 ≈ 0.5 fm, in comparison with the gluelump spectrum, extrapolated to the continuum limit [12] (circles, left-most data points). The gluelump spectrum has been shifted by an arbitrary constant to adjust the 1 +− state with the Π u and Σ − u potentials at short distance. In addition, we include the sum of the ground state (Σ + g ) potential and the scalar glueball mass m 0 ++ [21,22]. The lines are drawn to guide the eye.
at extremely short distances. We also remark that all potentials will diverge as r → 0. This does not affect our comparison with the gluelump results, since we have normalized them to the Π u /Σ − u potentials at the shortest distance available. (The gluelump values are plotted at r = 0 to simplify the figure.) On a qualitative level the short-distance data are very consistent with the expected degeneracies. From the figure we see that at r ≈ 2 r 0 ≈ 1 fm the spectrum of hybrid potentials displays the equi-distant band structure one would qualitatively expect from a string picture. Clearly this region, as well as the cross-over region to the short-distance behaviour r 0 < r < 2 r 0 , cannot be expected to be within the perturbative domain: at best one can possibly imagine perturbation theory to be valid for the left-most 2 data points. With the exception of the Π u , Π ′ u and Φ u potentials there are also no clear signs for the onset of the short distance 1/r behaviour with a positive coefficient as expected from perturbation theory. Furthermore, most of the gaps within multiplets of hybrid potentials, that are to leading order indicative of the size of the non-perturbative r 2 term, are still quite significant, even at r = 0.4 r 0 ; for instance the difference between the Σ − u and Π u potentials at this smallest distance is about 0.28 r −1 0 ≈ 110 MeV.

2.4
The difference between the Π u and Σ − u hybrids From the above considerations it is clear that for a more quantitative study we need lattice data at shorter distances. In this paper we have obtained these for the lowest two gluonic excitations, Π u and Σ − u (see Sec. 3). We display their differences in the continuum limit in Fig. 2. We see how these approach zero at small r, as expected from the short distance expansion. pNRQCD predicts that the next effects should be of O(r 2 ) (and renormalon-free). In fact, we can fit the lattice data rather well with a ∆E Πu−Σ + g = A Πu−Σ − u r 2 ansatz for short distances, with slope (see Fig. 2), where the error is purely statistical (lattice). This fit has been done using points r < ∼ 0.5 r 0 . By increasing the fit range to r < ∼ 0.8 r 0 the following result is obtained, indicating stability of the result Eq. (8) above. In order to estimate systematic errors one can add a quartic term: b r 4 (only even powers of r appear in the multipole expansion of this quantity). If the result is stable, our determination of A Πu−Σ − u should not change much. Actually this is what happens. If we fit up to r < ∼ 0.5 r 0 , we obtain the central value A Πu−Σu r 3 0 = 0.93 with a very small quartic coefficient, b r 5 0 = −0.05. If we increase the range to r < ∼ 0.8 r 0 , we obtain the same central value, A Πu−Σ − u r 3 0 = 0.93, but with a slightly bigger quartic term, b r 5 0 = −0.18. Introducing the quartic term enhances the stability of A Πu−Σ − u under variations of the fit range. From this discussion we conclude that the systematic error is negligible, in comparison to the error displayed in our result Eq. (8).
We remark that within the framework of static pNRQCD and to second order in the multipole expansion, one can relate the slope A Πu−Σ − u to gluonic correlators of QCD.

Lattice determination of hybrid potentials
We extract the hybrid potentials in two sets of simulations, using the Wilson gauge action on an isotropic lattice with volume 24 3 × 48 at β = 6.2 (a ≈ 0.14 r 0 ) as well as on three anisotropic lattices with spatial lattice spacings a σ ≈ 0.33, 0.23, 0.16 r 0 , respectively, with anisotropy a σ ≈ 4 a τ . The former result has been obtained in the context of the study of Ref. [4] (and has been published in Ref. [3]) while the simulation parameters, statistics and smearing of the latter runs are identical to those of Ref. [24]: (β, ξ 0 ) = (5.8, 3.1), (6.0, 3.2), (6.2, 3.25). The isotropic data are used as a consistency check and in Sec. 4.4 below, while we extrapolate the data obtained on the anisotropic lattices to the continuum limit. Some time was spent on improving the shape of the hybrid creation operators to optimize the overlap with the ground state [4]. The Π u potential has been determined on-axis as well  as along a plane-diagonal, r/a σ ∝ (1, 1, 0), while the Σ − u potential has only been obtained on-axis. Typically we achieved ground state overlaps of around 65 % for both potentials at β = 5.8 and between 85 % and 90 % at the larger two β values. Typical fit ranges for one-exponential fits to correlation functions for the Π u (Σ − u ) potential were 8 ≤ t/a τ ≤ 18, (9 ≤ t/a τ ≤ 14) at β = 5.8, 9 ≤ t/a τ ≤ 24 (11 ≤ t/a τ ≤ 21) at β = 6.0 and 13 ≤ t/a τ ≤ 30 (15 ≤ t/a τ ≤ 25) at β = 6.2. For all further details of the analysis we refer to Ref. [24] where potentials between sources in non-fundamental representations of SU(3) were extracted using exactly the same methods.
Subsequently, the potentials as well as differences between potentials have been extrapolated to the continuum limit. As one such example we display the difference between the Π u and the singlet potential in Fig. 3 for distances r ≤ r 0 . In this extrapolation we somewhat deviate from Ref. [24]: we follow Ref. [25] in removing the lattice artifacts to leading order in α s , by plotting the data as a function of the inverse lattice Coulomb propagator,  rather than of r. The lattice Coulomb propagator for the Wilson gauge action is given by, and agrees with the continuum 1/R-function up to O(a 2 /r 2 ) lattice artifacts. R = r/a denotes an integer valued three-vector and the Q i = q i a σ are dimensionless. For the Π u potential this procedure removes violations of rotational symmetry within the statistical errors and brings the plane-diagonal points in-line with the on-axis data. Unfortunately, we cannot perform a similar internal test for the Σ − u potential which we only determined for on-axis separations.
The next step involved fitting differences between hybrid potentials and Σ + g , ∆E H = E H − E Σ + g , for r ≥ 2a to the phenomenological interpolation, with parameters c i . We then extrapolated these interpolating curves to the continuum limit, assuming the leading order a 2 σ dependence. This was done separately for different pairs of two lattice spacings. The central value of the extrapolation is given by the result obtained from the a σ ≈ 0.33 r 0 and a σ ≈ 0.16 r 0 data sets. The error is estimated by the squared sum of the statistical error of the fine lattice data set and the difference between the above extrapolation and an extrapolation obtained from the a σ ≈ 0.23 r 0 and a σ ≈ 0.16 r 0 data sets. With decreasing r the interpolating fits become less well constrained and hence the latter systematic uncertainty increases. The resulting error band is depicted in Fig. 3. Reassuringly, the a σ ≈ 0.16 r 0 data are already in agreement with the continuum limit and the a σ ≈ 0.23 r 0 data agree within errors: the fine lattice data set effectively already corresponds to the continuum limit. The more precise isotropic reference data (a ≈ 0.14 r 0 ) are also close to the continuum limit. We also notice that the first three data points of the coarse lattice data by Juge et al. [1] (a σ ≈ 0.39 r 0 ) are compatible with our extrapolation. The same observations hold true for the Σ − u potentials. Rather than representing the continuum limit extrapolated potentials by error bands, in the remaining parts of this paper we add the difference between (finite a) interpolation and (continuum limit) extrapolation to the fine lattice data points and increase their errors by the systematic uncertainty involved in the extrapolation.

Static octet potential
We will discuss the octet potential in the OS (= "pole mass") scheme, compute the normalization constant of the renormalon and generalise the RS renormalon subtracted scheme [15] to this case. We will also discuss the structure of power divergences on the lattice and the analogous lattice scheme. Finally we discuss the running of the gluelump mass from one scale to another.

OS scheme for the octet potential
The octet potential in the case 1/r ≫ Λ QCD can be computed order by order in perturbation theory. Nevertheless, it is not an IR safe object [26]. Its perturbative expansion reads, where we have made explicit its dependence on the IR cutoff ν us and α s = α s (ν), where we define In what follows we will always identify α s with α MS . The first two coefficients V o,0 , V o,1 are known, as well as the leading-log terms of V o,3 [26] (for the renormalization-group improved expression see Ref. [27]). Note, however, that these leading logs are not associated to the first IR renormalon. For V o,2 there exists a preliminary computation [28], which we will use in what follows. V s,2 has been computed in Ref. [29]. For V o,3 , we will use the renormalon-based estimate that we obtain in Sec. 4.2 below (Table 2). Studying the convergence of perturbation theory of the octet potential in the OS scheme, conclusions similar to those in Ref. [16] are obtained. The poor convergence is demonstrated in Fig. 4, where we try two choices of the scale ν. In part (a) we use ν = ν i where, corresponds to the shortest distance r ′ for which the continuum limit extrapolated lattice potentials are available. In part (b) we vary ν = 1/r. Obviously the curves depicted in the two parts of the figure agree with each other at r = r ′ ≈ 0.15 r 0 . Note the difference in the vertical scale.

Static octet potential normalization constant
We define the Borel transform of the octet potential as follows, The behaviour of the perturbative expansion Eq. (13) at large orders is dictated by the closest singularity to the origin of its Borel transform, which happens to be located at t = 2π/β 0 . This singularity has two sources: one is an UV renormalon which cancels with the renormalon of twice the pole mass, the other is an IR renormalon that cancels with the UV renormalon of the gluelump energy. This result follows from the structure of the effective theory and the consequent factorization of the different scales in Eq. (3). Being more precise, the behaviour of the Borel transform of the static octet potential near the closest singularity to the origin [u = 1/2 where we define u = β 0 t/(4π)] reads, where by analytic term, we mean a function that is analytic up to the next IR renormalon at u = 3/2. This dictates the behaviour of the perturbative expansion at large orders to be, The structure of the renormalon is equal to the singlet one. This is due to the fact that the number of octet fields is conserved at leading order in the multipole expansion and that the mass (potential) does not renormalize at this order. Therefore the values of the coefficients b, c 1 , c 2 , . . . above are the same as for the case of the static potential and the pole mass and can be found in Refs. [30,31,15]. We display them here for ease of reference: and The only difference with respect to the static singlet potential is the value of N Vo . The cancellation of the renormalon in Eq. (3) requires, where N Λ is the normalization constant of the renormalon of the gluelump mass (B[Λ] reads the same as Eq. (18), with the replacement N Vo → N Λ ). Therefore, unlike in the static singlet potential case, we cannot fix N Vo from the knowledge of N m alone. Yet we will (approximately) determine N Vo from low orders in perturbation theory of the octet potential. Note also that N Λ is independent of H, the specific gluonic content of the gluelump, since it only depends on the high energy behaviour, which is universal. To leading non-trivial order one obtains, In analogy to Refs. [32,15,16] we define the new function, and try to approximately determine N Vo by using the first three coefficients of this series. In analogy to Refs. [15,16], we fix ν = 1/r and obtain [up to O(u 3 )| u=1/2 ], The convergence is rather good and, moreover, we have a sign alternating series. In fact, the scale dependence is becoming milder when we go to higher orders (see Fig. 5). Note that if the two loop coefficient V o,2 had been equal to that of the singlet case [29] (with colour factor We can now compute estimates for V o,n by using Eq. (19). These, as well as estimates for V s,n , are displayed in Table 2 for n f = 0. We can see that the exact results are reasonably well reproduced. Hence we feel confident that we are near the asymptotic regime dominated by the first IR renormalon and that for higher n our predictions will accurately approximate the exact results. Table 2: Values of V o,n with ν = 1/r: exact result (where available) and the estimate using Eq. (19). We also display the estimates of V s,n with ν = 1/r (extracted from Ref. [15]).
In order to avoid large corrections from terms depending on ν us , the predictions should be understood with ν us = 1/r and later-on one can use the renormalization group equations for the static potential [27] to keep track of the ν us dependence.

RS scheme for the octet potential
In Sec. 4.1 we have demonstrated the poor convergence of the perturbative expansion of the octet potential in the OS scheme. This bad behaviour is usually believed to be due to the singularities in the Borel transform of the perturbative expansion. Nevertheless, these singularities are fake since they cancel with singularities in the matrix elements. On the other hand this lack of convergence of perturbation theory arises because at higher orders in perturbation theory smaller and smaller momenta contribute to the short-distance matching coefficients of the effective theory. This clashes with the logic of scale separation in the EFT formalism. The solution advocated in Ref. [15] was to subtract this behaviour from the matching coefficients. At the practical level this was implemented by subtracting the Borel plane singularities of the matching coefficients. In Refs. [15,16] this has been worked out for the pole mass and the static singlet potential and we refer to these references for the definitions and further details. In particular Eq. (2) reads, where and (in the above equation we have already used the fact that the renormalon of the singlet potential cancels with the one of minus twice the pole mass) 3 , For the static hybrids, the spectrum reads, Obviously, we have to define the octet potential and the gluelump mass above. In the RS scheme the octet potential reads, where This specifies the gluelump mass which reads, Note that the potentials and Λ RS H depend on ν f which, in the context of pNRQCD, can be thought of as a matching scale between ultrasoft and soft physics. In what follows, we will set ν f = 2.5 r −1 0 . Results for different values of ν f can be obtained using the running on ν f , which is renormalon independent.
Analogously to the discussion of Ref. [16], we can study the convergence of the perturbative expansion in the RS scheme. In Fig. 6 we can see that the stability is greatly improved, compared to the OS scheme discussed in the previous section. No matter whether we choose to work with α s (ν i ) or α s (1/r), the expansions converge towards the same curve. In Fig. 7 we can also see that they agree with the continuum limit lattice data (we have to subtract an unknown constant for this comparison). In this figure the errors of E Πu (r) − E Πu (r ′ ) for r > r ′ are purely statistical while the (strongly correlated) systematic error of the continuum limit extrapolation is only displayed for the first data point The price we pay to obtain convergent expansions in α s for the potentials is the introduction of power-like terms (proportional to ν f , with logarithmic corrections). This behaviour very much resembles that of lattice regularization with a hard cut-off which we discuss below.
, two loops (dotted lines) and three loops (estimate) plus the leading single ultrasoft log (solid lines) compared with the non-perturbative continuum-limit results for E Πu (r) − E Πu (r ′ ) (symbols with error bars). For the scale of α s (ν), we set ν = ν i = 1/r ′ (stable behaviour at large distances) or ν = 1/r (diverging at large distances). A (small) constant C is arbitrarily adjusted to show agreement with the lattice data.

Lattice scheme for the octet potential
It is conceptionally illuminating also to consider the situation in lattice regularization. In this case, the inverse lattice spacing a −1 results in a hard UV cut-off of the gluon momenta. Feynman diagrams are UV finite and EFT matrix elements are manifestly renormalon-free as long as they are obtained in non-perturbative numerical simulations. The price paid is the existence of power divergences ∝ a −1 , which cannot be eliminated in the continuum limit.
The analogy with the previous sections can be made quite evident. In particular, all the quantities that we have defined in the OS and RS schemes can also be defined in a lattice scheme. There are some differences however. The lattice gluelump Λ L H (a) has a power divergence to start with (which can be traded in for a renormalon ambiguity when subtracted in perturbation theory). In this sense it is similar to Λ RS H (ν f ). While formally many expressions resemble those of the RS case, a −1 plays a slightly different rôle than ν f that separates soft from ultrasoft scales since a ≪ r ≪ ν −1 f . Another difference is that at finite lattice spacings the potentials remain finite as r → 0. In particular, we will see that gluelumps are the r → 0 limits of hybrid potentials (at finite lattice spacing), in perturbation theory as well as non-perturbatively. This should not be surprising since the r → 0 limit at finite lattice spacing corresponds to the situation r ≪ a. This means that the ultraviolet cutoff ∼ a −1 is much smaller than r −1 and that the dynamical degrees of freedom are only the ultrasoft ones. Actually, in this situation, ν f and 1/a play an analogous rôle.
Let us illustrate the above by first considering perturbation theory, before discussing the scale separation and how the lattice scheme translates into other schemes, at finite lattice spacings as well as in the continuum limit. For simplicity we will consider the Wilson discretization of the continuum action. In this case the "lattice Coulomb term" [1/R] L takes the form Eq. (11). For instance one can calculate the finite value, [1/0] L = 3.17 . . .. Using this notation, one finds the lattice results, where the "self energy" is given by, Note that unlike in dimensional regularization, by using a hard cut-off, such power divergencies appear naturally as part of the perturbative expansion. Eqs. (35) and (36) [33]. In pure gauge theory with Wilson action, the coefficients of the expansion of δm L stat read [34,18,33,35], α L = 3/(2πβ) denotes the lattice coupling at a scale a −1 which can be translated into other schemes such as MS by means of a perturbative computation, with [36], Let us now consider the singlet case. We have, where Λ QCD represents a generic non-perturbative scale like r −1 0 . The last two terms account for possible non-perturbative lattice artifacts, which vanish as a → 0. From the quarkonium energy E s (r) at r ≫ a, we can non-perturbatively obtain the heavy quark mass in a lattice scheme, By redefining, we can then achieve formal correspondence to Eqs. (26) and (2), respectively: where the above two equations are correct up to O(Λ 2 QCD a 2 ) and O(a 2 /r 2 ) lattice corrections. We can relate the heavy quark mass in the lattice scheme to the OS scheme, δm L stat (a) contains the same renormalon as m OS , such that Eq. (49) has good convergence properties when expanded in terms of α s . m L (a) is proportional to a −1 , with logarithmic, as well as O(a 2 ) lattice corrections. One can convert m L (a) order by order in perturbation theory into say m MS (ν), without renormalon ambiguity.
In the lattice scheme we also have E L Σ + g (0; a) = V s,L (0; a) = 0: the sources are "smeared out" on a scale a since the gluon, due to the UV cut-off, cannot resolve structures smaller than the lattice spacing. Consequently, the Coulomb term does not diverge as r → 0 but approaches a finite value in units of a. In perturbation theory, in the limit r → 0, the lattice [1/R] L term exactly cancels with 2δm L stat : the perturbative expansion of V s,L (r; a), Eq. (35) above, does not contain the renormalon associated with the pole mass. Nonperturbatively, in the limit r → 0 the Wilson loop becomes a time independent constant, such that E L Σ + g (0; a) = 0 too. As r > 0 the perturbative V s,L acquires a power term.
Next we consider the hybrid case. We can calculate the gluelump mass in perturbation theory 4 , where v 1 and v 2 are the same as for the case of δm L stat and can be found in Eqs. (38) and (39) above. Note that the O(α 3 s ) term is expected to be different and is not known at present. However, Λ H is related to the difference between V o and V s , such that any difference with respect to the v 3 of Eq. (40) above will be suppressed by a colour factor 1/N 2 c . The tree level expression for V o,L is displayed in Eq. (36). While the perturbative expansion of V s,L was unaffected by the renormalon of the pole mass, the one of V o,L contains the same renormalon as the expansion of δΛ L . For r ≫ a the renormalon-free combination V o,L (r; a) − δΛ L (a) plays the rôle of V o,RS (r; ν f ) in Eq. (30). At r = 0 we have, V o,L (0; a) = δΛ L (a) as well as the non-perturbative equality, We redefine, to achieve formal correspondence with Eqs. (3) and (6): Note , in analogy to Eq. (45). The combination, vanishes for r = 0 and is renormalon-free. The same holds true for (53) is not only valid for r > a but also for 5 r = 0. We have, Note that the above equation is only correct up to non-perturbative O(Λ 2 QCD a 2 ) contributions to Λ L H r 0 . Again Λ L H is renormalon-free but has a power divergence. By subtracting δΛ L (a) order by order in perturbation theory one can obtain an on shell Λ OS H , but at the price of a renormalon ambiguity. Note the similarity between the above equation and Eq. (33).  The open symbols correspond to the respective gluelumps, non-perturbatively (square with pentagon) and in lattice perturbation theory (diamond and circles). 5 Based on the results of Sec. 5.2 below as well as of Ref. [21], we know that the 1 +− glueball will become lighter than the gluelump Λ L B (a) around a < r c ≈ r 0 /7, when using the Wilson action. In fact we discussed a similar situation in Sec. 2.3 above, for the Σ +′ g potential. This limit is not yet relevant for the Π u and Σ − u potentials at the lattice spacings covered in this paper. In the case a < r c , Eq. (53) will still apply for r −1 c ≫ r −1 ≫ Λ QCD , however, Eq. (51) will become modified; it would apply to the first radial excitations in the hybrid channels rather than to the ground states, until finally around a ≈ r 0 /12 a continuum of two-glueball states is encountered.
In Fig. 8 we compare non-perturbative data on the splitting between hybrid potentials with respect to the ground state potential with the perturbative expectation. The data have been obtained by us on an isotropic lattice at β = 6.2 with lattice spacing a ≈ 0.137 r 0 . Both gaps, are plotted as a function of r/a [see Eq. (10)]. The differences are indicative of the size of the expected non-perturbative O(r 2 ) contributions. We compare the non-perturbative data to the perturbative expectation for V o,L (r; a) − V s,L (r; a). The latter perturbation theory will suffer from the same renormalon ambiguity as δΛ L (a) and the difference between perturbation theory and nonperturbative data corresponds to Λ OS B . The left-most points (open symbols) correspond to the Λ B gluelump, plotted at r/a = [1/0] L ≈ 0.315.
The evaluation was done both in terms of α s (a −1 ) and To simplify the figure we disregard the uncertainty in the determination of Λ MS = 0.602(48) r −1 0 [37]. At LO and NLO lattice perturbation theory results are available [33] (diamonds and squares). Since everything is plotted as a function of r/a = [1/R] −1 L all diamonds lie exactly on top of the dashed continuous r ≫ a curves while at small distances there are differences between the dashed-dotted NLO curves and the exact NLO results (circles). In addition we plot the r ≫ a limits to NNLO (dotted curves). The shapes of the perturbative curves remain qualitatively stable while the normalization is not convergent as the order of the expansion is increased and is also strongly affected by the scale of α s (ν). This behaviour reflects the presence of the renormalon of Λ OS B , quite similar to what we can see in Fig. 4a.
By comparing with the renormalon-free right hand side (rhs) of Eq. (53) a better convergence can be achieved. However, such a comparison is only possible up to O(α 2 s ) as we do not exactly know the O(α 3 s ) contribution to the counterterm δΛ L (a) in the lattice scheme. Instead we choose to demonstrate the quality of the perturbative expansion in Fig. 9 by adding global normalization constants to all curves in such a way that agreement is produced at r/a = √ 2. (We shall return to the question of renormalon cancellation in Λ L H (a) in Sec. 5.2 below.) Indeed the differences between NNLO and NLO are smaller than those between NLO and LO. Moreover, at higher orders the scale dependence is reduced. The ν = a −1 curves seem to describe the data better at small r while the ν = ν f curves seem to work better at larger r. Up to distances as big as r = r 0 ≈ 7.3 a the perturbative curves seem to have an accuracy better than the non-perturbative uncertainties, estimated by the difference E Σ − u (r) − E Πu (r). We mentioned above that while formally the lattice spacing a −1 appears in the same places in the lattice scheme as the scale ν f did in the RS scheme of Sec. 4.3, these two scales should not be confused with each other as a −1 > r −1 > ν f > Λ QCD . Conceptionally we have been discussing the situation in which the potentials are evaluated in perturbation theory at scales ν > ν f while Λ RS H is an ultrasoft matrix elements, associated to physics at scales smaller than ν f . The lattice encapsulates the same physical picture. For instance, to each finite order in perturbation theory 6 : the power contribution to the lattice mass δm L stat (whose perturbation theory is affected by the IR renormalon of the on shell mass) corresponds to the UV behaviour of the potentials while the power contribution to Λ L H (whose perturbative expansion has the UV renormalon of Λ OS H ) is associated with the low energy behaviour of V o,L . This is the same renormalon/power term structure as in the continuum OS/RS schemes.
For a ≪ r < Λ −1 QCD lattice effects become invisible and the formulae elaborated above apply under the replacement RS → L. To illustrate this quasi-continuum limit, we eliminate the a −1 dependence from the expressions altogether, which is straight forward: where Note that the running from one scale to another is renormalon-free. For the hybrid case we can directly write, where the combination V o,L − δΛ L replaces the V o,RS of Eq. (30). Finally, we mention that the situation r = 0 on the lattice resembles the r ≪ ν −1 f continuum situation. Unlike in the continuum, however, on the lattice, even at r = 0, all observables remain finite as a −1 provides us with a hard UV cut-off.

Scale dependence
As we have mentioned in the previous sections, the running of pole mass and gluelump energies with ν f , in the RS scheme, and with a, in the lattice scheme, is renormalon-free. Therefore, the functional dependence can be described by a convergent expansion in perturbation theory. Nevertheless, in order to achieve the renormalon cancellation, the same scale ν has to be used in the perturbative expansion. This produces large logs if the scales ν f and ν ′ f are widely separated and, eventually, some errors, if one works to finite order in perturbation theory. In the RS scheme, there exists a solution to this problem. Even though δm RS (ν f ) suffers from the renormalon ambiguity, the difference δm RS (ν f ) − δm RS (ν ′ f ) is renormalon-free. We can perform a resummation of δm RS (ν f ) with any prescription to avoid the singularity in the Borel plane since it will cancel in the difference. We will take here the Principal Value (PV) prescription, which yields, where and, denotes the incomplete Γ function. The second term in Eq. (62) corresponds to Λ MS , once introduced in the sum of Eq. (61). It cancels from the combination 7 , δm PV RS (ν f ) − δm PV RS (ν ′ f ), and we will not consider it any longer. The sum of Eq. (61) represents softer and softer singularities in the Borel plane. Therefore, we expect at least the difference δm PV RS (ν f ) − δm PV RS (ν ′ f ) to converge (although, obviously, we have no mathematical proof of this). Since the first three terms are known we can check if this actually happens. We can see that this is so with a high degree of confidence in Fig. 10. We can also compare −δm PV RS (ν f )+δm PV RS (ν ′ f ) with the corresponding difference, calculated at finite order in perturbation theory: We depict this comparison in Fig. 11, where we take ν = ν f to minimise one of the logs. We see how the finite order results approach the PV curve 8 , which we will use in what follows wherever we need the running. A similar behaviour holds if, instead of δm RS , we study δΛ RS . 7 One may wonder if this cancellation materializes itself in practice since we only know the first three terms of the series. However, we checked this numerically and the results turned out to be virtually indistinguishable. 8 For finite order computations we take α s with one, two, three etc. loop running according to the order in α s at which we work. If instead, we use α s with 4-loop running (the highest accuracy known until now) the convergence to the PV result is accelerated. with non-perturbative measurements of the running against which the finite order results can be tested. It is also possible to relate results in both schemes by perturbative renormalonfree expressions. We will investigate both, the running within the lattice scheme and the translation between both schemes in Secs. 5.2 6.1 and 6.3 below.

Phenomenological analysis of the gluelump spectrum
We will determine the lowest gluelump energy, Λ B , from two different observables in two different schemes: from the non-perturbative difference E Πu (r) − E Σ + g (r) in the continuum limit in the RS scheme as well as from gluelump energies Λ L B (a) obtained at finite lattice spacings in a lattice scheme. Lattice and RS scheme can be translated into each other and we find internal consistency. We finally present results on the whole gluelump spectrum and compare our findings to previous literature.
The situation discussed here is similar to the one encountered in the "binding energy" in static-light systems which we will address in Sec. 6 below. These mesons very much resemble gluelumps, with the only difference that the source is in the fundamental representation and screened by a light quark rather than by a gluonic operator.

Determination of Λ RS B from the static potentials
We intend to determine Λ B from the hybrid potentials. For this purpose we will use our n f = 0 lattice continuum limit data on ∆E Πu (r) = E Πu (r) − E Σ + g (r) as obtained in Sec. 3. Using this difference allows us to eliminate the power divergence that appears in lattice simulations of the potentials (or, in the continuum OS scheme, the renormalon associated with the pole mass). Note that the difference has a well defined continuum limit. It is also interesting to see that the large distance linear term is cancelled as well. At the same time, Λ B will still additively contribute to this combination, see Eq. (6). In order to extract this non-perturbative constant, the perturbative difference between octet and singlet potentials has to be subtracted. For a reliable determination, the perturbative series has to be well defined and show convergence. However, this is complicated by the contribution from the renormalon discussed above and can only be achieved in a scheme where such renormalon singularities are taken into account. We have worked out the RS scheme in Sec. 4.3, which is well suited for this purpose. We fit Λ B using the following equality (see Figs. 12 and 13 for the quality of the fit): where the non-perturbatively obtained left hand side (lhs) is renormalon-free but on the rhs the renormalon can be shifted between the two contributions, the ultrasoft matrix element Λ B and the soft Wilson coefficient V o − V s , at a given order of perturbation theory. This is why we have to specify the scheme, the RS scheme in our case, which we use to eliminate (or to reduce) this ambiguity.
Note that Λ B is the only fit parameter. Also note that the above value corresponds to the n f = 0 case. The errors of this determination stem from several sources (for the above fit we use lattice data up to distances of around 0.5 r 0 ): 1) "latt." denotes the statistical error of the fit: ±0.10.
2) "th." stands for the theoretical errors. We first consider the error due to the truncation of the perturbative series (higher orders in perturbation theory/scale dependence). We obtain a first estimate by performing the perturbative expansion in α s (ν i ) or in α s (1/r). This provides us with an estimate of neglected subleading logarithms. Actually, in both cases one and the same number, Λ RS B ≈ 2.25 r −1 0 , is obtained, which we take as our central value. The effects of higher orders in perturbation theory are estimated by considering the convergence of the determination of Λ RS B at each order in perturbation theory. Working with α s (ν i ), the series {2.43, 2.37, 2.28, 2.25} is obtained. This series seems to show convergence for the last terms. In any case, the corrections are small. Working with α s (1/r), the series {2.00, 2.40, 2.31, 2.25} is obtained. This series is clearly convergent although the corrections are larger than when using α s (ν i ) as the expansion parameter. To be conservative we will take the difference between the last two terms as the error made by truncating the perturbative series: ±0.06. There is also some source of error from the normalization constant of the renormalon of the singlet and octet potential. For the singlet potential (following Ref. [15]) we estimate a 10% error in N Vs , which produces a ±0.10 error. For the octet potential, the error is very small compared with other sources of error. Even if, conservatively, we consider the general shift produced by setting δV o,2 = 0 (note that this also accounts for the error in perturbation theory of the octet static potential) our result only changes by ≈ 0.01/0.02. We will neglect this error to avoid double counting. In the above analysis we have neglected non-perturbative effects. On general grounds they have the structure The B r 2 term is due to r·O † EO type contributions in the pNRQCD Lagrangian (see Ref. [9] for details). The other term in Eq. (67) is due to r · O † ES type contributions. This produces a perturbative mass gap. F is the convolution of a short distance and a long distance piece, depending on the ratio of V RS o − V RS s over the masses of the gluelumps. For the purpose of estimating the uncertainty it seems reasonable to keep the lowest order in this expansion. This is equivalent to having a quadratic contribution, If we introduce this term into the fit, we obtain r 0 Λ B ≈ 2.30 [working with α s (1/r)] with A Πu−Σ + g ≃ −0.4 r −3 0 . We take the difference as an indication of the error due to nonperturbative effects. By summing linearly all the above errors we obtain ±0.21.
3) "Λ MS ": this error is due to the uncertainty in Λ MS = [0.602 ± 0.48] r −1 0 [37]: ±0.08. We have performed the fit using lattice data within a window of inverse distances ranging from about ν i ≈ 2.6 GeV down to ν f ≈ 1 GeV. From the plots (see Figs. 12 and 13) one can actually see that the curves follow the lattice data up to values r < ∼ r 0 . This corresponds to very low energies (< 500 MeV). Being conservative, we will not use data determined at these low energies without a better understanding of the dynamics. Nonetheless, such a fit would actually produce very similar numbers to the ones quoted above. This is even more so if a quadratic term is included. In general, introducing more lattice points reduces the statistical errors ("latt."). Including a quadratic term will reduce the theoretical error on Λ B since some of the changes that occur when altering the order of perturbation theory can be absorbed into a variation of A Πu−Σ + g . However, the addition of a second fit parameter increases the statistical error and also the uncertainty due to Λ MS . We conclude that while the individual errors depend on the precise fitting details the total error remains remarkably stable.
One might ask whether, in addition to Λ B , a reliable value of A Πu−Σ + g can be obtained. This however would require more lattice data at short distances as well as a more detailed understanding of the r 2 renormalon of the static singlet potential.
We do not consider the Σ − u data in this section as we have already established in Sec. 2.4 above that the difference with respect to the Π u potential is proportional to r 2 to leading order. Hence we cannot obtain any independent new information on Λ B from these data, that have larger statistical errors.

Determination of Λ RS
There exists a direct determination of Λ L B (a) (the 1 +− or B gluelump) by Foster and Michael [12]. The numerical values are displayed in Table 3, where we used the same r 0 /a values as were used in this reference. It is clear from the discussion in Sec. 4.4 that these are perfectly sensible numbers if incorporated into a global scheme with renormalon cancellation, for instance, with the potentials also defined in the lattice scheme as in Sec. 4.4. In doing this we are able to independently determine Λ B in a different scheme. Consistency would require that after translating the lattice into the RS scheme the results should agree with each other. We will check this in this section.  (17) 3.55 (10) 3.25 (13) 3.16(13) Table 3: The inverse lattice spacing, the mass of the 1 +− gluelump Λ L B in the lattice scheme, as well as its conversion to the RS scheme to different orders in perturbation theory. NNNLO* stands for an estimate obtained neglecting 1/N 2 c corrections, for details see the text. In the last column, we state the values of Λ RS B (a −1 ) using Eq. (66) and the running according to the PV prescription, Eq. (61). The errors only incorporate the statistical uncertainties as well as the 8 % uncertainty in Λ M S r 0 , added in quadrature, but no estimates of "theoretical" errors.
The master formula that relates the lattice and the RS scheme reads (known up to NNLO), Both, Λ RS H and Λ L H have a power-like dependency on ν f and a −1 , respectively, but are renormalon-free, Λ L H exactly and Λ RS H within the precision of our estimation of the renormalon contribution. This implies that the combination δΛ L + δΛ RS does not contain a renormalon either if calculated in a consistent way: δΛ L (a) and δΛ RS (ν f ) contain one and the very same renormalon contribution (with negative relative sign). The sum of both terms, expanded in terms of α s has good convergence properties (using the same normalization point to enforce the renormalon cancellation at each order in perturbation theory). The explicit expression at NNLO reads, where the v i can be found in Eqs. (38) and ( Table 2. An estimate of the O(α 3 s ) term can be obtained from Eq. (78) below, under the replacements, C f → C A andṼ s,i → 2(Ṽ s,i −Ṽ o,i ). This estimate will be subject to O(1/N 2 c ) corrections to the coefficient v 3 .
In principle, ν f and a −1 need not be equal but we will take them similar to avoid large logs. The large numerical values of v 2 and b 1 are mainly due to contributions from lattice-specific tadpole diagrams that arise because the breaking of Lorentz symmetry becomes particularly evident at UV scales ≃ a −1 . This often results in badly convergent perturbative series when expanded in terms of α L (a). However, the convergence is vastly improved, once the series is re-expressed in terms of a more "physical" coupling like α s (a −1 ) = α L [1−b 1 α 2 L /(4π)+· · ·] (see e.g. Refs. [33,38]). This is also evident from Eq. (70) above as v 1 ≈ 3.17,  We can now translate the Λ L B values obtained by Foster and Michael [12] into the RS scheme. The results are shown in Table 3 and are also displayed in Fig. 14. "NLO" and "NNLO" refer to translating from the lattice scheme to the RS scheme via Eq. (70) to O(α s ) and O(α 2 s ), respectively 9 . Obviously, to leading order, Λ B is scheme independent. "NNNLO*" stand for an estimate obtained assuming that the NNNLO contribution to δΛ L is equal to the NNNLO contribution to δm L with the replacement of the overall factor C f → C A . This is correct up to O(1/N 2 c ) effects. Finally, the conversion from the lattice to the RS scheme has been performed using the 4-loop running of α s at ν = a −1 = ν f . This accelerates the convergence to the RS results. If, instead, we use the n-loop running of α s that is consistent with the order of the calculation, we still see convergence but with, in the NLO and NNLO cases, larger corrections. This is mainly due to the fact that within the present window of energies the values obtained for α s (ν) from Λ MS from a one-or two-loop running are significantly different from those from the three-loop running (which is close to four-loop). The lattice prediction of Λ MS that we use as an input applies to very high energies, such that it is important to run α s down to ν ≥ 2.5 r −1 0 as precisely as possible. Within present errors we can fit the data with straight lines but there will be logarithmic corrections and, in the gluelump data Λ L B r 0 , additional O(a 2 ) = O(ν −2 f ) lattice artifacts. The figure reveals that at the lattice spacings investigated these are tiny, relative to the linear slope. Except for these lattice corrections the running of Λ L B is non-perturbatively accurate. Needless to say that the power dependence on a −1 is universal for all gluelumps, such that gluelump mass splittings have a well defined continuum limit, which is also confirmed in Ref. [12].
In lattice perturbation theory we can calculate the "running" of the gluelump data to O(α  Table 3 and perform the running with NNNLO* accuracy, we obtain the dashed line that joins the "LO" RS(= L) points. We can see that this parametrization is quite close to the non-perturbatively evaluated data. Moreover, there is overall convergence, with higher order terms being numerically smaller in the lattice scheme. We will discuss this in more detail in Sec. 6 below, in the context of the static-light binding energy Λ L , which has a similar perturbative expansion, up to an overall factor C f /C A , see , in fact we already find agreement at the NNLO level. In Secs. 6.1 and in particular 6.3 below we will analyse the running of the binding energy of static-light mesons in more detail, see also Fig. 18.
We obtain an independent second prediction for Λ RS B from the gluelump data. The statistical errors are smaller in the gluelump case than those we encountered from the continuum potentials. In a first step we obtain the fit parameter, from a global NNNLO* fit, where we have chosen ν = ν f = 7.32 r −1 0 . We can then convert this result into, using the PV running in the RS scheme. This compares well with the result from the potentials, Eq. (66). The errors displayed in Eq. (71) above are due to the following sources: 1) "latt." is the sum of the statistical error (±0.03) and the error encountered when varying the fit range (i.e. excluding the left-most data point): ±0.01.
2) "th." is the sum of perturbative and non-perturbative errors. As perturbative errors we take the difference between NNLO and NNNLO* results (±0.20) as well as a 10 % uncertainty in N Vs −N Vo (±0.18). To investigate possible non-perturbative effects we include an a 2 term into the fit. We estimate an additional ±0.04 uncertainty from this source. Adding these three errors linearly results in ±0.42.
Whereas the statistical error is smaller in this determination than the one of Eq. (66) and the uncertainty due to the error of Λ MS is comparable in size, the systematics are less well under control, which is reflected in the large theoretical error. First of all, for the lattice gluelumps we only have the perturbative result to O(α 2 s ) with an estimate of the O(α 3 s ) term while in Sec. 5.1 above we knew the O(α 3 s ) results and have an estimate of the O(α 4 s ) terms. Furthermore, as the previous analysis was based on observables with a well defined continuum limit, we circumvented the problem of disentangling the a −1 "running" of Λ B r 0 from O(Λ 2 QCD a 2 ) lattice artifacts. With gluelump data on more, and in particular finer, lattice spacings the latter disadvantage (which at present is however not the dominant one) can in principle be overcome. In conclusion, it is nice to observe perfect agreement between the two predictions, which enhances our confidence in the methods applied and adds further credibility to our error estimates.

Higher gluelump excitations
Now that we have fixed the energy of the lightest gluelump, we can quote absolute values for the remaining gluelump spectrum using the results of Foster and Michael [12]. We display our predictions in Table 4 where the errors correspond to the sum of the individual uncertainties, added linearly. The dominant uncertainty is that of Λ B , as the mass differences between the different gluelumps have been determined with very good accuracy. Needless to say that these results are scheme and scale dependent. The quoted numbers refer to the RS scheme with ν f = 2.5 r −1 0 ≈ 1 GeV. With the information presented in this paper they can be run to different scales. For ease of reference we also converted these values into GeV units (using r −1 0 = 394 MeV). However, we note that one should add a scale uncertainty of about 10 % to them to account for the fact that all results have only been obtained in the quenched approximation.
Note that the gluelump operators can be represented in terms of gluonic fields [9,39]. In general one and the same gluelump can be created by infinitely many different adjoint operators H. Within each channel we display (one of) the lowest dimensional such choice(s) in the table. The basic building blocks are the covariant derivative D i (with J P C = 1 −+ ,  Table 4: Absolute values for the gluelump masses in the continuum limit in the RS scheme at ν f = 2.5 r −1 0 ≈ 1 GeV, in r 0 units and in GeV. Note that an additional uncertainty of about 10 % should be added to the last column to account for the quenched approximation. We also display examples of creation operators H for these states. The curly braces denote complete symmetrization of the indices. dimension 1), the chromomagnetic field B i (1 +− , dimension 2) and the chromoelectric field E i (1 −− , dimension 2). The curl of the electric field has the quantum numbers of the magnetic field, such that on the lattice all states can be created by operators that are local in time. Furthermore, D · B and D · E can be eliminated, the first because it is identically zero, using the Jacobi identity, the second by applying the equations of motion. One example: the lowest dimensional operator that creates the 3 +− state is D {i D j B k} , where the curly braces denote the sum over all 10 symmetric permutations of the indices. This includes three terms D i j D j B j = 0 such that indeed there remain only seven independent operators to create this seven dimensional representation. Also note that D {i B j} and D {i E j} each only contain five independent operators, consistent with J = 2 etc..
It is interesting to see that the level ordering roughly corresponds to the lowest dimension of the creation operator, once the equations of motion are used to eliminate the E field [9]. This makes the E field "heavier" than a B field, increasing its dimension by one. The 3 −− gluelump (two derivatives and one E which corresponds to dimension five, after substituting E) is not included into the table as no controlled continuum limit extrapolation was possible. However, its mass at fixed finite lattice spacing is in the same ball park as that of the other dimension five states, 4 −− and 1 −+ , in support of this naïve operator counting picture.

Comparison with previous results
We shall relate our results to previous determinations of the gluelump masses. All these suffer from the problem of obtaining the global constant and, in none of these, the scheme was clearly defined, such that they need not yield the same results that we obtain.
In Ref. [39] the gluelumps were studied within a string model. One general feature of this approach is the excess of predicted states. This seems to be a problem of this model since it does not appear to be compatible with QCD, or more precisely with its realization for this kinematical regime: pNRQCD [9] (see also the discussion in Ref. [40]). The prediction of this model, Λ B (n f = 0) = 1.87 GeV, is by a factor of two larger than our result. In Ref. [41], the same value for the electric and magnetic correlation length is obtained: Λ E (n f = 0) = Λ B (n f = 0) = 0.90(5)(10) GeV, from lattice simulations using the cooling method. The number for Λ B coincides with ours. However, the splitting between chromoelectric and -magnetic correlators is unaccounted for. From the results of Foster and Michael one would then assign a systematic error of the order of this splitting ≈ 400 MeV: clearly a better conceptional understanding of how "cooling" removes short distance fluctuations, without destroying essential infrared physics would be useful. On the other hand it is comforting that numbers similar to our results are obtained in this approach, that is also meant to subtract the perturbative contributions from the low energy matrix element.
In Ref. [42] a sum rule analysis of the electric and magnetic correlator was made. The main result was Λ E (n f = 0) = (1.9 ± 0.5) GeV. It should be noted that the value of Λ MS on the lattice is now smaller by 5%, compared to the value used in this analysis. Taking this into account we find this result compatible with ours [1.25(16) GeV], within errors. Moreover, in this analysis, evidence for Λ E > Λ B was reported.
In Ref. [43], an MIT bag model calculation was used to obtain the gluelump spectrum. No errors were assigned to this evaluation. The value of Λ B is by about 500 MeV larger than ours and quite consistent with the sum rules evaluation. The same holds true for Λ E , however, for the higher excitations the agreement with the results of Foster and Michael is less convincing.
In Ref. [11], lattice correlation functions that are needed to calculate relativistic corrections to the static potential were used in order to check the validity of the stochastic vacuum model in the Gaussian approximation. Under this assumption, which was to some extent tested in this reference, these correlation functions could be related to gluonic field strength correlators and upper limits for the gluelump masses were obtained: Λ B (n f = 0) ≤ 1.64 (16) GeV and Λ E (n f = 0) ≤ 1.04(15) GeV, respectively: the ordering of the gluelumps is wrong, however, the upper limits quoted are in no contradiction to our results (or indeed to a different ordering).
In Ref. [44], a constituent quark model was used. The results roughly agree (within a 200 -300 MeV error) with the splittings predicted by Michael and Foster and the hybrid spectrum at short distances (see Ref. [40] for some criticism of this evaluation). For the lightest gluelump they obtain Λ B ≈ 1.4 GeV.
We have seen how different determinations of Λ B result in values ranging from less than one GeV up to nearly two GeV. These numbers are all scheme dependent. This may explain the huge differences between different results. Our result provides strong constraints on vacuum models. Furthermore, the RS scheme provides a unified framework to study the non-perturbative effects in an unambiguous and model independent way.

Static-light systems
The situation discussed above very much resembles the one that one encounters in heavylight mesons in the static limit. In this case, the adjoint source is replaced by a fundamental source which is not screened by gluonic fields but by a light Dirac quark instead. (A light Higgs scalar in the fundamental representation would be an alternative possibility.) In these systems the binding energy Λ of the 1 vector heavy-light mesons, once 1/m b corrections and the spin of the heavy quark are taken into account) plays a rôle similar to that of the Λ B discussed above. The experimental mass of the B meson M B can be factorized into, where both Λ and m b depend on scheme and scale. In the literature (see e.g. Ref. [18]) the binding energy in the lattice scheme is referred to as E(a) = Λ L (a), which is renormalon-free but has an a −1 power divergence. For the Wilson action and n f = 0 this δm L stat (a) power term is known to O(α 3 s ) in perturbation theory [Eqs. (37) - (40)]. Subtracting this perturbative result introduces renormalons.
It is also possible to define the binding energy in an entirely non-perturbative renormalonfree and power-term free way, for instance by subtracting the energy of a temporal Schwinger line in Coulomb or Landau gauge [45]. In fact the same can be achieved in the case of the lowest gluelump mass, either by subtracting the energy of an adjoint Schwinger-line in a fixed gauge (see also Ref. [19]) or by subtracting the on-shell mass of an adjoint Polyakov-Wilson line, encircling a compactified lattice dimension. From an EFT point of view however one would like to combine a non-perturbative low energy result with a perturbative calculation at high energies. For instance to quote a value for the b quark mass in the MS scheme, the UV renormalon of the binding energy is required to cancel the IR renormalon of the OS mass and hence a perturbative subtraction is essential: the renormalon of the expansion of the power divergence is the same as the one that is encountered in the conversion from the OS mass into the MS mass. This procedure has been implemented in the past in calculations of the b quark mass from lattice simulations in the static limit [18].
The b quark mass has also been obtained in perturbative QCD in the RS scheme at ν f = 2 GeV from the Υ(1S) system using EFTs [15]. Subtracting this value from the spin-averaged mass of the B meson yields, This number is different from the value quoted in Ref. [15] 10 , since here we have performed the running to ν f = 1 GeV using the PV prescription and not included O(1/m b ) corrections into the fit (these two effects partially compensate each other). Using the PV prescription allows us to perform the log resummation for the renormalon related terms. However, the result strongly depends on the value of Λ MS . Eq. (75) has been obtained from the physical Υ and B systems, not in the quenched approximation. The scale r −1 0 = 394 ± 20 MeV [20,46,47] is also obtained from Υ phenomenology. Re-expressed in terms of r 0 we get, In what follows we will extract Λ RS from lattice data of static-light mesons. After addressing the b quark mass we will conclude with a more detailed study of the running in the lattice and RS schemes, using precision data from the static potential within an energy range, 2 < ∼ r 0 ν f = r 0 /a < ∼ 15.

RS
We will use Eq. (76) as our starting point for the n f = 0 situation. In order to compare with lattice results in the quenched approximation we will employ the n f = 0 running of Λ RS (ν f ) and keep in mind that on top of the errors stated above one might expect an additional 10 % quenching error.   [51], the static-light binding energy Λ L = E [48,34,49,50] in the lattice scheme, as well as its conversion to the RS scheme to different orders in perturbation theory. In the last column, we state the values of Λ RS B (a −1 ) using the PV running, Eq. (61), of the result Eq. (76) in the RS scheme. The errors only incorporate the statistical uncertainties as well as the 8 % uncertainty in Λ M S r 0 [37], added in quadrature. The values in the last column, which have been obtained from the physical Υ(1S) and B meson masses, have additional errors inherited from Eq. (76), which, however, will only result in an overall upward or downward shift and will not affect their differences.
Λ L (a) has been calculated on a variety of lattice spacings by different collaborations [48,34,49,50]. The main source of uncertainty in these determinations is the extrapolation to zero light quark mass. We used the r 0 /a values from the interpolation of Ref. [51] to assign the scale 11 . The results are displayed in Table 5 and are roughly consistent with each other, with the exception of the coarsest lattice point r 0 ≈ 2.93 a that corresponds to β = 5.7. Here the raw data of Ref. [50] are more accurate but the chiral extrapolation of Ref. [34] should be better controlled. We multiply the values obtained for Λ L B r 0 of Ref. [12] (that are displayed in Table 3) by the colour factor C f /C A . At β = 5.7, 6.0 and 6.2, respectively, we obtain the numerical values 2.37(4), 3.11(2) and 3.72 (2). The corresponding values in Table 5 read 2.45(6)|2.22(4), 3.28(6) and 3.83(8)|3.87 (11) where for both, β = 5.7 and β = 6.2, two independent determinations exist. The qualitative agreement is remarkable: not only the perturbative expansions of δΛ and δm stat are dominated by terms that are proportional to the respective Casimirs of the gauge group representation of the static source but also the non-perturbative values  Figure 15: Perturbative running of the binding energy Λ in the lattice scheme, in comparison with lattice data, starting at the smallest available lattice spacing. The NNNLO error band incorporates the error due to the uncertainty in Λ MS [37], and the statistical error.
themselves. In fact also in the RS scheme the result Eq. (66) is close to the value displayed in Eq. (76), multiplied by C A /C f = 9/4. Similar to the discussion in Sec. 5.2 above, we can translate the results from the lattice scheme into the RS scheme. The master formula in this case is very similar to Eq. (69) and reads (known to NNNLO), with where,  and the coefficientsṼ s,1 andṼ s,2 can be found in Table 2. The coefficients v i and b i can be found in Eqs. (38) - (40) and Eqs. (42) and (43), respectively.
Eqs. (77) and (78) also relate results obtained at different lattice spacings to each other, To illustrate this we display the Λ L (a) values of Table 5 in Fig. 15, together with the expected running, starting at the finest, i.e. right-most, lattice point at LO, NLO, NNLO and NNNLO. The NNNLO error band contains both, the statistical error and that due to the uncertainty in α s (a). The running is done in each order in a self-consistent way to the given order in α s , according to Eq. (78) (without theṼ s,i terms). We used ν = 9.76 r −1 0 and the initial value, α s (ν), was calculated from Λ MS r 0 using the four loop running. We observe convergence and moreover the series is sign alternating. To NNNLO, except for the lower lying of the two r 0 /a ≈ 2.93 data points, there is no contradiction between data and the expectation. However, the points of Ref. [34] have a slightly more pronounced slope such that the r 0 /a ≈ 2.9, 4.5 and 6.3 points (β = 5.7, 5.9 and 6.1) lie below the curve while the r 0 /a ≈ 8.5 point (β = 6.3) lies somewhat above.
We also display the data of the table in Fig. 16, in analogy to Fig. 14 but disregard the LO result that is already displayed in Fig. 15. The size of δm L stat (a)−δm RS (a −1 ) increases linearly in a −1 , with logarithmic corrections: at coarse lattice spacings there might be significant perturbative O(α 4 s ) and non-perturbative O(a 2 /r 2 0 ) corrections affecting the slope of this function while at fine lattice spacings the slope can be determined accurately but the δmdifference itself becomes large. An accurate conversion between the two schemes can therefore neither be obtained at extremely fine nor at very coarse lattice spacings. Setting ν = ν f , the difference between NLO and NNLO translation is minimised for 3 < ∼ r 0 ν f < ∼ 4 while that between NNLO and NNNLO is minimal for 7.5 < ∼ r 0 ν f < ∼ 9, where the widths of these bands are determined by our uncertainty in the value of Λ MS r 0 .
We choose to translate the lattice scheme results into the RS scheme by means of a global NNNLO fit to the r 0 /a > 5, i.e. β ≥ 6.0 data, expanded in terms of α s (ν = 9.76 r −1 0 ), where we set ν f = ν. The result reads, Note that Λ RS is the only fit parameter.
The dashed curves in Fig. 16 correspond to such an NNNLO fit to the LO results, subsequently transformed in the same way as the data points to NLO, NNLO and NNNLO. The error band corresponds to the result Eq. (81) above, without the theoretical error, run to different energies, using the PV prescription, Eq. (61): unlike the band displayed in Fig. 14 above, this is not the result of an independent determination. We would also have found agreement with the result Eq. (76), but only within the large theoretical errors of this un-quenched determination.
At high order in the perturbative expansion and at high energies one would expect the slope of the non-perturbative running in the lattice scheme, translated into the RS scheme, to approach that of the running within the RS scheme. Discarding the four data points of Ref. [34], Fig. 16 nicely confirms this expectation. We will investigate this running with higher accuracy in Sec. 6.3 below.
The errors of the determination Eq. (81) above stem from the following sources: 1) "latt." is the sum of the statistical error (±0.03) and the error encountered when varying the fit range a −1 min r 0 = 4.48, 5.37, 6.32 (±0.05): ±0.08. 2) "th." is the sum of perturbative and non-perturbative errors. As perturbative error we take the difference between NNLO and NNNLO results. Varying the fit range as above this difference never exceeds ±0.04. We also study the error due to the uncertainty of N Vs obtaining ±0.06. To investigate possible non-perturbative effects we include an a 2 term into the fit. We estimate an additional ±0.08 uncertainty. Adding these three errors linearly results in ±0. 18.
from the value Eq. (81). This n f = 0 result compares reasonably well with the phenomenological n f = 4 value of Eq. (76) above and its error is of a comparable size.

Comment on the b quark mass
We cannot resist the temptation to obtain a value for the RS scheme bottom quark mass, using Eq. (74) and our quenched result Eq. (82). We obtain, We feel that 1 GeV might be a more natural scale to obtain an n f = 4 prediction from the quenched model than 4 GeV and hence prefer the central value of Eq. (84). After all, the quenched model has been adjusted to reproduce low energy QCD phenomenology and indeed Eqs. (76) and (82) agree with each other within errors. However, as discussed above and as indicated by the 150 MeV difference from using a different perturbative running, such predictions have to be consumed with some caution. Eq. (84) demonstrates the precision that can be achieved in lattice simulations of static-light mesons with sea quarks to NNNLO. Obviously, the "latt." error can systematically be reduced. Note that, with NNNLO perturbative results, the dominant theoretical uncertainty (apart from the sea quark content) is due to 1/m b corrections.

The running of Λ from the static lattice potential
To leading order, the singlet static energy E s is the sum of twice the heavy quark mass and the singlet potential, Eq. (2), while M B is the sum of the quark mass and the binding energy Λ, Eq. (7). Consequently, in the OS (RS) schemes V s contains twice the leading renormalon (power term) of Λ. In QCD with sea quarks this is also evident from the large distance behaviour, where E s (r) will approach 2M B .
In the lattice scheme, the non-perturbative energy E L Σ + g differs from E s by twice the quark mass, Eq. (45), and contains the same power term as the static-light energy Λ L (times two).
One can explicitly verify this in perturbation theory. In QCD with sea quarks E L Σ + g (r) will approach 2Λ L for r > ∼ r c , where r c denotes the distance associated with "string breaking" and is implicitly defined by, E s (r c ) = 2M B . We find the static potential [52,46] E Σ + g to exceed the values of 2Λ L of Ref. [48] at r > r c = (2.25 ± 0.15) r 0 , a distance that is statistically indistinguishable from the value r c ≈ 2.3 r 0 , obtained in simulations with n f = 2 light sea quarks [3,13].   (87), as well as its conversion to the RS scheme to different orders in perturbation theory. The errors only incorporate the statistical uncertainties of the E Σ + g (r 0 ) data, as well as the 8 % uncertainty in Λ M S r 0 [37], added in quadrature. The overall error due to the uncertainty in ∆, which does not affect the running of Λ L , is not displayed.
In what follows we will investigate the running of, as a function of a −1 . The static lattice potential E Σ + g (r 0 )/2 can be determined more precisely than Λ L : in terms of computer time it is cheaper to obtain with the same statistical error and, since no chiral extrapolation is involved, with virtually no systematic uncertainties. However, we do not know the absolute normalization ∆. We re-analyse the lattice potentials of Refs. [52,46], to correctly account for the propagation of the uncertainty of r 0 into the combination r 0 E Σ + g (r 0 ). By matching the lattice potential E L Σ + g (r 0 )/2 = (2.856 ± 0.014) r −1 0 to Λ L = (3.844 ± 0.065) r −1 0 at β = 6.2 (r 0 /a ≈ 7.3), where we have two independent results for the latter quantity [48,49], we obtain For ease of comparison with Sec. 6.1, we display the resulting Λ L pot (a) in Table 6 as well as in the figures. The additional uncertainty due to the error in ∆ should be kept in mind. We display Λ L pot in Table 6, together with conversions into the RS scheme, according to Eqs. (77) and (78). The data are also depicted in Fig. 17 (full diamonds), together with the results from the static-light energies Λ L (open diamonds). Except for the four data points of Ref. [34] at r 0 /a ≈ 2.9, 4.5, 6.3 and 8.5, whose slope is somewhat incompatible with the results from the other references as well as with perturbation theory (as we already noticed in Sec. 6.1 above), we find agreement between the non-perturbative running of E Σ + g (r 0 )/2 and that of Λ L , down to the lowest scales. This need not be so since in principle the results may differ by O(a 2 /r 2 0 ) lattice terms. We also compare this running with the expectation from the value Λ MS ≈ 0.602 r −1 0 (solid line), where we use the normalization suggested by the central value of Eq. (81), Λ RS (ν f = 9.76 r −1 0 ) ≈ 1.70 r −1 0 . As can be seen there is no contradiction between the lattice data and NNNLO perturbation theory down to scales as low as 2 r −1 0 and as high as 15 r −1 0 . This agreement is quantifiable: a one-parameter NNNLO fit to the a −1 > 5 r Including all available data results in χ 2 /N DF = 6.91/9 with Λ RS (9.76 r −1 0 ) = (1.70 ± 0.01) r −1 0 . The errors of the above examples are purely statistical. The uncertainties in ∆ and Λ MS as well as theoretical errors are unaccounted for. If we go to NNLO we obtain the χ 2 /N DF values of 16.3 (all data points), 23.0 (a −1 > 5 r −1 0 ) and 6.7 (a −1 > 9 r −1 0 ). Also, the predicted value of Λ RS (9.76 r −1 0 )r 0 becomes somewhat unstable, ranging from 1.76 (all data points), 1.79 (a −1 > 5 r −1 0 ), up to 1.91 (a −1 > 9 r −1 0 ): within the accuracy of the data it is essential to go to at least NNNLO in perturbation theory.
In Fig. 17 and Table 6 we have also displayed the results, translated into the RS scheme to different orders in perturbation theory. In Fig. 18 we focus on this comparison. This figure very much resembles Fig. 16, only that now the error bars are smaller as we discard the error of ∆, which will only affect the overall value of Λ but not the running with the scale. The dashed lines correspond to NNNLO perturbation theory in the lattice scheme with Λ MS = 0.602 r −1 0 , and the central value of Eq. (88) as normalization point. This running perfectly agrees with the data down to very low energies. As already observed in Sec. 6.1 above, we also find nice convergence for a −1 > ∼ 3 r −1 0 , as the order of the perturbation theory is increased. The error band corresponds to the PV prescription of the running in the RS scheme 13 with Λ RS (9.76 r −1 0 ) = (1.70 ± 0.04) r −1 0 , run to different scales, using Λ MS = (0.602 ± 0.048) r −1 0 . Note that the errors that we display in this case are only due to the uncertainty in α s , with all other error sources of Eq. (81) (as well as the uncertainty of ∆) ignored.
We find excellent agreement between data and the predicted running. In fact, one can in principle determine α s from the logarithmic corrections to the a −1 running of the binding energy: in dedicated lattice simulations of the short distance static potential tremendous statistical accuracy can be achieved and tiny lattice spacings are accessible [25]. Even using our static singlet potentials [52,46] that are less accurate than those of this recent reference, a two-parameter NNNLO fit to the a −1 > 5 r −1 0 data yields, Λ MS = (0.590 ± 0.036) r −1 0 13 At ν f ≫ ν = 9.76 r −1 0 we find some differences between the NNNLO running in the lattice scheme (dashed black line) and the PV prediction (error band), due to large logs in the difference Eq. (80), where we have not attempted a log resummation.  In conclusion, we have demonstrated that the running of the binding energy in the lattice scheme can be reproduced with incredible accuracy in NNNLO perturbation theory, in terms of α s . This accuracy is possible since, unlike in the case of the binding energy itself, there is no leading renormalon contribution to its running. Down to energies of about 1 GeV we do not see any sign of a break-down of perturbation theory or evidence of significant non-perturbative contributions to the running. We have also confirmed that the theoretical errors estimated in Eqs. (81) and (82) are indeed conservative.

Gluinonium and other related issues
We already mentioned that gluelumps are interesting in the context of bound states including heavy adjoint particles, such as gluinos of SUSY models (even if it is quite likely that they will decay before any kind of hadronization takes place). In this case, to leading order in HGET (Heavy Gluino Effective Theory), the gluino mass can be obtained from the relation, in a scheme of choice that can then be converted into the mass in say the MS scheme mg(mg), analogously to the discussion of Sec. 6 above. We will limit most of the discussion below to the RS and OS schemes but translation into lattice schemes is straight forward.
MG denotes the mass of the lightest (spin-averaged) glueballino. Note that in this context the gluelump energy Λ B plays the same rôle as the binding energy Λ did for heavy-light mesons. We have Λ RS δmg ,RS = −δΛ RS in the glueballino case corresponds to the δm RS of heavy-light mesons. We can also write down the above equations in the lattice scheme in which case, using the same conventions as in other parts of this paper, δmg ,L = δΛ L .
In addition to glueballinos one can imagine gluino-gluino bound states: gluinonia, Γ. Their dynamics is dictated by the following Lagrangian, at leading order in 1/mg and in the multipole expansion. This is analogous to Eq. (1), replacing the static sources in the fundamental by static sources in the adjoint representation. This means that there will be further multiplets in Eq. (91) that we will not consider in this paper. The singlet potential between two adjoint sources V A,s (r) has been calculated in perturbation theory to O(α 2 s ) and the corresponding energy E L A (r; a) was determined in lattice simulations (see e.g. Ref. [24]). Up to lattice artifacts ∝ a 2 we can write, where the normalization of E L A (r) with respect to E A (r) can be obtained from the gluinonium spectrum. Obviously, lim while for the bottomonium energy in QCD with sea quarks one obtains (up to 1/m b corrections and neglecting radial and gluonic excitations), In combining Eq. (94) with Eqs. (26) and (30) we obtain the important equality, where we have used the fact that MG = mg ,RS + Λ RS B and E B ∈ {E Πu , E Σ − u }. The effect of δm RS cancels from the combination E B − E s and δΛ RS from E A + 2E H . Since the glueballino mass is a physical observable this implies that, up to O(r 2 ) corrections, the combination V A,s (r) + 2[V o (r) − V s (r)] is scale independent and free of renormalon and power contributions: the UV renormalon of V o is cancelled by the UV renormalon of V s while the leading IR renormalon of V o , which we studied in this paper, is cancelled by one half of the UV renormalon of V A,s . In fact, to O(α 2 s ) this combination explicitly vanishes and the O(α 3 s ) term is suppressed by a colour factor 1/N 2 c . In the above equation E B (r) corresponds to the Π u or Σ − u hybrid levels. For a general E H (r) we would have to replace the MG on the rhs by the mass of the excited glueballino in the respective channel. At r → ∞ the rhs of Eq. (97) will approach 2MG, see Eqs. (95) and (96).
We wish to compare our expectation with lattice data. This can either be done after an extrapolation of these to the continuum limit or at finite lattice spacing in the lattice scheme. Re-expressing Eq. (97) in terms of the energy levels as determined in the lattice scheme [E A (r) = E L A (r; a) + 2mg ,L (a) etc.], and using the conventions of Sec. 4.4 above, this amounts to, Both lhs and rhs of the above equation are explicitly free of a −1 power terms (and of leading renormalons). In fact the rhs vanishes in perturbation theory, to at least O(α 3 s /N 2 c ). As indicated in the equation, in general there will be non-perturbative O(a 2 /r 2 ) as well as O(Λ 2 QCD a 2 ) lattice artifacts, in addition to the O(r 2 ) corrections from higher terms in the multipole expansion.
The above combination is extremely interesting as for small r there should only be a quadratic but no linear term. At r > ∼ 2 r 0 the adjoint string will break and the lhs of the equation will approach zero like 1/r. In the intermediate region 0.5 r 0 < r < 2 r 0 one would expect two non-perturbative contributions, a linear term from the slope of E A (r), with an effective string tension [24], σ eff = [3.09 ± 0.10]r −2 0 , as well as a 1/r term that dominantly originates from the difference between static hybrid and singlet potentials and whose coefficient will approach 2π as r → ∞, in an effective string model expectation. In fact for r ≈ r 0 one would expect this 1/r term still to dominate over the linear term.
We wish to compare this expectation to numerical data. Unfortunately, on isotropic lattices where we know the gluelump mass in the lattice scheme we did not compute the adjoint potential while on our anisotropic data sets all potentials, singlet, adjoint and hybrid are available but the gluelump mass is unknown. In Fig. 19 we display the combination continuum limit, see Sec. 3 and Ref. [24]. Note that there is an additional 1 % overall error on ordinate and abscissa due to the conversion from lattice units into units of r 0 that we do not display.
From Eq. (98) we would expect the combination shown to approach the gluelump energy in the lattice scheme, Λ L B (a), as r → 0. We see that the approach towards this limit is remarkably flat. In fact, excluding the r > 0.9 r 0 data, which are clearly in the nonperturbative regime anyway, we are unable to resolve deviations of the data from a constant. Note that the units on the ordinate, 0.2 r −1 0 ≈ 80 MeV, are quite small. A linear plus quadratic fit, to r < 0.5 r 0 data yields, A purely phenomenological fit to the same functional form for all distances results in, while in a physically completely unmotivated funnel parametrisation, 2Λ L B + e/r + k r, we obtain, the r dependence is so weak that on the 1 % error level of the lattice data we are unable to discriminate between different parametrizations. However, we can determine the gluelump mass rather precisely, Λ L B (a σ , a τ ) = (7.75±0.05±0.07) r −1 0 , where the second error reflects the uncertainty in the lattice determination of r 0 /a τ . In fact the same can be done for the a σ ≈ 0.23 r 0 and a σ ≈ 0.33 r 0 data sets. The respective results read, Λ L B = (6.71 ± 0.04 ± 0.09) r −1 0 and Λ L B = (5.75 ± 0.10 ± 0.05) r −1 0 , respectively. The data are in agreement with a linear slope in a −1 but, unfortunately, at present we only know the NLO perturbation theory for the anisotropic case. After subtracting twice these gluelump energies, we find scaling of the coarse lattice data with the results depicted in the figure, within error bars of comparable size.
In particular from the fit to the funnel type parametrisation we see that the data leave little room for perturbation theory style short-distance Coulomb terms. This is in agreement with our expectation. However, miraculously there is also no evidence for a quadratic term in the r < 0.9 r 0 data and in fact we can set the limits 0.46 > c r 3 0 > −0.18 for such a contribution, from the r < 0.7 r 0 data. We believe that the smallness of this term is accidental as had we replaced the Π u by the Σ − u potential, it would certainly be present, see Sec. 2.4. One can however speculate that there might be a cancellation of r 2 effects and that Π u does not receive a significant r 2 contribution in the multipole expansion. This issue should be addressed in future theoretical and numerical studies with enhanced accuracy.
The observed slope at larger distances (≈ 0.84 r −2 0 ) is much smaller than that of the adjoint potential in this region (≈ 3.09 r −2 0 ), in agreement with our expectation that the 1/r contribution from the difference 2(E Πu − E Σ − u ) cannot be neglected. There is no evidence of a linear non-standard short-distance term for r < 0.9 r 0 , at least not of the size expected in various models [53]. A possible explanation of the absence of such a term from our calculation of a quantity that vanishes to low orders in perturbation theory would be that α s (q) itself receives O(1/q 2 ) corrections (see Refs. [53]). We remark to this end that α s is not a physical observable. In the MS scheme it is perturbatively defined. The difference between α s and any non-perturbative generalisation of this coupling, that would allow inclusion of such singularities, will necessarily not be universal but depend on the prescription used. However, we are investigating a physical observable here that is scheme independent.
Combinations of different potentials that lead to renormalon and low order perturbation theory cancellations are certainly an arena worthwhile to study for a determination of higher order terms in the multipole expansion and for testing the validity of the standard operator product expansion picture. As we shall detail below many such combinations exist.
There are also hybrid excitations in the adjoint channel. The perturbation theory in this case is richer than for potentials between fundamental sources as 8 ⊗ 8 = 1 ⊕ 8 ⊕ 8 ⊕ 10 ⊕ 10 * ⊕ 27: in addition to singlet and octet, we have another octet, two decuplet fields and a 27-let which have to be included into Eq. (91). Consequently, adjoint hybrid potentials cannot only have the octet perturbative expansion but some will correspond to decuplets and others to 27-lets. Note that the decuplet representation is not self-adjoint but has vanishing triality as it should be.
The renormalon of the octet potential between adjoint sources is the same as in the fundamental case but the decuplet and 27-plet adjoint potentials contain new renormalons, which are related to those of the singlet potentials between colour charges in these respective representations. This exactly resembles the situation discussed above where the adjoint singlet potential contains the same renormalon as the fundamental octet potential. In fact one can define an infinite tower of states with different renormalons following this construction, a theoretically interesting enterprise but not likely to be of much direct phenomenological use.
The inclusion of the octet states of Eq. (91) is necessary for any consistent pNRQCD calculation of gluino pair production near threshold at NLO [54]. At NNLO the decuplet and 27-plet fields will also play a rôle. In fact such contributions, depending on the mass of the gluino (and on its existence), might be of bigger importance than in the case of tt production because there are more of them. This is an exciting and very clean-cut situation since v 2 and r −1 are bigger by a factor ∼ C A /C f than for quarkonia, such that all "soft" physics is clearly and extremely safely within the perturbative domain.
Let us finally mention that Λ B , the binding energy of the lightest glueballino, determines the size of the splitting between the adjoint singlet potential and the lowest adjoint hybrid potential at short distances, the latter of which, unfortunately, has never been determined in lattice simulations. This is very different from the case of fundamental sources where binding energies of heavy-light systems, Λ, are much smaller than the gaps between ground state and hybrid excitations. In "hadrinos", that contain stable adjoint sources, gluonic excitations would hence play a very prominent rôle and simple constituent-gluino models might fail terribly. Unfortunately, in nature we do not encounter such particles. It would however be most entertaining to confirm this expectation in lattice simulations.

Conclusions
We report compelling evidence that for distances around 1 GeV −1 the gluonic excitations of the static potential are in the short distance regime.
We are able to obtain a value for the lowest lying mass Λ B of the bilocal gluonic correlation functions with well controlled uncertainties, by fitting to the difference between the Π u and Σ + g potentials. The RS scheme result reads, Note that one should also add an extra error of order 10 % due to quenching to these numbers. With the information presented in this paper ν f can be run to different scales (see Fig. 14). We also obtain values for the masses of other gluelumps, listed in Table 4, as well as for the non-perturbative slope A Πu−Σ − u = 0.92 +0.53 −0.52 r −3 0 of the quadratic difference between the lowest two hybrid potentials.
In order to state sensible numbers for Λ B , the scheme for the renormalon cancellation has to be specified. Otherwise, very different numbers can be obtained, as we can see from a comparison of the result in the lattice and the RS schemes. One can translate from one scheme into the other in a renormalon-free way, order by order in perturbation theory and check whether both results are consistent with each other. We have been able to confirm this. If we use the gluelump results from Foster and Michael [12] at finite lattice spacing we obtain, Λ RS B (2.5 r −1 0 ) = 2.31 ± 0.04(latt.) ± 0.33(th.) +0.18 −0.19 (Λ MS ) r −1 0 , which is perfectly compatible with the result Eq. (103) above, albeit with slightly larger errors.
We also investigate the binding energy of heavy-light mesons in the static limit and to NNNLO in the conversion. We arrive at similar conclusions. For the binding energy we obtain the n f = 0 value, Λ RS (ν f = 2.5 r 0 ) = [1.17 ± 0.08(latt.) ± 0.13(th.) ± 0.09(Λ MS )] r −1 0 , which is in good agreement with the phenomenological value, obtained from the experimental Υ(1S) and B meson systems [15], Λ RS = [0.92 ± 0.22(th.) +0.15 −0.11 (Λ MS )] r −1 0 . We have demonstrated the internal consistency of our approach. Lattice predictions for Λ L B and Λ L at different lattice spacings have been studied. We have shown that the perturbative series, Eq. (80), relating Λ L B (a) and Λ L (a) with Λ L B (a ′ ) and Λ L (a ′ ), respectively, in the lattice scheme is free of renormalon singularities and has nice convergence properties, as indicated by the consistency with the non-perturbatively obtained values. In particular this means that from the knowledge of Λ L B and Λ L at a given lattice spacing values at different lattice spacings can accurately be predicted. We have studied the conversion of lattice predictions for Λ L B and Λ L into the RS scheme. This conversion is also dictated by a perturbative series which is free of renormalons. We have verified that the values in the lattice scheme indeed approach the results in the RS scheme with a convergent pattern and, remarkably, the ν f -scale dependence predicted by the RS scheme is reproduced, within errors. We remark that for the ν f -scale running it is possible to obtain a resummed nonperturbative expression in which the renormalon is cancelled and at the same time the log resummation is performed. We stress that the RS scheme used here is designed to smoothly converge to (MS-style dimensional regularization) perturbation theory at low orders in α s ; after all, the renormalon effect only sets in asymptotically at large orders in perturbation theory. Different values for Λ B and for Λ can be obtained in other schemes but only at the inconvenience of having large corrections to "standard" perturbation theory at low orders. In this sense we consider our approach "natural"; the RS scheme incorporates salient features of both, dimensional and lattice regularization. The approach readily benefits from results computed in the MS scheme, the scheme in which perturbative quantities are usually known to the highest order. On the other hand, by subtracting renormalons we encounter explicit power divergencies, which is exactly what one obtains with a hard lattice cut-off too.
Our model independent non-perturbative predictions can directly be incorporated into perturbative calculations, within effective field theories, or exploited in the context of QCD vacuum models or calculations based on non-local condensates. Obvious phenomenological applications in the context of EFTs are pNRQCD in the kinematic domain mv 2 < ∼ Λ QCD < mv, translating glueballino masses into RS or MS gluino masses within HGET (Heavy Gluino Effective Theory), or converting heavy-light meson masses into quark masses within HQET. We observe that Λ B ≈ (C A /C f )Λ ≈ m G /2, where m G denotes the mass of the lightest glueball. The first similarity is not necessarily surprising since there are technical parallels between Λ B , which corresponds to the binding energy of an adjoint source, and Λ, the energy of a fundamental source. We do not intend to advocate a constituent gluon picture. Nevertheless, it may seem reasonable that the binding energy of the glue to an adjoint source has about half the size of the energy of an entirely gluonic state. It should however be noted that the latter is an unambiguously defined state in the physical spectrum while for the binding energy Λ B we necessarily encounter the scheme and scale dependence that we discussed.
We have also investigated the scenario of gluinonia and other excitations in non-fundamental channels. While gluinos might not exist in nature and certainly do not form light bound states, such that phenomenological applications are limited, from a theoretical and conceptual point of view the existence of this part of the spectrum is very interesting. The inclusion of such potentials allows one to identify many combinations in which renormalons and other un-wanted contributions vanish, opening up a window to the study of non-perturbative short distance physics.