Perturbative heavy quarkonium spectrum at next-to-next-to-next-to-leading order

We compute the energy levels of some of the lower-lying heavy quarkonium states perturbatively up to O(alpha_s^5*m) and O(alpha_s^5*m*log[alpha_s]). Stability of the predictions depends crucially on the unknown 4-loop pole-MSbar mass relation. We discuss the current status of the predictions with respect to the observed bottomonium spectrum.

During the past decade, spectroscopy of heavy quarkonium states (in particular the bottomonium states) has provided an important testing ground of perturbative QCD. On the one hand, we have at our disposal relatively many terms of the perturbative expansions. On the other hand, the system in question is a small-size (compared to the typical QCD scale Λ −1 QCD ) color-singlet object, from which large part of infra-red (IR) degrees of freedom decouple. In fact, the discovery of a cancellation of the IR renormalons in the energy levels of heavy quarkonium states led to a drastic improvement in the predictability of the energy levels within perturbative QCD [1]. We observe that stability and convergence properties of the perturbative series for the energy levels are fairly good, even in comparison to other observables of a heavy quarkonium, such as production cross sections, transition rates, or partial decay widths. Important applications of the spectroscopy are precise determinations of the heavy quark masses from the lowest-lying energy levels. The bottom and charm quark masses have been determined, and in the future the top quark mass is expected to be determined accurately in this way. (See [2] for reviews.) The full O(α 4 s m) corrections to the energy levels were computed in [3]. Analyses of the bottomonium spectrum, which incorporate the renormalon cancellation and the perturbative corrections up to this order, have shown that the gross structure of the bottomonium spectrum, including the levels of the n = 1, 2 and some of the n = 3 states (n is the principal quantum number), is reproduced reasonably well within the estimated perturbative uncertainties [4].
During the subsequent years, our understanding on the energy of a heavy quarkonium system based on perturbative QCD has been advanced. Stability and agreement with experimental data of the predictions for the energy levels are predominantly determined by the prediction for the static QCD potential V QCD (r). After canceling the renormalon in V QCD (r) (in various schemes), perturbative predictability improves and the predictions for V QCD (r) agree with lattice computations or typical phenomenological potentials in the relevant distance range [5,6]. Furthermore, by increasing the order of the perturbative expansion, the range of r, where convergence and agreement are seen, extends to larger r [13]. The details, however, depend on the schemes adopted for canceling the renormalon.
Taking these features into account, some improvements of the predictions for finer structures of the bottomonium spectrum have been examined. Including all the known terms of V QCD (r) in the zeroth-order Hamiltonian, the fine and hyperfine splittings as well as the splittings between S-and P -wave levels have been computed in a specific organization of perturbative series [7]. This prescription enables to incorporate the effects of the rise of V QCD (r) at larger r on the wave functions, and this results in a better agreement of the above splittings with the experimental data. (See also [8].) [26] In the meantime, computations of the O(α 5 s m) and O(α 5 s m log α s ) corrections to the energy levels have made progress. Development of effective field theories, such as potential non-relativistic QCD (pNRQCD) [9] or velocity non-relativistic QCD [10], enabled systematic computations of the higher-order corrections, by separating the different kinematical regions involved in the corrections.
Within pNRQCD, the corrections consist of two parts, the contributions from the potential region and ultrasoft (US) region. The next-to-next-to-next-to-leading order (NNNLO) Hamiltonian, which dictates the contributions from the potential region, was computed in [11] (besides the 3-loop corrections to V QCD (r), a 3 , which were computed later in [12,13]). It is a straightforward (but cumbersome) computation to obtain the energy levels of the Hamiltonian in perturbative expansions analytically. The contributions from the US region contain, besides the part calculable analytically, a QCD analogue of the Bethe logarithm for the Lamb shift in QED. The QCD Bethe logarithm for each state can be written as a oneparameter integral of elementary functions [14]. Up to now, the O(α 5 s m log α s ) correction for a general state labeled by the quantum numbers (n, l, s, j) was computed in [15], while the O(α 5 s m) and O(α 5 s m log α s ) corrections for a general S-wave state (n, j) were computed in [16]. We have confirmed these results. (See also [17]   Since the analytic expressions plus integral forms are too lengthy to be presented here, and since one-parameter integrals need to be evaluated numerically for individual (n, l)'s in any case, we present the results numerically for some (n, l)'s. (The full formula and the derivation will be presented elsewhere.) In particular, we present the corrections necessary for all the observed bottomonium states whose masses are listed in [18] and which lie below the threshold for decays into two B mesons (2M B = 10.558 GeV).
We consider a bound-state composed of a quark (with the pole mass m pole ) and its anti-quark. The energy of the bound-state X, identified by (n, l, s, j), is given by Here, C F = 4/3 denotes the color factor; α s (µ) denotes the strong coupling constant in the theory with n l active flavors only, renormalized at the renormalization scale µ, and defined in the modified-minimal-subtraction (MS) scheme; P i (L) denotes an i-th-degree polynomial of L. α s (µ) obeys the renormalization-group (RG) equation where β i represents the (i + 1)-loop coefficient of the beta function. The only part of P i (L) that is not determined by the RG equation for E X is c i ≡ P i (0). [27] For i = 3, we have Our results of c 3 are listed in Tab. I, given as functions of (s, j), n l and log[α s (µ)] for fixed (n, l)'s. [28] Here, We neglect the masses of n l light quarks. The nonlogarithmic terms of the P -and D-wave levels are new.
Using the NNNLO results we compute the energies of the observed bottomonium states and compare them with the experimental data. We follow the prescription used in the analyses [4]. [29] To cancel the renormalons, we express the pole mass in terms of the quark mass defined in the MS scheme (MS mass) as where m ≡ m MS (m MS ) denotes the MS mass renormalized at the MS mass scale. The 4-loop constant d 3 , which is needed for our analysis, is not known yet. Up to now there exist some estimates of d 3 [19][20][21][22]. We adopt the estimate [22], obtained from perturbative stability of the energy of a static quark pair in the following manner. The upper bound of the estimate is determined by requiring stability of the perturbative prediction for E tot (r) ≡ 2m pole + V QCD (r) at NNNLO at least up to the same r as NNLO. In particular, as the value of d 3 exceeds its estimated upper bound, the perturbative prediction for E tot (r) becomes unstable rapidly. The lower bound of the estimate is obtained by requiring that the difference between the NNLO and NNNLO predictions for E tot (r) be within an O(Λ 3 QCD r 2 ) uncertainty. When d 3 is chosen within the estimated range, qualitatively the prediction for E tot (r) becomes stable and the series exhibits a reasonably convergent behavior.
After rewriting E X in terms of m and α s (µ) via eq. (8) and the solution to eq. (3), we expand E X in α s (µ). To  . E (i) X denotes the i-th order term of the ε-expansion. E pert Numerical values except in the second and last columns are in GeV. The last row represents the spin-averaged 3Pj energy for j = 0, 1, 2 with the weight factor 2j + 1. make the cancellation of the renormalons explicit, we reorder the series in the so-called "ε-expansion scheme" [23]. See [4] for details. We set α s (M Z ) = 0.1184, n l = 4, and m is fixed such that the energy of the Υ(1 3 S 1 ) state coincides with the experimental value. We vary the value of d 3 within (0.95 +0.01 −0.05 )×d large-β0 3 , which is the stability range of E tot (r) [30], where d large-β0 3 ≈ 1324.49 for n l = 4 [19].
[31] We find that practically the stability of the perturbative prediction for E X is determined by the stability of the perturbative prediction for E tot (r). In fact, the scale dependence and convergence properties of E X are qualitatively similar to those of E tot (r). In Fig. 1 we show the scale dependence of E X for the χ b (2 3 P 0 ) state, which is one of the highest states predicted reliably at NNLO. The NNNLO prediction is stable if d 3 is within the above range. If , an optimal value with respect to the stability of E tot (r). Then, for each state X, we fix the scale µ = µ X by demanding stability of E X against variation of the scale: This scale exists for all the bottomonium states considered here. The convergence behaviors of the perturbative expansions are reasonable. This means that the predictability range extends to higher levels compared to the NNLO predictions. We list the perturbative predictions and the experimental data in Tab tom quark MS mass, fixed on the Υ(1 3 S 1 ) state, is given by [33] m In Fig. 2 we compare the experimental data and the predicted bottomonium spectrum. We see a reasonable agreement for the gross structure of the spectrum. (See, however, the comments on uncertainties below.) If d 3 is within the range (0.90-0.96) × d large-β0 3 , convergence properties and stability of the predictions are qualitatively similar to those listed in Tab. II, although the level of agreement with the experimental data varies. As we raise (reduce) the value of d 3 , level spacings among different states increase (decrease). Variations are larger for higher states. If d 3 is raised above 0.96 d large-β0 3 , the extremum scale µ X disappears quickly from higher levels. If d 3 > 1.2 d large-β0 3 , µ X no longer exists even for the Υ(1 3 S 1 ) state.
In principle, we can estimate uncertainties of the predictions for E X originating from various sources, similarly to the analyses [4], for each given value of d 3 . With respect to the uncertainties, we can discuss agreement or disagreement with the experimental data. It would eventually lead to, for instance, quantification of nonperturbative contributions to individual energy levels. Since, however, stability of the predictions for E X depends crucially on d 3 , we consider the current status to be too premature to do this quantitatively. Here, we briefly comment on uncertainties. (i) Dependence on α s (M Z ): If we vary α s (M Z ) within the current uncertainties ±0.0007 [18], variation of the energy levels [after fixing m on Υ(1 3 S 1 )] is fairly small and minor as compared to uncertainties from other sources. (ii) Non-zero charm mass effects: Although a full account of non-zero charm mass effects in loops requires a separate analysis of its own, the analysis at NNLO indicates that the level spacings among higher levels increase due to the decoupling of the charm quark. Phenomenologically this indicates that a smaller value of d 3 may be favored for a better agreement with the experimental data after inclusion of these effects. (iii) Higher-order effects on level splittings: As already explained, inclusion of higher-order effects of V QCD (r) in the bound-state wave functions increases the fine and hyperfine splittings as well as the S-P splittings, which improves the agreement with the experimental data. We note that concerning the former splittings all the NNNLO corrections included in the present analysis have already been included in [7], so that the differences from the our results stem only from higher-order effects. Let us comment on non-perturbative contributions to the bottomonium energy levels. In general there are two ways to compute a physical quantity whose major contributions come from UV region. One way is to compute thoroughly within perturbative QCD. The other way is to compute by factorizing UV and IR contributions within a Wilsonian low energy effective theory. In the former computation, there are well-established prescriptions to estimate uncertainties of the prediction within perturbative QCD. Empirically estimates of perturbative uncertainties are approximated well by IR renormalons, in the case that IR renormalons turn out to be large. In the latter computation, UV contributions are encoded in the Wilson coefficients, which are free from IR renormalons and have small uncertainties once higher-order corrections are known, while IR contributions are included in non-perturbative matrix elements. The correspondence of the two computations is that IR part of the former computation is replaced by the matrix elements of the latter, and that the residual UV contributions of the former equals the Wilson coefficients of the latter. Thus, the uncertainties by IR renormalons in the former computation are replaced by the non-perturbative matrix elements in the latter computation. Our computation in this paper corresponds to the former type of computation. A meaningful and consistent comparison between the two types of computations would be to compare perturbative uncertainties (IR renormalons) with direct evaluation of the leading non-perturbative matrix elements. [34] We have computed the quarkonium energy levels perturbatively. In particular, the US corrections, which first appear at NNNLO, are computed perturbatively. IR contributions from the scale of order Λ QCD in these computations give rise to uncertainties (IR renormalons) of order Λ 3 QCD a 2 X . (a X denotes the typical size of the quarkonium state X.) The above estimates (iv) of perturbative uncertainties of our predictions are consistent with this estimate. Within pNRQCD, non-perturbative (IR) contributions and UV contributions can be factorized [9]. The former are given by matrix-elements of non-local gluon condensates; the latter are given by the Wilson coefficients, which can be predicted reliably by perturbative QCD, free from IR renormalons. The leading-order non-perturbative contribution is estimated to be of order Λ 3 QCD a 2 X from dimensional analysis. Thus, the renormalon uncertainty can be absorbed into a non-perturbative matrix element with the same order of magnitude. [35] The analyses in [4,7] confirm consistency of the bottomonium spectrum at NNLO with this relation. A similar feature is confirmed also for the static potential at NNLO and at NNNLO in [6]. Namely, it has been confirmed that the magnitudes of perturbative uncertainties are of order Λ 3 QCD a 2 X or Λ 3 QCD r 2 , and that the perturbative predictions are consistent with the experimental data or with lattice computations within the estimated uncertainties. Unfortunately, up to now, there exists no direct evaluation of the leading-order non-local gluon condensate by lattice simulations or by other methods. A qualitatively new aspect of our present analysis consists in the perturbative evaluation of the US corrections, which includes the perturbative evaluation of the non-local gluon condensates. The convergence of the perturbative expansions of the energy levels (within order Λ 3 QCD a 2 X uncertainties) observed at NNNLO seems to indicate that the perturbative evaluation of the gluon condensates provides reasonable order-of-magnitude estimates ∼ Λ 3 QCD a 2 X . However, a definite conclusion cannot be drawn until we know the precise value of d 3 .
Perhaps a well-known estimate of non-perturbative contributions to the energy levels of heavy quarkonium states is the Voloshin-Leutweyler formula expressed in terms of the local gluon condensate ∼ n 6 α s G a µν (0)G a µν (0) /(m 3 α 4 s ) [25]. As shown in [9], the non-local gluon condensates in pNRQCD can be ex-pressed by the local gluon condensates in the case that the time scale of US gluons T US ∼ a X /(C A α s ) is much smaller than 1/Λ QCD , namely, in the case that a X is extremely small (≪ C A α s /Λ QCD ). If, in addition, the wave functions of the quarkonium states can be approximated by the Coulomb wave functions, we obtain the non-perturbative contributions as given by the Voloshin-Leutweyler formula. Neither of these conditions, however, are met by the bottomonium states, especially by the excited states. As shown by series of studies on heavy quarkonium states in perturbative QCD, the bottomonium states lie in the intermediate-distance region, where deviation of the static potential from the Coulomb potential by the higher-order QCD corrections is essential and where the US time scale is not very much smaller than 1/Λ QCD . The Voloshin-Leutweyler formula is theoretically interesting but applicable only to hypothetical ultra-heavy quarkonium states, which lie in a deep part of the Coulomb potential. Inapplicability of the formula to the bottomonium states is signaled by an uncontrollably rapid increase (proportional to n 6 ) of the formula with the principal quantum number n. In fact, already for n ∼ 2 -4, the formula gives numerically unrealistically large contributions. This n 6 behavior results from a combination of (1) r 2 T US ∼ a 3 X behavior of the coefficient of the local gluon condensate (in contrast to r 2 ∼ a 2 X behavior of the non-local condensate) and (2) a rapid increase of the radius of the Coulomb state with n, a X ∝ n 2 , since the potential becomes flat as r increases; note that, it is the remediation of this unphysical behavior of the potential that has been essential in reproducing the gross structure of the bottomonium energy levels within perturbative QCD. Thus, such a rapid n-dependence cannot appear for the bottomonium states.
The scales µ X for the n ≥ 2 states listed in Tab. II are small and the corresponding values of α s (µ X ) are large. Hence, one may question validity of the perturbative predictions. Generally validity of a perturbative QCD prediction is examined through stability against variation of scales, convergence of perturbative series, comparison with lattice computations, and ultimately comparison with the experimental data. A common feature observed today in various (well-behaved) observables of perturbative QCD is as follows. The range of the energy scale where a prediction is stable becomes wider as we include higher-order terms of the perturbative series. The range extends not only in the UV direction but also in the IR direction. The level of stability and convergence of perturbative predictions depend on the observables and the typical scales involved in the observables. The stability range of the static potential and (consequently) that of the spectrum, after cancellation of the leading-order renormalons, have extended to surprisingly long-distance (IR) region and higher states, respectively. Concerning limitation of these perturbative predictions, we believe that the predictions are evidently invalid at r > ∼ Λ −1 MS ≈ 1 fm, where the string-breaking phenomenon takes place, and equivalently, above the BB threshold in the case of the energy levels. On the other hand, in order to judge at which r or the energy level the perturbative predictions break down before entering this nonperturbative regime, we have no other criteria than to apply the above general prescriptions to examine validity of the perturbative predicitons. In this analysis we have presented a first examination of the entire bottomonium spectrum at NNNLO.
The current status of the perturbative prediction for the bottomonium spectrum may be summarized as follows. We expect that stability of E tot (r) = 2m pole + V QCD (r) is realized, to a certain extent, as a result of decoupling of IR contributions due to a general property of the gauge theory. Nevertheless, the present status of the perturbative prediction for the bottomonium spectrum is practically determined by a fine-level cancellation between 2m pole and V QCD (r) and depends sensitively on the precise value of d 3 . If d 3 is tuned to stabilize E tot (r) optimally, we observe a reasonable agreement between the predictions and experimental data within estimated perturbative uncertainties.