Higher-order condensate corrections to $\Upsilon$ masses, leptonic decay rates and sum rules

With the recent completion of NNNLO results, the perturbative description of the $\Upsilon$ system has reached a very high level of sophistication. We consider the non-perturbative corrections as an expansion in terms of local condensates, following the approach pioneered by Voloshin and Leutwyler. The leading order corrections up to dimension eight and the potential NLO corrections at dimension four are computed and given in analytical form. We then study the convergence of the expansion for the masses, the leptonic decay rates and the non-relativistic moments of the $\Upsilon$ system. We demonstrate that the condensate corrections to the $\Upsilon(1S)$ mass exhibit a region with good convergence, which allows us to extract $\overline{m}_b(\overline{m}_b) = 4214\pm37\,(\text{pert.})\,_{-22}^{+20}\,(\text{non-pert.})\text{ MeV}$, and show that non-perturbative contributions to the moments with $n\approx10$ are negligible.


Introduction
In recent years, the accuracy of the perturbative description of the bottomonium system has been extended to next-to-next-to-next-to-leading order (NNNLO). The full spectrum [1,2] 1 , the leptonic decay rate of the Υ(1S) [4] and the non-relativistic moments of the total e + e − → bbX cross section [5] have been determined and perturbation theory is well behaved. This has important phenomenological implications. For instance, some of the most precise determinations of the bottom-quark mass rely on the comparison of the perturbative expressions for the non-relativistic moments [5][6][7][8] or the masses of the Υ(1S) and η b (1S) resonances [9][10][11] and recently also the n = 2 states [12] with their experimental values.
Bottomonium can be treated as a non-relativistic system where the bottom-quark velocity v is of the order of the strong-coupling constant: v ∼ α s (m b v) 1. There is a large hierarchy between the dynamical scales m b (hard), m b v (soft) and m b v 2 (ultrasoft) of the system. The perturbative calculations [1][2][3][4][5] were performed using the effective theory potential non-relativistic QCD (PNRQCD) [13][14][15][16], where the hard and the soft scale have been integrated out and the only dynamical modes left are potential bottom quarks ψ and anti-bottom quarks χ, with energy and momentum of the order m b v 2 and m b v, respectively, as well as ultrasoft gluons and light quarks. The Lagrangian for perturbative calculations up to NNNLO takes the form where the coupling to the ultrasoft gluon field in the bottom-quark bilinear parts has been multipole expanded in the spatial components [17], the third line describes the interactions through spatially non-local potentials, which are given in [16,18], and the ultrasoft Lagrangian is a copy of the QCD Lagrangian which only contains the ultrasoft gluon and light quark fields.
Purely perturbative calculations within PNRQCD are valid when the ultrasoft scale m b v 2 is much larger than the QCD scale Λ QCD . This is certainly the case for top quarks, which are studied in [19,20] 2 , but is questionable in the bottomonium sector. Assuming the hierarchy holds, non-perturbative corrections can be incorporated in terms of local vacuum condensates as a power series in (Λ QCD /(m b v 2 )) 2 , following the approach of Voloshin and Leutwyler [23][24][25][26]. In this work, we compute higher-order corrections in this approach and assess the convergence of the series.
In the limit Λ QCD m b v 2 , the gluon field in the PNRQCD Lagrangian can be split into two parts The superscripts denote the ultrasoft and the non-perturbative gluon field with momentum of the order m b v 2 and Λ QCD , respectively. All couplings of the nonperturbative component to other modes must be multipole-expanded because the non-perturbative field with a large wavelength of the order 1/Λ QCD cannot resolve the dynamics of the potential bottom quarks or of the ultrasoft gluons. A convenient gauge choice for the non-perturbative gluon field is given by Fock-Schwinger gauge which removes the coupling of the bottom quarks to the A np 0 field. The leading non-perturbative contribution in the PNRQCD Lagrangian then takes the form of a chromoelectric dipole term L non-perturbative = ψ † (−g s x · E np (0, 0) + . . . ) ψ + χ † (−g s x · E np (0, 0) + . . . ) χ, (4) and is of the order m 2 b v 2 Λ 2 QCD because E np ∼ Λ 2 QCD and the strong coupling at the QCD scale is counted as order one. This implies that the chromoelectric dipole coupling to the non-perturbative gluon field is suppressed by v(Λ QCD /(m b v 2 )) 2 with respect to the leading order Lagrangian. The time-dependent terms in the multipole expansion and the expanded couplings between the non-perturbative and ultrasoft modes in L ultrasoft are not required for the leading order condensate corrections. Their relevance at higher orders is assessed in Section 3, where we discuss the NLO QCD corrections to the leading term in the Voloshin-Leutwyler approach.
The condensate corrections to the considered observables can be extracted from the non-relativistic Green function at the origin where E = √ s−2m b is the non-relativistic energy of the system and the Hamiltonian has the formĤ =Ĥ bb +Ĥ np +Ĥ D + . . . .
The bottomonium Hamiltonian follows from (1) and is given by a perturbative series in g s ∼ v 1/2 : with the color-singlet and color-octet projectors where the color indices are assigned in the same way as in the potential term in (1). The LO HamiltonianĤ bb,0 is of the order m b v 2 . The non-perturbative dynamics at the scale Λ QCD are described by the HamiltonianĤ np which is of the order Λ QCD . The leading interaction between the bottomonium and non-perturbative sector is given by the chromoelectric dipole term with ξ A abcd = T A ab δ cd + δ ab T A cd when the color indices are again assigned in the same way as in the potential term in (1), which is of the order Λ 2 QCD /(m b v). Assuming Λ QCD m b v 2 the interactionĤ D and the non-perturbative HamiltonianĤ np can therefore both be treated as perturbations and the physical state factorizes into the product of a bottom-antibottom state |0 bb at zero spatial separation and the non-perturbative vacuum state |0 np . The expansion of the Green function (5) in powers of Λ QCD then takes the form is the perturbative part of the Green function and we adopted the notation of [27]: The propertiesĤ np |0 np = 0 and 0|g s (E np ) A i |0 np = 0 have been used to remove insertions ofĤ np that are not in between insertions ofĤ D and single insertions of H D . Terms with an odd number ofĤ np insertions between the twoĤ D insertions vanish, because they can be related to the vacuum expectation values of operators with odd numbers of Lorentz indices by using Lorentz invariance, see [23].
The first term in (11) is the purely perturbative part. The sum contains the leading non-perturbative contributions, which are proportional to vacuum expectation values of operators of even dimensions and are suppressed by v 2 (Λ QCD /(m b v 2 )) 4,6,8,... with respect to the perturbative expression. The contributions of dimension four and six are shown in Figure 1. The extra suppression factor v 2 is present because terms without at least two insertions ofĤ D vanish. The dimension-four correction contains the gluon condensate αs π G 2 and has been studied in [4,[23][24][25][26]. The dimension-six correction to the masses and leptonic decay rates of the Υ(N S) resonances has been calculated in [27].
Due to the extra suppression factor v 2 , smallness of the dimension-four contribution is not sufficient to demonstrate the convergence of the expansion in Λ QCD /(m b v 2 ) and the calculation of higher-order condensate corrections is necessary to gain more insight. We compute the leading corrections up to dimension eight in Section 2. The NLO potential corrections to the dimension-four condensate contribution are determined in Section 3. The size of the condensate corrections to observables in the Υ system is discussed in Section 4. We conclude in Section 5.
2 Leading order condensate corrections of dimensions four, six and eight The leading order condensate corrections are finite and can be computed in four dimensions in position space. Inserting spatial integrations, the dimension-four contribution in (11) takes the form The integrals can be evaluated using the known representations of the LO Green function G (1,8) 0 , where the superscript indicates whether the bottomonium state is in a color singlet (1) or octet (8) configuration. It is convenient to decompose the Green function in terms of partial waves where l is the quantum number of the angular momentum of the bottom pair and P l (z) are the Legendre polynomials. We use an integral representation from [28], valid for r < r and a sum representation from [24,29], with the Laguerre polynomials We have defined the variables In the following we take N c = 3 and use the variable λ ≡ λ (1) = −8λ (8) . We use the integral representation (16) for the color-singlet Green functions and the sum representation (17) for the color-octet Green function. The angular integrals in (14) project out the S-wave component of the color-singlet Green functions and the Pwave component of the color-octet Green function. We obtain where The coefficients H c contain the remaining integrations and read The sum in (20) yields where ψ and ψ 1 are the polygamma functions of order 0 and 1, respectively. The condensate corrections to the S-wave energy levels E N and the wave functions at the origin |ψ N (0)| 2 can be obtained from the expansion of (24) for λ near positive integer values N as described e.g. in [20,30,31]. The results are given in Appendix A. The same strategy can be applied for the calculation of the dimension six and eight condensate corrections. Again, the angular integrals project out the S-wave component of the color-singlet Green functions and the P-wave components of the color-octet Green functions. We find and where and the coefficients K c read Since K c (a, b) is only non-vanishing for |a − b| ≤ 1 the multiple sums in (24) and (25) are reduced to a single sum, which can be solved in terms of polygamma functions.
The lengthy results are available as ancillary files with the arXiv version of this article. The dimension six and eight contributions to the energy levels and wave functions are given in Appendix A.

Dimension four contribution at NLO: Potential contributions
The NLO corrections to the dimension-four condensate contribution involve an insertion of the NLO Coulomb potential as shown in Figure 2 and ultrasoft loops as shown in Figure 3. The upper panel of Figure 3 shows the diagrams with ultrasoft gluon loops, where the gluon coupling to the color-octet state originates from the leading term g s A us 0 (t, 0) in the multipole expansion. The equivalent coupling to the color singlet state vanishes because the ultrasoft gluons cannot resolve the spatial separation of the bottom-antibottom state and the net color charge vanishes in the singlet state. 3 The diagram in the lower panel of Figure 3 shows the contribution from the light-quark condensate qq with q = u, d, s which is also counted as dimension four because, due to chirality suppression, the quark condensate qq only appears together with one power of the light quark mass m q which is of the order Λ QCD . There is a number of other effects that could possibly contribute at that order: • An α s correction to the Wilson coefficient of the chromoelectric dipole operator (4). The Wilson coefficient was found to be trivial up to O(α 2 s ) in [33].
terms in the multipole expansion (4) of the gluon coupling to bottom quarks in the spatial components. They are identical to the multipole expansion of the coupling to the ultrasoft gluon field and were determined in [34], where they are denoted as h SO . There is no NLO contribution from these terms because they either have vanishing tree level Wilson coefficients and are thus suppressed by an additional power of α s ∼ v or involve the chromomagnetic instead of the chromoelectric field, which only yields a vanishing condensate 0 E A i B B j 0 = 0 at NLO. 3 See also [32] for a more formal argument based on a field transformation. • Contrary to the ultrasoft gluon-bottom coupling, the interactions of the nonperturbative gluon field must also be multipole expanded in the time component. The expansion of the A 0 component is trivial due to our gauge choice (3) and already the linear term t(∂ 0 A i )(0, 0) in the expansion of the spatial component is only relevant at higher powers. Thus, no contributions of this type need to be considered at NLO.
The potential corrections are determined below, whereas the ultrasoft contribution is postponed to future work. The NLO correction to the Coulomb potential is given by where the color factors are given by where n f is the number of massless quarks. Denoting the contribution from the left (right) diagram in Figure 2 by DVD (DDV), we find The first triple insertion function takes the form where The full u-dependence of (32) is not needed here. To evaluate (30) we only need the value and the first derivative at u = 0. We obtain The derivative of (32) at zero can be solved by applying the methods used for the Coulomb triple insertion in [30]. This yields where L λ = ln(λµ/(m b α s C F )) and The second triple-insertion function yields where and H(u, k) is defined as in [30]. Also here, we only need the value and the first derivative at u = 0: The infinite sums in (31) and (36) converge quickly and can be truncated with negligible uncertainty at s i ∼ 30 for the numerical evaluation of the Green function. The contributions to the energy levels and wave functions from the potential corrections can be extracted by expanding (31) and (36)

Phenomenology of condensate corrections
The size of non-perturbative corrections to the moments and to the properties of the Υ resonances has been strongly disputed for various reasons. First, the assumption Λ QCD m b v 2 is questionable and is certainly only valid for a limited number of observables in the Υ system. Here, we perform an unbiased analysis of the expansion in terms of local condensates and assess the validity based on its convergence. The breakdown of this expansion is a clear indication that the above assumption is inappropriate.
Furthermore, the numerical values of the local condensates are very uncertain. The condensate O 0 is proportional to the gluon condensate and we will use the standard value αs π G 2 SVZ = 0.012 GeV 4 from [35] below, unless indicated otherwise. We note however, that significantly larger values have also been obtained in the literature, see e.g. [36,37]. Clearly, the situation is even more uncertain for the higher-dimensional condensates. Since our main objective is the assessment of the convergence properties, we rely on naive rescaling The value of O 0 is scale independent since the gluon condensate αs π G 2 is not renormalized. We neglect the scale dependence of the higher-dimensional condensates which is very weak compared to that of the coefficients which contain large powers of α s . The estimate O naive 1 is in good agreement with the result of [27], where an expression for O 1 in terms of the dimension-six gluon condensate G 3 and the quark-condensate qq has been derived based on the factorization hypothesis. The analysis of [27] also shows that O 1 is only weakly scale dependent.
In addition, the corrections to the masses and leptonic decay rate depend strongly on the renormalization scale, because large powers of α s appear in the ratios (21) and (26). The fact that different powers of α s appear in the contributions of different dimensions also complicates the assessment of the convergence and different conclusions have been drawn based on different scale choices. We distinguish the scale µ c , used in the condensate corrections, from the renormalization scale µ in the perturbative contribution. The main motivation for the calculation of the potential corrections to the dimension-four contribution has been to gain more insight into the appropriate scale choice for µ c by considering the convergence of the perturbative series. We note that the potential corrections contain all logarithms ln(µ c ) that are required to cancel the µ c dependence of the dimension-four contribution at NLO. The ultrasoft correction must therefore be free of logarithms ln(µ c ) and is less scale dependent, which justifies performing this analysis based on incomplete NLO corrections. Scales below 0.8 GeV are not considered below, because the value of α s and perturbation theory in general become unreliable in this regime.

The Υ(1S) mass
First, we briefly review the status of the purely perturbative prediction for the mass of the Υ(1S) resonance. We use QQbar Threshold [20,38] in the PS mass scheme [39] with the input value m PS b = 4.532 +0.013 −0.039 GeV from [5,6]. The effects of a nonzero charm-quark mass are included up to NNLO [5] using the mass m c (3 GeV) = 993 MeV from [40,41]. The default values of QQbar Threshold are taken for the strong coupling α s (m Z ) = 0.1184 ± 0.0010 and all other parameters, and QED corrections are taken into account with NNLO accuracy. The result is shown in the top panel of Figure 4. We observe that the convergence is best for scales that are considerably larger than the soft scale µ ∼ m b α s (µ)C F . This has motivated the authors of [11] to choose a central scale of 5.35 GeV, which is significantly larger than that of [9,10,12] (2.5 and 1.9 GeV, respectively) and leads to a much smaller estimate for the perturbative uncertainty. Here, we choose a central scale of 3 GeV such that a variation by factors 1/2 and 2 covers the choices of [9][10][11][12]. The perturbative expansion takes the form In addition to the perturbative uncertainty from scale variation, we also take into account the parametric uncertainty from the bottom-quark PS mass and we use the size of the charm-quark mass effects up to NNLO as an estimate for the missing NNNLO correction. The parametric uncertainty from the strong coupling is small in the PS mass scheme and is neglected. The condensate corrections with the values of (40) are shown in the lower panel of Figure 4. At the considered orders, the mass scheme is ambiguous and we use the one-loop pole mass in the condensate contribution. In the PS scheme, the condensate contributions are slightly enhanced and the convergence is slightly worsened but the overall conclusions are unchanged. We observe that the potential corrections to the dimension-four condensate contribution stabilize the behaviour under scale variation and show a clear preference for rather small scales around µ c = 1.2 GeV, which we take as the central value. The condensate contribution takes the form The grey band in Pole scheme  [9,10], KMS'15 [11], MO'17 [12]) and non-relativistic sum rules (HRS'12 [8], BMPR'14 [5,6]). The bottom-quark masses obtained from lower orders in pure perturbation theory, while retaining all known condensate contributions, are shown as well. The order of the PS-MS mass relation has been correlated with the order in perturbation theory.
ultrasoft NLO correction. Combining the perturbative and condensate contributions we find where we have symmetrized the uncertainty from variation of the renormalization scale µ, by taking the maximum of the positive and negative error. The perturbative uncertainty is obtained by adding the errors from variation of µ and α s as well as our estimate of higher-order charm-quark mass effects in quadrature. The variation of the scale µ c and of the values of the condensates is combined into the non-perturbative uncertainty. The result (44) is converted to the MS scheme at NNNLO [42,43] using QQbar Threshold. We distinguish the scale µ m used in the conversion, which is set to m PS b , and estimate the uncertainty through variation of µ m by factors of 1/2 and 2 and symmetrization as described above. We find The result shows good convergence and agrees with other recent determinations of m b from the data on the Υ system as shown in Figure 5. In conclusion, our analysis demonstrates that the determination of the bottom-quark mass from the Υ(1S) mass is possible with a total uncertainty of the order of ±45 MeV. It should however be noted that this approach to the determination of the bottom-quark mass is on a less sound footing theoretically than the extraction based on non-relativistic moments with n ≈ 10, which are discussed in Section 4.4.

The Υ(2S) mass
We repeat the above discussion for the Υ(2S) mass. The scale dependence of the perturbative result is shown in Figure 6. Since the soft scale is lower for the n = 2 states, we reduce the central scale to 2 GeV, where the perturbative series takes the form M pert Υ(2S) (2 GeV) = (9 534 + 198 + 154 + 116) MeV.
As the plot shows, the convergence is rather slow, independently of the choice of scale. We also note that the charm-mass effects at NNLO are +39 MeV and significantly larger than for the Υ(1S) mass (+8 MeV). As we argued in [5], the charm-mass effects are a measure for the IR sensitivity of an observable. Thus, the significantly larger value is an indication that the non-perturbative correction should be considerably larger and less convergent for the Υ(2S) mass than for the the Υ(1S) mass.
Turning to the condensate corrections, which are shown in the lower panel of Figure 6, we can confirm this expectation. The expansion already breaks down for µ c = 0.8 GeV, where the individual contributions are MeV. (47) At lower scales, the use of perturbation theory cannot be justified. Thus, while we cannot rule out the convergence of the local condensate expansion unambiguously Pole scheme due to the large uncertainties of the O i , clearly no reliable prediction for the nonperturbative contribution can be obtained like this. A more promising approach to the Υ(2S) mass is to assume the hierarchy Λ QCD ∼ m b v 2 m b v. Then, the ultrasoft contribution takes the form of a non-local condensate instead of a perturbative correction [15,44,45]. This implies that the leading non-perturbative correction is of the order ∆M non-perturbative which is formerly of NNLO, and the conclusion that the local condensate expansion breaks down is equivalent to the statement that the Υ(2S) system is outside the radius of convergence for the presently unknown function ρ. In this scenario, the perturbative NNNLO results, which contain the perturbative evaluation of the ultrasoft contribution cannot be used and we have to resort to the NNLO expressions. The result for the Υ(2S) mass reads where the estimate for the non-perturbative contributions follows from the assumption that the function ρ in (48)

The Υ(1S) → l + l − decay width
The perturbative NNNLO result for the leptonic decay width of the Υ(1S) resonance has been obtained in [4]. Here, we repeat their analysis including charm-mass effects up to NNLO, which increase the leptonic width by 0.03 keV. The scale dependence is shown in Figure 7 and we adopt 3.5 GeV as the central scale. The perturbative series stabilizes at NNNLO but falls short of the experimental value Γ exp (Υ(1S) → l + l − ) = 1.340 ± 0.018 keV by about 20%. Following [4] we determine the scale uncertainty from variation between 3 and 10 GeV. The other input parameters are varied as above.  The condensate contributions are shown in the lower panel of Figure 7. Using the same central scale µ c = 1.2 GeV as for the Υ(1S) mass, we obtain (51) Focusing first on the leading-order contributions, we see that the expansion converges and yields a contribution of 0.27 keV that closes the difference between the perturbative and the experimental value. Compared to the Υ(1S) mass the expansion breaks down at a smaller scale around 1.6 GeV.
However, with the addition of the potential corrections to the dimension four contribution, the agreement is destroyed. The potential correction already exceeds the LO term at the scale 0.7 GeV and becomes twice as large at 0.9 GeV. This apparent breakdown of the perturbative series makes it impossible to give a reliable estimate of the non-perturbative contribution. However, it is conceivable that the large potential corrections are compensated by the missing ultrasoft correction, thus stabilizing the perturbative expansion of the dimension-four contribution. Therefore, no definite conclusions about the validity of the local condensate expansion for Γ(Υ(1S) → l + l − ) can be drawn without a calculation of the full NLO corrections to the dimension-four contribution.

The non-relativistic moments
The moments M n of the normalized inclusive bb production cross section in e + e − collisions with the center-of-mass energy s, are defined as The normalized cross section is related to the bottom-quark contribution Π b to the photon vacuum polarization by the optical theorem R b (s) = 12π Im Π b (s + i ). The contour C must be closed around s = 0 without crossing the branch cut for real s ≥ M 2 Υ(1S) . The perturbative contributions to the moments up to NNNLO have been discussed in detail in [5]. The non-perturbative corrections can be determined by inserting the condensate contribution to the cross section where E = √ s − 2m b and c (1) v = −8C F is the hard matching coefficient of the vector current, into (53). Following the discussion in [5] we choose not to expand the prefactor 1/s around 1/(4m 2 b ). Contrary to the perturbative contribution, we cannot split the condensate corrections into a resonance and continuum part, since both are separately divergent [5]. The total corrections to the moments are however well-defined and can be computed numerically using the representation of the moments (53) involving contour integration or, in principle, analytically by taking derivatives at q 2 = 0.
The scale dependence of the dimension-four contribution are shown in Figure 8. Results are given in the pole mass scheme using the same inputs as given above. We refrain from using the PS or other threshold mass schemes, because the perturbative expansion in these schemes becomes unstable in large regions of the scale µ c . This can be traced back to the appearance of large powers of λ in the expression for the Green function (24), which are expanded in the PS mass scheme as where is the N i LO contribution to the PSpole mass relation. 4 This is reminiscent of the destabilization of the NLO correction to the gluon condensate contribution [46] to the relativistic moments in the MS scheme [47].
The top panel of Figure 8 shows the leading order result. Although the contribution is proportional to R 4 ∝ α −6 s (µ c ), its absolute value decreases for larger scales µ c . Given that the condensate corrections to the Υ(1S) mass and leptonic decay rate become unstable for scales larger than about 2 GeV and 1.6 GeV, respectively, this behaviour must be caused by very pronounced cancellations between the contribution from the Υ(1S) resonance and the remaining resonances and the continuum (rest), which was pointed out in [5]. For the tenth, sixteenth and twenty-fourth moment, this cancellation is effective at the level of one part in 139, 52 and 20 at the scale µ c = 2 GeV and at one part in 1530, 659 and 297 for µ c = 10 GeV and the growth of the degree of the cancellation for higher scales dominates over the growth of the factor α −6 s (µ c ). While this qualitative behaviour is expected due to the reduced infrared  sensitivity of the moments compared to the properties of the Υ(1S) resonance, the extent of the cancellations and the resulting smallness of the corrections is rather surprising, especially for larger values of n 16 where power counting predicts a breakdown of the expansion in powers of (nΛ QCD /m b ).
The results including the potential NLO corrections are shown in the lower panel of Figure 8. Above 3 GeV the corrections do not exceed the size of about ±20% for the considered moments. This is due to even more pronounced cancellations within the potential corrections which are effective up to about one part in 10 4 (n = 10), 5 · 10 4 (n = 16) and 2 · 10 3 (n = 24) at µ c = m b . The corrections mainly have the effect of stabilizing the scale dependence at lower scales µ c 2 GeV, such that we find good behaviour of the dimension-four contribution at partial NLO over the considered range of scales between 1 and 10 GeV.
We can try to assess the convergence of the condensate expansion based only on the dimension-four results. Compared to the perturbative result, they are of the relative order 1/n × (nΛ QCD /m b ) 4 , where the extra factor of 1/n accounts for the v 2 -suppression from the two insertions of the dipole operator. Thus, we expect a breakdown of the condensate expansion in nΛ QCD /m b when the dimension-four contribution is of the relative size 1/n. From the lower plot in Figure 8 we deduce that this point is reached in the ballpark of n ≈ 20, where the condensate contribution is of the size of -4% of the experimental moment at its peak, which is compatible with the expectation from the power counting argument.
In Figure 9, the relative contributions of dimension six (upper panel) and eight (lower panel) are shown. Both are significantly smaller than our expectation based on the putative breakdown of the expansion around n ≈ 20, which would imply that the dimension six and eight corrections are both of the order 1/n ≈ 0.05. This smallness is the result of cancellations between the contribution from the Υ(1S) resonance and the rest that are even stronger than at dimension four. Explicitly, they are at the level of about one part in 3 · 10 5 (n = 10), 5 · 10 4 (n = 16) and 10 4 (n = 24) at dimension six and about one part in 10 8 (n = 10), 10 7 (n = 16) and 2 · 10 6 (n = 24) at dimension eight. We believe that the reason for this behaviour is the off-shellness of the moments which are defined as derivatives of the vacuum polarization function at q 2 = 0, far away from the physical cut at s ≥ M 2 Υ(1S) . This off-shellness effectively acts as an IR cutoff and suppresses higher-dimensional corrections, which probe the IR regime. On the other hand, the properties of the Upsilon resonances, that we discussed above, are on-shell quantities and the higher-dimensional condensate contributions do not appear suppressed with respect to our expectations from power counting.
From the point of view of the convergence of the condensate expansion, it appears  that the moments can be described reliably up to values of n much larger than 20. However, as pointed out in [5], the validity of quark-hadron duality must be questioned when the moment is completely saturated by lowest state. This is the case for the higher values considered here, where the relative contribution of the Υ(1S) to the experimental moments amounts to 95% for n = 20 and 97% for n = 24 [5]. By the term 'violation of quark-hadron duality' we refer to contributions which have a trivial Taylor expansion and are, therefore, not captured by the condensate expansion, like e.g. exponential terms of the form exp(−m b /(nΛ QCD )). Behaviour that is consistent with the presence of such contributions has been observed in the 't Hooft model [48]. 5 However, the size of these contributions in four-dimensional QCD is difficult to quantify and we do not attempt this here. We note, however, that the exponential terms originate from coherent soft fluctuations [49], e.g. from contributions where the off-shellness is distributed among many soft lines carrying momenta of the order Λ QCD , which pushes the bottom pair close to its mass shell. It is conceivable that such an effect does not experience a similar suppression from the effective IR cutoff as the higher-dimensional condensate contributions.
In the range nΛ QCD ∼ m b the exponential exp(−m b /(nΛ QCD )) is of order one and we cannot exclude that duality violation effects are relevant at the high accuracy we require for reliable determinations of the bottom-quark mass. We conclude that, in practice, the range of moments is limited by our knowledge of the validity of quarkhadron duality and not by the convergence of the condensate expansion and advise that moments with n 16 are not used for determinations of the bottom-quark mass. On the other hand, for n ≈ 10 duality-violating effects are exponentially suppressed and the condensate expansion provides a reliable determination of the non-perturbative effects. Our results given in Figure 8 and 9 show that the condensate contributions in this region are in the subpercent range and can safely be neglected compared to the perturbative uncertainties.

Conclusions
We have determined the leading order condensate corrections to the Υ(N S) masses, leptonic decay rates and sum rules up to and including dimension eight. In addition the potential NLO corrections to the dimension-four contribution have been computed, which allows us to assess the preferred scale choice in the condensate corrections. Our results suggest that the expansion is well behaved for the mass of the Υ(1S), but breaks down for the higher states. The former observation has been used to determine the bottom-quark mass with the results given in (44) and (45).
The leading order condensate corrections to Γ(Υ(1S) → l + l − ) have a small window of convergence for 0.8 GeV µ c 1.6 GeV, where they lead to good agreement with the experimental value, but the partial NLO corrections to the dimension-four contribution exceed the leading order correction and cause us to question the perturbative stability. Thus, a final verdict for the leptonic decay rate of the Υ(1S) is only possible once the missing ultrasoft correction has been calculated.
Last but not least, we have considered the non-relativistic moments (53). We find extremely good convergence of the higher-dimensional condensate contributions which clearly shows that non-perturbative contributions to moments with n ≈ 10 are negligible. On the other hand, we cannot unambiguously exclude the possibility of relevant violations of of quark-hadron duality for n 16 despite the surprising smallness of the dimension six and eight corrections. Thus, the non-relativistic moments with n ≈ 10 remain the theoretically cleanest approach for determinations of the bottom-quark mass from the Υ system.

A Condensate corrections to the energy levels and wave functions
We give the results for the condensate corrections to the energy levels and the wave functions at the origin of the S-wave bottomonium states. The contributions are parametrized as where the leading order expressions are given by the perturbative corrections of relative order α i s are e in agreement with the results from [24,25]. For the dimension-six corrections to the energy levels and wave functions, we find The correction to the energy levels is identical to the result of [27]. Our result for the wave function correction however differs from the one in [27] in the coefficients in the square bracket that multiply powers of N , while the constant term is in agreement. Numerically the difference is tiny, dropping from 3 permille for N = 1 to 1.8 permille for N = 10. Our dimension-eight results read The NLO corrections to the dimension-four condensate contributions have the form where the potential terms are where L N = ln(N µ/(mα s C F )), S 1 (x) = x k=1 k −1 is the analytic continuation of the harmonic number to non-integer values and S i = N k=1 k −i without an explicit argument is the N th harmonic number of rank i. The ultrasoft corrections e us N and f us N are currently unknown.